<p><b>xylar@lanl.gov</b> 2012-01-18 16:23:55 -0700 (Wed, 18 Jan 2012)</p><p>BRANCH COMMIT<br>
<br>
* Added land_ice_lifev module to interface with LifeV finite element solver<br>
* LifeV calls are currently commented out<br>
<br>
* Modified land_ice_vel subroutines to call land_ice_lifev subroutines<br>
* Separated per block and full domain initializaiton in land_ice_vel (needed by LifeV)<br>
<br>
* Added a call to land_ice_vel_solve in the time integrator, for testing<br>
<br>
Code compliles, but I need an initial condition (or the code for generating one) to test the code<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice/mpas/src/core_land_ice/Makefile
===================================================================
--- branches/land_ice/mpas/src/core_land_ice/Makefile 2012-01-18 22:43:25 UTC (rev 1389)
+++ branches/land_ice/mpas/src/core_land_ice/Makefile 2012-01-18 23:23:55 UTC (rev 1390)
@@ -4,6 +4,7 @@
mpas_land_ice_test_cases.o \
mpas_land_ice_time_integration.o \
mpas_land_ice_vel.o \
+ mpas_land_ice_lifev.o \
mpas_land_ice_global_diagnostics.o
all: core_land_ice
@@ -15,8 +16,10 @@
mpas_land_ice_time_integration.o: mpas_land_ice_vel.o
-mpas_land_ice_vel.o:
+mpas_land_ice_vel.o: mpas_land_ice_lifev.o
+mpas_land_ice_lifev.o:
+
mpas_land_ice_global_diagnostics.o:
mpas_land_ice_mpas_core.o: mpas_land_ice_global_diagnostics.o mpas_land_ice_test_cases.o mpas_land_ice_time_integration.o mpas_land_ice_vel.o
Added: branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_lifev.F
===================================================================
--- branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_lifev.F (rev 0)
+++ branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_lifev.F 2012-01-18 23:23:55 UTC (rev 1390)
@@ -0,0 +1,377 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! land_ice_lifev
+!
+!> \MPAS land-ice velocity driver
+!> \author Xylar Asay-Davis
+!> \date 18 January 2012
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the routines for interfacing with LifeV
+!>
+!
+!-----------------------------------------------------------------------
+
+module land_ice_lifev
+
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_dmpar
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: land_ice_lifev_init, &
+ land_ice_lifev_finalize, &
+ land_ice_lifev_block_init, &
+ land_ice_lifev_solve
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ ! halo exchange arrays
+ integer, dimension(:), pointer :: sendCellsArray, &
+ recvCellsArray, &
+ sendVerticesArray, &
+ recvVerticesArray, &
+ sendEdgesArray, &
+ recvEdgesArray
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine land_ice_lifev_init
+!
+!> \brief Initializes velocity solver
+!> \author Xylar Asay-Davis
+!> \date 18 January 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes the ice velocity solver.
+!
+!-----------------------------------------------------------------------
+
+ subroutine land_ice_lifev_init(domain, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (domain_type), intent(inout) :: domain
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: nCells, nEdges, nVertices, nCellsSolve, nEdgesSolve, nVerticesSolve, nVertLevels
+
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnCell, verticesOnEdge
+
+ real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell
+
+ logical, dimension(:), pointer :: verticesMask
+
+ logical :: keepVertex
+
+ err = 0
+
+ ! Note: the following code assumes that blocklist contains a single block
+ ! Significant modification will be needed when multiple blocks per processor
+ ! are supported.
+
+ !extract data from domain
+ nCellsSolve = domain % blocklist % mesh % nCellsSolve
+ nEdgesSolve = domain % blocklist % mesh % nEdgesSolve
+ nVerticesSolve = domain % blocklist % mesh % nVerticesSolve
+
+ nCells = domain % blocklist % mesh % nCells
+ nEdges = domain % blocklist % mesh % nEdges
+ nVertices = domain % blocklist % mesh % nVertices
+ nVertLevels = domain % blocklist % mesh % nVertLevels
+
+ cellsOnEdge => domain % blocklist % mesh % cellsOnEdge % array
+ cellsOnVertex => domain % blocklist % mesh % cellsOnVertex % array
+ verticesOnCell => domain % blocklist % mesh % verticesOnCell % array
+ verticesOnEdge => domain % blocklist % mesh % verticesOnEdge % array
+
+ xCell => domain % blocklist % mesh % xCell % array
+ yCell => domain % blocklist % mesh % yCell % array
+ zCell => domain % blocklist % mesh % zCell % array
+
+ !allocate(verticesMask(nVertices))
+
+ !verticesMask = .true. ! all vertices are part of the domain, since thickness isn't necessarily valid here
+
+ !call lifev_init(dminfo % comm, nVerteces)
+
+ !build send and receive arrays using exchange_list
+
+ call array_from_exchange_list(domain % blocklist % parinfo % verticesToSend, sendVerticesArray)
+ call array_from_exchange_list(domain % blocklist % parinfo % verticesToRecv, recvVerticesArray)
+ call array_from_exchange_list(domain % blocklist % parinfo % cellsToSend, sendCellsArray)
+ call array_from_exchange_list(domain % blocklist % parinfo % cellsToRecv, recvCellsArray)
+ call array_from_exchange_list(domain % blocklist % parinfo % edgesToSend, sendEdgesArray)
+ call array_from_exchange_list(domain % blocklist % parinfo % edgesToRecv, recvEdgesArray)
+
+ !call lifeV and set the grid of the velocity solver
+
+ !call lifev_set_velocity_solver_grid(nCells, nEdges, nVertices, nCellsSolve, nEdgesSolve, nVerticesSolve, &
+ ! cellsOnEdge, cellsOnVertex, verticesOnCell, verticesOnEdge, verticesMask, xCell, yCell, zCell, &
+ ! sendCellsArray, recvCellsArray, sendEdgesArray, recvEdgesArray, sendVerticesArray, recvVerticesArray)
+
+ !deallocate(verticesMask)
+
+ !call lifev_initialize_L1L2_solver(nVertLevels)
+
+ !--------------------------------------------------------------------
+
+ end subroutine land_ice_lifev_init
+
+!***********************************************************************
+!
+! routine land_ice_lifev_finalize
+!
+!> \brief Initializes velocity solver
+!> \author Xylar Asay-Davis
+!> \date 18 January 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes the ice velocity solver.
+!
+!-----------------------------------------------------------------------
+
+ subroutine land_ice_lifev_finalize(domain, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (domain_type), intent(inout) :: domain
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ err = 0
+
+ deallocate(sendCellsArray, &
+ recvCellsArray, &
+ sendVerticesArray, &
+ recvVerticesArray, &
+ sendEdgesArray, &
+ recvEdgesArray)
+
+ !--------------------------------------------------------------------
+
+ end subroutine land_ice_lifev_finalize
+
+!***********************************************************************
+!
+! routine land_ice_lifev_block_init
+!
+!> \brief Initializes velocity solver
+!> \author Xylar Asay-Davis
+!> \date 18 January 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes each block of the ice velocity solver.
+!
+!-----------------------------------------------------------------------
+
+ subroutine land_ice_lifev_block_init(mesh, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ mesh !< Input: mesh information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ err = 0
+
+ !--------------------------------------------------------------------
+
+ end subroutine land_ice_lifev_block_init
+
+!***********************************************************************
+
+ subroutine land_ice_lifev_solve(mesh, state, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ mesh !< Input: mesh information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ type (state_type), intent(in) :: &
+ state !< Input/Output: state information
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+ real (kind=RKIND), dimension(:), pointer :: &
+ thck, lsrf, beta
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ unorm
+ real (kind=RKIND), dimension(:,:,:), pointer :: &
+ tracers
+
+ integer :: index_temp
+
+
+
+ err = 0
+
+ thck => state % thck % array
+ lsrf => state % lsrf % array
+ beta => state % beta % array
+ tracers => state % tracers % array
+ unorm => state % unorm % array
+
+ index_temp = state % index_temp
+
+ ! call lifev_velocity_solver_L1L2(thck, lsrf, beta, tracers(index_temp,:,:), unorm)
+
+ !--------------------------------------------------------------------
+
+ end subroutine land_ice_lifev_solve
+
+
+ ! private subroutines
+
+ subroutine array_from_exchange_list(exlist, array)
+
+ implicit none
+
+ type (exchange_list), pointer, intent(in) :: exlist
+ type (exchange_list), pointer :: listPtr
+
+ integer, dimension(:), pointer, intent(out) :: array
+ integer :: dataSize, offset, i
+
+ dataSize = 1 !in first position we will store the size of the vector
+ listPtr => exlist
+ do while (associated(listPtr))
+ dataSize = dataSize + (listPtr % nlist) + 2
+ listPtr => listPtr % next
+ end do
+
+ allocate(array(dataSize))
+
+ array(1) = dataSize;
+ offset = 2;
+ listPtr => exlist
+ do while (associated(listPtr))
+ array(offset) = listPtr % procID
+ offset = offset + 1
+ array(offset) = listPtr % nlist
+ do i=1,listPtr % nlist
+ array(i+offset) = listPtr % list(i) -1
+ end do
+ offset = offset + listPtr % nlist + 1
+ listPtr => listPtr % next
+ end do
+
+ end subroutine array_from_exchange_list
+
+!***********************************************************************
+
+end module land_ice_lifev
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Modified: branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_mpas_core.F
===================================================================
--- branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_mpas_core.F 2012-01-18 22:43:25 UTC (rev 1389)
+++ branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_mpas_core.F 2012-01-18 23:23:55 UTC (rev 1390)
@@ -19,6 +19,7 @@
use mpas_configure
use mpas_grid_types
use land_ice_test_cases
+ use land_ice_vel, only: land_ice_vel_init
implicit none
@@ -28,7 +29,9 @@
real (kind=RKIND) :: dt
type (block_type), pointer :: block
+ integer :: err
+
if (.not. config_do_restart) call setup_land_ice_test_case(domain)
!
@@ -38,6 +41,8 @@
call simulation_clock_init(domain, dt, startTimeStamp)
+ call land_ice_vel_init(domain, err)
+
block => domain % blocklist
do while (associated(block))
call mpas_init_block(block, block % mesh, dt)
@@ -115,7 +120,7 @@
use mpas_rbf_interpolation
use mpas_vector_reconstruction
use land_ice_time_integration
- use land_ice_vel, only: land_ice_vel_init
+ use land_ice_vel, only: land_ice_vel_block_init
implicit none
@@ -139,7 +144,7 @@
block % state % time_levs(1) % state % uReconstructMeridional % array &
)
- call land_ice_vel_init(mesh, err_tmp)
+ call land_ice_vel_block_init(mesh, err_tmp)
err = ior(err, err_tmp)
if(err == 1) then
@@ -346,12 +351,15 @@
subroutine mpas_core_finalize(domain)
use mpas_grid_types
+ use land_ice_vel, only: land_ice_vel_finalize
implicit none
type (domain_type), intent(inout) :: domain
integer :: ierr
+ call land_ice_vel_finalize(domain, ierr)
+
call mpas_destroy_clock(clock, ierr)
end subroutine mpas_core_finalize
Modified: branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_time_integration.F
===================================================================
--- branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_time_integration.F 2012-01-18 22:43:25 UTC (rev 1389)
+++ branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_time_integration.F 2012-01-18 23:23:55 UTC (rev 1390)
@@ -27,6 +27,8 @@
type (block_type), pointer :: block
+ integer :: err
+
! if (trim(config_time_integration) == 'RK4') then
! call land_ice_rk4(domain, dt)
! else
@@ -35,7 +37,12 @@
! stop
! end if
+ ! a temporary test that calls the velocity solver on the state at the current time
+
block => domain % blocklist
+ call land_ice_vel_solve(block % mesh, block % state % time_levs(1) % state, err)
+
+ block => domain % blocklist
do while (associated(block))
block % state % time_levs(2) % state % xtime % scalar = timeStamp
block => block % next
Modified: branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_vel.F
===================================================================
--- branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_vel.F 2012-01-18 22:43:25 UTC (rev 1389)
+++ branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_vel.F 2012-01-18 23:23:55 UTC (rev 1390)
@@ -16,6 +16,7 @@
use mpas_grid_types
use mpas_configure
+ use land_ice_lifev
implicit none
private
@@ -34,6 +35,8 @@
!--------------------------------------------------------------------
public :: land_ice_vel_init, &
+ land_ice_vel_finalize, &
+ land_ice_vel_block_init, &
land_ice_vel_solve
!--------------------------------------------------------------------
@@ -52,15 +55,15 @@
! routine land_ice_vel_init
!
!> \brief Initializes velocity solver
-!> \author William Lipscomb
-!> \date 10 January 2012
+!> \author Xylar Asay-Davis
+!> \date 18 January 2012
!> \version SVN:$Id$
-!> \details
+!> \details
!> This routine initializes the ice velocity solver.
!
!-----------------------------------------------------------------------
- subroutine land_ice_vel_init(grid, err)
+ subroutine land_ice_vel_init(domain, err)
!-----------------------------------------------------------------
!
@@ -68,8 +71,7 @@
!
!-----------------------------------------------------------------
- type (mesh_type), intent(in) :: &
- grid !< Input: grid information
+ type (domain_type), intent(inout) :: domain
!-----------------------------------------------------------------
!
@@ -93,22 +95,27 @@
err = 0
-! cellsOnEdge => grid % cellsOnEdge % array
-! edgesOnEdge => grid % edgesOnEdge % array
-! nEdgesSolve = grid % nEdgesSolve
+ call land_ice_lifev_init(domain, err)
-! do iEdge=1,grid % nEdgesSolve
-! cell1 = cellsOnEdge(1,iEdge)
-! cell2 = cellsOnEdge(2,iEdge)
-! end do
!--------------------------------------------------------------------
end subroutine land_ice_vel_init
!***********************************************************************
+!
+! routine land_ice_vel_finalize
+!
+!> \brief Initializes velocity solver
+!> \author Xylar Asay-Davis
+!> \date 18 January 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes the ice velocity solver.
+!
+!-----------------------------------------------------------------------
- subroutine land_ice_vel_solve(grid, thck, unorm, err)
+ subroutine land_ice_vel_finalize(domain, err)
!-----------------------------------------------------------------
!
@@ -116,20 +123,109 @@
!
!-----------------------------------------------------------------
+ type (domain_type), intent(inout) :: domain
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ err = 0
+
+ call land_ice_lifev_finalize(domain, err)
+
+ !--------------------------------------------------------------------
+
+ end subroutine land_ice_vel_finalize
+
+!***********************************************************************
+!
+! routine land_ice_vel_block_init
+!
+!> \brief Initializes velocity solver
+!> \author William Lipscomb
+!> \date 10 January 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes each block of the ice velocity solver.
+!
+!-----------------------------------------------------------------------
+
+ subroutine land_ice_vel_block_init(mesh, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
type (mesh_type), intent(in) :: &
- grid !< Input: grid information
+ mesh !< Input: mesh information
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- thck !< Input: Thickness at cell center
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
!-----------------------------------------------------------------
!
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ err = 0
+
+ call land_ice_lifev_block_init(mesh, err)
+
+ !--------------------------------------------------------------------
+
+ end subroutine land_ice_vel_block_init
+
+!***********************************************************************
+
+ subroutine land_ice_vel_solve(mesh, state, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ mesh !< Input: mesh information
+
+ !-----------------------------------------------------------------
+ !
! input/output variables
!
!-----------------------------------------------------------------
- real (kind=RKIND), dimension(:,:,:), intent(in) :: &
- unorm !< Input/Output: Normal velocity at cell edge
+ type (state_type), intent(in) :: &
+ state !< Input/Output: state information
!-----------------------------------------------------------------
!
@@ -147,15 +243,8 @@
err = 0
-! cellsOnEdge => grid % cellsOnEdge % array
-! edgesOnEdge => grid % edgesOnEdge % array
-! nEdgesSolve = grid % nEdgesSolve
+ call land_ice_lifev_solve(mesh, state, err)
-! do iEdge=1,grid % nEdgesSolve
-! cell1 = cellsOnEdge(1,iEdge)
-! cell2 = cellsOnEdge(2,iEdge)
-! end do
-
!--------------------------------------------------------------------
end subroutine land_ice_vel_solve
</font>
</pre>