subroutine test_MOABMeshToMesh(itrp, csrv, rc)
logical, intent(out) :: itrp
logical, intent(out) :: csrv
integer, intent(out) :: rc
integer :: localrc
type(ESMF_Mesh) :: srcMesh
type(ESMF_Mesh) :: dstMesh
type(ESMF_Field) :: srcField
type(ESMF_Field) :: dstField
type(ESMF_Field) :: xdstField
type(ESMF_Field) :: srcAreaField, dstAreaField
type(ESMF_Field) :: srcFracField, dstFracField
type(ESMF_RouteHandle) :: routeHandle
type(ESMF_ArraySpec) :: arrayspec
type(ESMF_VM) :: vm
real(ESMF_KIND_R8), pointer :: srcFarrayPtr(:), dstFarrayPtr(:), xdstFarrayPtr(:)
real(ESMF_KIND_R8), pointer :: srcAreaPtr(:), dstAreaPtr(:)
real(ESMF_KIND_R8), pointer :: srcFracPtr(:), dstFracPtr(:)
integer :: clbnd(1),cubnd(1)
integer :: i1,i2,i3
real(ESMF_KIND_R8) :: x,y,z
integer :: localPet, petCount
real(ESMF_KIND_R8) :: srcmass(1), dstmass(1), srcmassg(1), dstmassg(1)
real(ESMF_KIND_R8) :: maxerror(1), minerror(1), error
real(ESMF_KIND_R8) :: maxerrorg(1), minerrorg(1), errorg
real(ESMF_KIND_R8) :: errorTot, errorTotG
integer, pointer :: nodeIds(:),nodeOwners(:)
real(ESMF_KIND_R8), pointer :: nodeCoords(:)
integer, pointer :: elemIds(:),elemTypes(:),elemConn(:),elemMask(:)
integer :: numNodes
integer :: iconn,inode
integer :: numQuadElems,numTriElems
integer :: numPentElems,numHexElems,numTotElems
integer :: numElemConn
! result code
integer :: finalrc
! Init to success
rc=ESMF_SUCCESS
itrp=.true.
csrv=.true.
! Don't do the test is MOAB isn't available
#ifdef ESMF_MOAB
! get pet info
call ESMF_VMGetGlobal(vm, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_VMGet(vm, petCount=petCount, localPet=localpet, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! If we don't have 1 or 4 PETS then exit successfully
if ((petCount .ne. 1) .and. (petCount .ne. 4)) then
rc=ESMF_SUCCESS
return
endif
! Turn on MOAB
call ESMF_MeshSetMOAB(.true., rc=localrc)
if (localrc .ne. ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
!!!! Setup source mesh !!!!
if (petCount .eq. 1) then
! Set number of nodes
numNodes=9
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/1,2,3,4,5,6,7,8,9/)
! Allocate and fill node coordinate array.
! Since this is a 2D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(2*numNodes))
nodeCoords=(/0.0,0.0, & ! node id 1
1.0,0.0, & ! node id 2
2.0,0.0, & ! node id 3
0.0, 1.0, & ! node id 4
1.0, 1.0, & ! node id 5
2.0, 1.0, & ! node id 6
0.0, 2.0, & ! node id 7
1.0, 2.0, & ! node id 8
2.0, 2.0 /) ! node id 9
! Allocate and fill the node owner array.
! Since this Mesh is all on PET 0, it's just set to all 0.
allocate(nodeOwners(numNodes))
nodeOwners=0 ! everything on PET 0
! Set the number of each type of element, plus the total number.
numQuadElems=3
numTriElems=2
numTotElems=numQuadElems+numTriElems
! Allocate and fill the element id array.
allocate(elemIds(numTotElems))
elemIds=(/1,2,3,4,5/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numTotElems))
elemTypes=(/ESMF_MESHELEMTYPE_QUAD, & ! elem id 1
ESMF_MESHELEMTYPE_TRI, & ! elem id 2
ESMF_MESHELEMTYPE_TRI, & ! elem id 3
ESMF_MESHELEMTYPE_QUAD, & ! elem id 4
ESMF_MESHELEMTYPE_QUAD/) ! elem id 5
! Allocate and fill the element connection type array.
! Note that entries in this array refer to the
! positions in the nodeIds, etc. arrays and that
! the order and number of entries for each element
! reflects that given in the Mesh options
! section for the corresponding entry
! in the elemTypes array.
allocate(elemConn(4*numQuadElems+3*numTriElems))
elemConn=(/1,2,5,4, & ! elem id 1
2,3,5, & ! elem id 2
3,6,5, & ! elem id 3
4,5,8,7, & ! elem id 4
5,6,9,8/) ! elem id 5
else if (petCount .eq. 4) then
! Setup mesh data depending on PET
if (localPET .eq. 0) then !!! This part only for PET 0
! Set number of nodes
numNodes=4
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/1,2,4,5/)
! Allocate and fill node coordinate array.
! Since this is a 2D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(2*numNodes))
nodeCoords=(/0.0, 0.0, & ! node id 1
1.0, 0.0, & ! node id 2
0.0, 1.0, & ! node id 4
1.0, 1.0 /) ! node id 5
! Allocate and fill the node owner array.
allocate(nodeOwners(numNodes))
nodeOwners=(/0, & ! node id 1
0, & ! node id 2
0, & ! node id 4
0/) ! node id 5
! Set the number of each type of element, plus the total number.
numQuadElems=1
numTriElems=0
numTotElems=numQuadElems+numTriElems
! Allocate and fill the element id array.
allocate(elemIds(numTotElems))
elemIds=(/1/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numTotElems))
elemTypes=(/ESMF_MESHELEMTYPE_QUAD/) ! elem id 1
! Allocate and fill the element connection type array.
! Note that entry are local indices
allocate(elemConn(4*numQuadElems+3*numTriElems))
elemConn=(/1,2,4,3/) ! elem id 1
else if (localPET .eq. 1) then !!! This part only for PET 1
! Set number of nodes
numNodes=4
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/2,3,5,6/)
! Allocate and fill node coordinate array.
! Since this is a 2D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(2*numNodes))
nodeCoords=(/1.0,0.0, & ! node id 2
2.0,0.0, & ! node id 3
1.0, 1.0, & ! node id 5
2.0, 1.0 /) ! node id 6
! Allocate and fill the node owner array.
allocate(nodeOwners(numNodes))
nodeOwners=(/0, & ! node id 2
1, & ! node id 3
0, & ! node id 5
1/) ! node id 6
! Set the number of each type of element, plus the total number.
numQuadElems=0
numTriElems=2
numTotElems=numQuadElems+numTriElems
! Allocate and fill the element id array.
allocate(elemIds(numTotElems))
elemIds=(/2,3/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numTotElems))
elemTypes=(/ESMF_MESHELEMTYPE_TRI, & ! elem id 2
ESMF_MESHELEMTYPE_TRI/) ! elem id 3
! Allocate and fill the element connection type array.
allocate(elemConn(4*numQuadElems+3*numTriElems))
elemConn=(/1,2,3, & ! elem id 2
2,4,3/) ! elem id 3
else if (localPET .eq. 2) then !!! This part only for PET 2
! Set number of nodes
numNodes=4
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/4,5,7,8/)
! Allocate and fill node coordinate array.
! Since this is a 2D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(2*numNodes))
nodeCoords=(/0.0,1.0, & ! node id 4
1.0,1.0, & ! node id 5
0.0,2.0, & ! node id 7
1.0,2.0 /) ! node id 8
! Allocate and fill the node owner array.
! Since this Mesh is all on PET 0, it's just set to all 0.
allocate(nodeOwners(numNodes))
nodeOwners=(/0, & ! node id 4
0, & ! node id 5
2, & ! node id 7
2/) ! node id 8
! Set the number of each type of element, plus the total number.
numQuadElems=1
numTriElems=0
numTotElems=numQuadElems+numTriElems
! Allocate and fill the element id array.
allocate(elemIds(numTotElems))
elemIds=(/4/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numTotElems))
elemTypes=(/ESMF_MESHELEMTYPE_QUAD/) ! elem id 4
! Allocate and fill the element connection type array.
allocate(elemConn(4*numQuadElems+3*numTriElems))
elemConn=(/1,2,4,3/) ! elem id 4
else if (localPET .eq. 3) then !!! This part only for PET 3
! Set number of nodes
numNodes=4
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/5,6,8,9/)
! Allocate and fill node coordinate array.
! Since this is a 2D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(2*numNodes))
nodeCoords=(/1.0,1.0, & ! node id 5
2.0,1.0, & ! node id 6
1.0,2.0, & ! node id 8
2.0,2.0 /) ! node id 9
! Allocate and fill the node owner array.
allocate(nodeOwners(numNodes))
nodeOwners=(/0, & ! node id 5
1, & ! node id 6
2, & ! node id 8
3/) ! node id 9
! Set the number of each type of element, plus the total number.
numQuadElems=1
numTriElems=0
numTotElems=numQuadElems+numTriElems
! Allocate and fill the element id array.
allocate(elemIds(numTotElems))
elemIds=(/5/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numTotElems))
elemTypes=(/ESMF_MESHELEMTYPE_QUAD/) ! elem id 5
! Allocate and fill the element connection type array.
allocate(elemConn(4*numQuadElems+3*numTriElems))
elemConn=(/1,2,4,3/) ! elem id 5
endif
endif
! Create Mesh structure in 1 step
srcMesh=ESMF_MeshCreate(parametricDim=2,spatialDim=2, &
coordSys=ESMF_COORDSYS_CART, &
nodeIds=nodeIds, nodeCoords=nodeCoords, &
nodeOwners=nodeOwners, elementIds=elemIds,&
elementTypes=elemTypes, elementConn=elemConn, &
rc=rc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Array spec for fields
call ESMF_ArraySpecSet(arrayspec, 1, ESMF_TYPEKIND_R8, rc=rc)
! Create source field
srcField = ESMF_FieldCreate(srcMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="source", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Create source area field
srcAreaField = ESMF_FieldCreate(srcMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="source_area", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Create source frac field
srcFracField = ESMF_FieldCreate(srcMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="source_frac", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Load test data into the source Field
! Should only be 1 localDE
call ESMF_FieldGet(srcField, 0, srcFarrayPtr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! set interpolated function
iconn=1
do i1=1,numTotElems
! Loop through nodes in elem
! to compute point in center
x=0.0
y=0.0
do i2=1,elemTypes(i1)
inode=2*(elemConn(iconn)-1)
x=x+nodeCoords(inode+1)
y=y+nodeCoords(inode+2)
iconn=iconn+1
enddo
x=x*(1.0/REAL(elemTypes(i1),ESMF_KIND_R8))
y=y*(1.0/REAL(elemTypes(i1),ESMF_KIND_R8))
! Set source function
srcFarrayPtr(i1) = 20.0+x+y
enddo
! deallocate node data
deallocate(nodeIds)
deallocate(nodeCoords)
deallocate(nodeOwners)
! deallocate elem data
deallocate(elemIds)
deallocate(elemTypes)
deallocate(elemConn)
!!!!!!!!!!!!!!! Setup Destination Mesh !!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Creates the following mesh on
! 1 or 4 PETs. Returns an error
! if run on other than 1 or 4 PETs
!
! Mesh Ids
!
! 2.0 7 ------- 8 -------- 9
! | | |
! | 3 | 4 |
! | | |
! 1.0 4 ------- 5 -------- 6
! | | |
! | 1 | 2 |
! | | |
! 0.0 1 ------- 2 -------- 3
!
! 0.0 1.0 2.0
!
! Node Ids at corners
! Element Ids in centers
!
!!!!!
!
! The owners for 1 PET are all Pet 0.
! The owners for 4 PETs are as follows:
!
! Mesh Owners
!
! 2.0 2 ------- 2--------- 3
! | | |
! | 2 | 3 |
! | | |
! 1.0 0 ------- 0 -------- 1
! | | |
! | 0 | 1 |
! | | |
! 0.0 0 ------- 0 -------- 1
!
! 0.0 1.0 2.0
!
! Node Owners at corners
! Element Owners in centers
!
! Setup mesh info depending on the
! number of PETs
if (petCount .eq. 1) then
! Fill in node data
numNodes=9
!! node ids
allocate(nodeIds(numNodes))
nodeIds=(/1,2,3,4,5,6,7,8,9/)
!! node Coords
allocate(nodeCoords(numNodes*2))
nodeCoords=(/0.0,0.0, &
1.0,0.0, &
2.0,0.0, &
0.0,1.0, &
1.0,1.0, &
2.0,1.0, &
0.0,2.0, &
1.0,2.0, &
2.0,2.0 /)
!! node owners
allocate(nodeOwners(numNodes))
nodeOwners=0 ! everything on proc 0
! Fill in elem data
numTotElems=4
!! elem ids
allocate(elemIds(numTotElems))
elemIds=(/1,2,3,4/)
!! elem types
allocate(elemTypes(numTotElems))
elemTypes=ESMF_MESHELEMTYPE_QUAD
!! elem conn
allocate(elemConn(numTotElems*4))
elemConn=(/1,2,5,4, &
2,3,6,5, &
4,5,8,7, &
5,6,9,8/)
else if (petCount .eq. 4) then
! Setup mesh data depending on PET
if (localPet .eq. 0) then
! Fill in node data
numNodes=4
!! node ids
allocate(nodeIds(numNodes))
nodeIds=(/1,2,4,5/)
!! node Coords
allocate(nodeCoords(numNodes*2))
nodeCoords=(/0.0,0.0, &
1.0,0.0, &
0.0,1.0, &
1.0,1.0/)
!! node owners
allocate(nodeOwners(numNodes))
nodeOwners=(/0,0,0,0/) ! everything on proc 0
! Fill in elem data
numTotElems=1
!! elem ids
allocate(elemIds(numTotElems))
elemIds=(/1/)
!! elem type
allocate(elemTypes(numTotElems))
elemTypes=ESMF_MESHELEMTYPE_QUAD
!! elem conn
allocate(elemConn(numTotElems*4))
elemConn=(/1,2,4,3/)
else if (localPet .eq. 1) then
! Fill in node data
numNodes=4
!! node ids
allocate(nodeIds(numNodes))
nodeIds=(/2,3,5,6/)
!! node Coords
allocate(nodeCoords(numNodes*2))
nodeCoords=(/1.0,0.0, &
2.0,0.0, &
1.0,1.0, &
2.0,1.0/)
!! node owners
allocate(nodeOwners(numNodes))
nodeOwners=(/0,1,0,1/)
! Fill in elem data
numTotElems=1
!! elem ids
allocate(elemIds(numTotElems))
elemIds=(/2/)
!! elem type
allocate(elemTypes(numTotElems))
elemTypes=ESMF_MESHELEMTYPE_QUAD
!! elem conn
allocate(elemConn(numTotElems*4))
elemConn=(/1,2,4,3/)
else if (localPet .eq. 2) then
! Fill in node data
numNodes=4
!! node ids
allocate(nodeIds(numNodes))
nodeIds=(/4,5,7,8/)
!! node Coords
allocate(nodeCoords(numNodes*2))
nodeCoords=(/0.0,1.0, &
1.0,1.0, &
0.0,2.0, &
1.0,2.0/)
!! node owners
allocate(nodeOwners(numNodes))
nodeOwners=(/0,0,2,2/)
! Fill in elem data
numTotElems=1
!! elem ids
allocate(elemIds(numTotElems))
elemIds=(/3/)
!! elem type
allocate(elemTypes(numTotElems))
elemTypes=ESMF_MESHELEMTYPE_QUAD
!! elem conn
allocate(elemConn(numTotElems*4))
elemConn=(/1,2,4,3/)
else
! Fill in node data
numNodes=4
!! node ids
allocate(nodeIds(numNodes))
nodeIds=(/5,6,8,9/)
!! node Coords
allocate(nodeCoords(numNodes*2))
nodeCoords=(/1.0,1.0, &
2.0,1.0, &
1.0,2.0, &
2.0,2.0/)
!! node owners
allocate(nodeOwners(numNodes))
nodeOwners=(/0,1,2,3/)
! Fill in elem data
numTotElems=1
!! elem ids
allocate(elemIds(numTotElems))
elemIds=(/4/)
!! elem type
allocate(elemTypes(numTotElems))
elemTypes=ESMF_MESHELEMTYPE_QUAD
!! elem conn
allocate(elemConn(numTotElems*4))
elemConn=(/1,2,4,3/)
endif
endif
! XMRKX
! Create Mesh structure in 1 step
dstMesh=ESMF_MeshCreate(parametricDim=2,spatialDim=2, &
coordSys=ESMF_COORDSYS_CART, &
nodeIds=nodeIds, nodeCoords=nodeCoords, &
nodeOwners=nodeOwners, elementIds=elemIds,&
elementTypes=elemTypes, elementConn=elemConn, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_MeshSetMOAB(.false., rc=localrc)
if (localrc .ne. ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Array spec
call ESMF_ArraySpecSet(arrayspec, 1, ESMF_TYPEKIND_R8, rc=rc)
! Create dest. field
dstField = ESMF_FieldCreate(dstMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="dest", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Create dest. area field
dstAreaField = ESMF_FieldCreate(dstMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="dest_area", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Create dest. frac field
dstFracField = ESMF_FieldCreate(dstMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="dest_frac", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Create exact dest. field
xdstField = ESMF_FieldCreate(dstMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="xdest", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Init destination field to 0.0
! Should only be 1 localDE
call ESMF_FieldGet(dstField, 0, dstFarrayPtr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Init destination field to 0.0
dstFarrayPtr=0.0
! Init exact destination field
! Should only be 1 localDE
call ESMF_FieldGet(xdstField, 0, xdstFarrayPtr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! set interpolated function
iconn=1
do i1=1,numTotElems
! Loop through nodes in elem
! to compute point in center
x=0.0
y=0.0
do i2=1,elemTypes(i1)
inode=2*(elemConn(iconn)-1)
x=x+nodeCoords(inode+1)
y=y+nodeCoords(inode+2)
iconn=iconn+1
enddo
x=x*(1.0/REAL(elemTypes(i1),ESMF_KIND_R8))
y=y*(1.0/REAL(elemTypes(i1),ESMF_KIND_R8))
! Set source function
xdstFarrayPtr(i1) = 20.0+x+y
enddo
! For now, Easy set interpolated function
!xdstFarrayPtr=1.0
! deallocate node data
deallocate(nodeIds)
deallocate(nodeCoords)
deallocate(nodeOwners)
! deallocate elem data
deallocate(elemIds)
deallocate(elemTypes)
deallocate(elemConn)
#if 0
call ESMF_MeshWrite(srcMesh,"srcMesh")
call ESMF_MeshWrite(dstMesh,"dstMesh")
#endif
!!! Regrid forward from the A grid to the B grid
! Regrid store
call ESMF_FieldRegridStore( &
srcField, &
dstField=dstField, &
routeHandle=routeHandle, &
regridmethod=ESMF_REGRIDMETHOD_CONSERVE, &
! COMMENT THESE OUT UNTIL THAT PART IS WORKING
! dstFracField=dstFracField, &
! srcFracField=srcFracField, &
unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, &
rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Do regrid
call ESMF_FieldRegrid(srcField, dstField, routeHandle, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldRegridRelease(routeHandle, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
#if 1
! Get the integration weights
call ESMF_FieldRegridGetArea(srcAreaField, &
rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Get the integration weights
call ESMF_FieldRegridGetArea(dstAreaField, &
rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
#endif
! Check if the values are close
minerror(1) = 100000.
maxerror(1) = 0.
error = 0.
errorTot=0.0
dstmass = 0.
! get dst Field
call ESMF_FieldGet(dstField, 0, dstFarrayPtr, computationalLBound=clbnd, &
computationalUBound=cubnd, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! get exact destination Field
call ESMF_FieldGet(xdstField, 0, xdstFarrayPtr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! get dst area Field
call ESMF_FieldGet(dstAreaField, 0, dstAreaPtr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! get frac Field
call ESMF_FieldGet(dstFracField, 0, dstFracptr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! destination grid
!! check relative error
do i1=clbnd(1),cubnd(1)
! This is WRONG, shouldn't include Frac
! dstmass = dstmass + dstFracptr(i1,i2)*dstAreaptr(i1)*fptr(i1)
! Instead do this
dstmass = dstmass + dstAreaptr(i1)*dstFarrayPtr(i1)
! If this destination cell isn't covered by a sig. amount of source, then compute error on it.
! (Note that this is what SCRIP does)
!if (dstFracptr(i1) .lt. 0.999) cycle
! write(*,*) i1,"::",dstFarrayPtr(i1),xdstFarrayPtr(i1)
if (xdstFarrayPtr(i1) .ne. 0.0) then
error=ABS(dstFarrayPtr(i1) - xdstFarrayPtr(i1))/ABS(xdstFarrayPtr(i1))
errorTot=errorTot+error
if (error > maxerror(1)) then
maxerror(1) = error
endif
if (error < minerror(1)) then
minerror(1) = error
endif
else
error=ABS(dstFarrayPtr(i1) - xdstFarrayPtr(i1))/ABS(xdstFarrayPtr(i1))
errorTot=errorTot+error
if (error > maxerror(1)) then
maxerror(1) = error
endif
if (error < minerror(1)) then
minerror(1) = error
endif
endif
enddo
srcmass(1) = 0.
! get src pointer
call ESMF_FieldGet(srcField, 0, srcFarrayPtr, computationalLBound=clbnd, &
computationalUBound=cubnd, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! get src Field
call ESMF_FieldGet(srcAreaField, 0, srcAreaptr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! get frac Field
call ESMF_FieldGet(srcFracField, 0, srcFracptr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
do i1=clbnd(1),cubnd(1)
! srcmass(1) = srcmass(1) + srcFracptr(i1)*srcAreaptr(i1)*srcFarrayPtr(i1)
srcmass(1) = srcmass(1) + srcAreaptr(i1)*srcFarrayPtr(i1)
enddo
srcmassg(1) = 0.
dstmassg(1) = 0.
call ESMF_VMAllReduce(vm, srcmass, srcmassg, 1, ESMF_REDUCE_SUM, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_VMAllReduce(vm, dstmass, dstmassg, 1, ESMF_REDUCE_SUM, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_VMAllReduce(vm, maxerror, maxerrorg, 1, ESMF_REDUCE_MAX, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_VMAllReduce(vm, minerror, minerrorg, 1, ESMF_REDUCE_MIN, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! return answer based on correct flags
csrv = .false.
if (ABS(dstmassg(1)-srcmassg(1))/srcmassg(1) < 10E-10) csrv = .true.
itrp = .false.
if (maxerrorg(1) < 10E-2) itrp = .true.
! Uncomment these calls to see some actual regrid results
if (localPet == 0) then
write(*,*) "=== MOAB Mesh ==="
write(*,*) "Conservation:"
write(*,*) "Rel Error = ", ABS(dstmassg(1)-srcmassg(1))/srcmassg(1)
write(*,*) "SRC mass = ", srcmassg(1)
write(*,*) "DST mass = ", dstmassg(1)
write(*,*) " "
write(*,*) "Interpolation:"
write(*,*) "Max Error = ", maxerrorg(1)
write(*,*) "Min Error = ", minerrorg(1)
write(*,*) "Avg Error = ", (maxerrorg(1) + minerrorg(1))/2
write(*,*) " "
endif
! Destroy the Fields
call ESMF_FieldDestroy(srcField, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldDestroy(dstField, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldDestroy(srcAreaField, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldDestroy(dstAreaField, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldDestroy(srcFracField, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldDestroy(dstFracField, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldDestroy(xdstField, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Free the meshes
call ESMF_MeshDestroy(srcMesh, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_MeshDestroy(dstMesh, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
#endif
! rc, itrp, csrv init to success above
end subroutine test_MOABMeshToMesh