From c1f9118cda26f4a0a085174eebe6d2bd97cc5bf2 Mon Sep 17 00:00:00 2001 From: Robert Oehmke Date: Fri, 18 Oct 2024 12:16:21 -0600 Subject: [PATCH 1/4] Connect new transpose routehandle argument down through interface levels. --- .../Field/src/ESMF_FieldRegrid.F90 | 14 ++++++++++-- .../Mesh/include/ESMCI_MeshCap.h | 5 +++-- .../Mesh/include/ESMCI_Mesh_Regrid_Glue.h | 5 +++-- src/Infrastructure/Mesh/src/ESMCI_MeshCap.C | 10 +++++---- .../Mesh/src/ESMCI_Mesh_Regrid_Glue.C | 5 +++-- .../Regrid/interface/ESMCI_Regrid_F.C | 6 +++-- src/Infrastructure/Regrid/src/ESMF_Regrid.F90 | 22 ++++++++++++++++--- 7 files changed, 50 insertions(+), 17 deletions(-) diff --git a/src/Infrastructure/Field/src/ESMF_FieldRegrid.F90 b/src/Infrastructure/Field/src/ESMF_FieldRegrid.F90 index 34846258e4..5938139998 100644 --- a/src/Infrastructure/Field/src/ESMF_FieldRegrid.F90 +++ b/src/Infrastructure/Field/src/ESMF_FieldRegrid.F90 @@ -379,6 +379,7 @@ subroutine ESMF_FieldRegridStoreNX(srcField, dstField, keywordEnforcer, & routehandle, & factorList, factorIndexList, & weights, indices, & ! DEPRECATED ARGUMENTS + transposeRoutehandle, & srcFracField, dstFracField, & dstStatusField, & unmappedDstList, & @@ -410,6 +411,7 @@ subroutine ESMF_FieldRegridStoreNX(srcField, dstField, keywordEnforcer, & integer(ESMF_KIND_I4), pointer, optional :: factorIndexList(:,:) real(ESMF_KIND_R8), pointer, optional :: weights(:) ! DEPRECATED ARG integer(ESMF_KIND_I4), pointer, optional :: indices(:,:) ! DEPRECATED ARG + type(ESMF_RouteHandle), intent(inout), optional :: transposeRoutehandle type(ESMF_Field), intent(inout), optional :: srcFracField type(ESMF_Field), intent(inout), optional :: dstFracField type(ESMF_Field), intent(inout), optional :: dstStatusField @@ -463,7 +465,10 @@ subroutine ESMF_FieldRegridStoreNX(srcField, dstField, keywordEnforcer, & ! \item[8.6.0] Added argument {\tt vectorRegrid} to enable the user to turn on vector regridding. This ! functionality treats an undistributed dimension of the input Fields as the components of a vector and ! maps it through 3D Cartesian space to give more consistent results (especially near the pole) than -! just regridding the components individually. +! just regridding the components individually. +! +! \item[8.8.0] Added argument {\tt transposeRoutehandle} to enable the user to retrieve +! a routeHandle containing the transpose of the regrid sparse matrix. ! ! \end{description} ! \end{itemize} @@ -660,8 +665,11 @@ subroutine ESMF_FieldRegridStoreNX(srcField, dstField, keywordEnforcer, & ! The {\tt factorIndexList} array is allocated by the method and the user is responsible for deallocating it. ! \item [{[weights]}] ! \apiDeprecatedArgWithReplacement{factorList} -! \item [{[indices]}] +! \item [{[indices]}] ! \apiDeprecatedArgWithReplacement{factorIndexList} +! \item [transposeRoutehandle] +! A routeHandle to the transpose of the regrid sparse matrix. The +! transposed operation goes from {\tt dstField} to {\tt srcField}. ! \item [{[srcFracField]}] ! The fraction of each source cell participating in the regridding. Only ! valid when regridmethod is {\tt ESMF\_REGRIDMETHOD\_CONSERVE} or {\tt regridmethod=ESMF\_REGRIDMETHOD\_CONSERVE\_2ND}. @@ -1207,6 +1215,7 @@ subroutine ESMF_FieldRegridStoreNX(srcField, dstField, keywordEnforcer, & unmappedaction, localIgnoreDegenerate, & srcTermProcessing, pipeLineDepth, & routehandle, tmp_indices, tmp_weights, & + transposeRoutehandle, & unmappedDstList, & localCheckFlag, & localrc) @@ -1239,6 +1248,7 @@ subroutine ESMF_FieldRegridStoreNX(srcField, dstField, keywordEnforcer, & unmappedaction, localIgnoreDegenerate, & srcTermProcessing, pipeLineDepth, & routehandle, & + transposeRoutehandle=transposeRoutehandle, & unmappedDstList=unmappedDstList, & checkFlag=localCheckFlag, & rc=localrc) diff --git a/src/Infrastructure/Mesh/include/ESMCI_MeshCap.h b/src/Infrastructure/Mesh/include/ESMCI_MeshCap.h index 5bc2a27f5f..6ca8b729c7 100644 --- a/src/Infrastructure/Mesh/include/ESMCI_MeshCap.h +++ b/src/Infrastructure/Mesh/include/ESMCI_MeshCap.h @@ -285,8 +285,9 @@ namespace ESMCI { int *extrapNumInputLevels, int *unmappedaction, int *_ignoreDegenerate, int *srcTermProcessing, int *pipelineDepth, - ESMCI::RouteHandle **rh, int *has_rh, int *has_iw, - int *nentries, ESMCI::TempWeights **tweights, + ESMCI::RouteHandle **rh, int *has_rh, + int *has_iw,int *nentries, ESMCI::TempWeights **tweights, + ESMCI::RouteHandle **trh, int *has_trh, int *has_udl, int *_num_udl, ESMCI::TempUDL **_tudl, int *has_statusArray, ESMCI::Array **statusArray, int *checkFlag, diff --git a/src/Infrastructure/Mesh/include/ESMCI_Mesh_Regrid_Glue.h b/src/Infrastructure/Mesh/include/ESMCI_Mesh_Regrid_Glue.h index e03bf900a9..201bd2b2a8 100644 --- a/src/Infrastructure/Mesh/include/ESMCI_Mesh_Regrid_Glue.h +++ b/src/Infrastructure/Mesh/include/ESMCI_Mesh_Regrid_Glue.h @@ -62,8 +62,9 @@ void ESMCI_regrid_create(Mesh **meshsrcpp, ESMCI::Array **arraysrcpp, ESMCI::Poi int *extrapNumInputLevels, int *unmappedaction, int *_ignoreDegenerate, int *srcTermProcessing, int *pipelineDepth, - ESMCI::RouteHandle **rh, int *has_rh, int *has_iw, - int *nentries, ESMCI::TempWeights **tweights, + ESMCI::RouteHandle **rh, int *has_rh, + int *has_iw, int *nentries, ESMCI::TempWeights **tweights, + ESMCI::RouteHandle **trh, int *has_trh, int *has_udl, int *_num_udl, ESMCI::TempUDL **_tudl, int *has_statusArray, ESMCI::Array **statusArray, int *checkFlag, diff --git a/src/Infrastructure/Mesh/src/ESMCI_MeshCap.C b/src/Infrastructure/Mesh/src/ESMCI_MeshCap.C index 6fadeeaa7a..a6cf4ea4b0 100644 --- a/src/Infrastructure/Mesh/src/ESMCI_MeshCap.C +++ b/src/Infrastructure/Mesh/src/ESMCI_MeshCap.C @@ -1576,8 +1576,9 @@ void MeshCap::regrid_create( int *extrapNumInputLevels, int *unmappedaction, int *_ignoreDegenerate, int *srcTermProcessing, int *pipelineDepth, - ESMCI::RouteHandle **rh, int *has_rh, int *has_iw, - int *nentries, ESMCI::TempWeights **tweights, + ESMCI::RouteHandle **rh, int *has_rh, + int *has_iw, int *nentries, ESMCI::TempWeights **tweights, + ESMCI::RouteHandle **trh, int *has_trh, int *has_udl, int *_num_udl, ESMCI::TempUDL **_tudl, int *has_statusArray, ESMCI::Array **statusArray, int *checkFlag, @@ -1668,8 +1669,9 @@ void MeshCap::regrid_create( extrapNumInputLevels, unmappedaction, _ignoreDegenerate, srcTermProcessing, pipelineDepth, - rh, has_rh, has_iw, - nentries, tweights, + rh, has_rh, + has_iw, nentries, tweights, + trh, has_trh, has_udl, _num_udl, _tudl, has_statusArray, statusArray, checkFlag, diff --git a/src/Infrastructure/Mesh/src/ESMCI_Mesh_Regrid_Glue.C b/src/Infrastructure/Mesh/src/ESMCI_Mesh_Regrid_Glue.C index 982698e777..44404b1daf 100644 --- a/src/Infrastructure/Mesh/src/ESMCI_Mesh_Regrid_Glue.C +++ b/src/Infrastructure/Mesh/src/ESMCI_Mesh_Regrid_Glue.C @@ -99,8 +99,9 @@ void ESMCI_regrid_create( int *extrapNumInputLevels, int *unmappedaction, int *_ignoreDegenerate, int *srcTermProcessing, int *pipelineDepth, - ESMCI::RouteHandle **rh, int *has_rh, int *has_iw, - int *nentries, ESMCI::TempWeights **tweights, + ESMCI::RouteHandle **rh, int *has_rh, + int *has_iw, int *nentries, ESMCI::TempWeights **tweights, + ESMCI::RouteHandle **trh, int *has_trh, int *has_udl, int *_num_udl, ESMCI::TempUDL **_tudl, int *_has_statusArray, ESMCI::Array **_statusArray, int *_checkFlag, diff --git a/src/Infrastructure/Regrid/interface/ESMCI_Regrid_F.C b/src/Infrastructure/Regrid/interface/ESMCI_Regrid_F.C index deb33b93e9..8a1de18c87 100644 --- a/src/Infrastructure/Regrid/interface/ESMCI_Regrid_F.C +++ b/src/Infrastructure/Regrid/interface/ESMCI_Regrid_F.C @@ -66,6 +66,7 @@ extern "C" void FTN_X(c_esmc_regrid_create)(MeshCap **meshsrcpp, int *srcTermProcessing, int *pipelineDepth, ESMCI::RouteHandle **rh, int *has_rh, int *has_iw, int *nentries, ESMCI::TempWeights **tweights, + ESMCI::RouteHandle **trh, int *has_trh, int *has_udl, int *_num_udl, ESMCI::TempUDL **_tudl, int *has_statusArray, ESMCI::Array **statusArray, int *checkFlag, @@ -103,8 +104,9 @@ MeshCap::regrid_create(meshsrcpp, arraysrcpp, plsrcpp, extrapNumInputLevels, unmappedaction, _ignoreDegenerate, srcTermProcessing, pipelineDepth, - rh, has_rh, has_iw, - nentries, tweights, + rh, has_rh, + has_iw, nentries, tweights, + trh, has_trh, has_udl, _num_udl, _tudl, has_statusArray, statusArray, checkFlag, diff --git a/src/Infrastructure/Regrid/src/ESMF_Regrid.F90 b/src/Infrastructure/Regrid/src/ESMF_Regrid.F90 index 3ace9bbee7..636a4166e1 100644 --- a/src/Infrastructure/Regrid/src/ESMF_Regrid.F90 +++ b/src/Infrastructure/Regrid/src/ESMF_Regrid.F90 @@ -165,6 +165,7 @@ subroutine ESMF_RegridStore(srcMesh, srcArray, srcPointList, src_pl_used, & pipelineDepth, & routehandle, & indices, weights, & + transposeRoutehandle, & unmappedDstList, & checkFlag, & rc) @@ -194,6 +195,7 @@ subroutine ESMF_RegridStore(srcMesh, srcArray, srcPointList, src_pl_used, & integer, intent(inout), optional :: srcTermProcessing integer, intent(inout), optional :: pipelineDepth type(ESMF_RouteHandle), intent(inout), optional :: routehandle + type(ESMF_RouteHandle), intent(inout), optional :: transposeRoutehandle integer(ESMF_KIND_I4), pointer, optional :: indices(:,:) real(ESMF_KIND_R8), pointer, optional :: weights(:) integer(ESMF_KIND_I4), pointer, optional :: unmappedDstList(:) @@ -239,7 +241,7 @@ subroutine ESMF_RegridStore(srcMesh, srcArray, srcPointList, src_pl_used, & ! \end{description} !EOPI integer :: localrc - integer :: has_rh, has_iw, nentries + integer :: has_rh, has_trh, has_iw, nentries type(ESMF_TempWeights) :: tweights integer :: has_udl, num_udl type(ESMF_TempUDL) :: tudl @@ -262,7 +264,9 @@ subroutine ESMF_RegridStore(srcMesh, srcArray, srcPointList, src_pl_used, & endif ! Next, we require that the user request at least something - if (.not.(present(routehandle) .or. present(indices))) then + if (.not.(present(routehandle) .or. & + present(transposeRoutehandle) .or. & + present(indices))) then localrc = ESMF_RC_ARG_BAD if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, & ESMF_CONTEXT, rcToReturn=rc)) return @@ -282,6 +286,10 @@ subroutine ESMF_RegridStore(srcMesh, srcArray, srcPointList, src_pl_used, & if (present(indices)) has_iw = 1 if (present(unmappedDstList)) has_udl = 1 + ! Record if transpose routehandle is present + has_trh = 0 + if (present(transposeRoutehandle)) has_trh = 1 + if (present(unmappedaction)) then localunmappedaction=unmappedaction else @@ -366,6 +374,7 @@ subroutine ESMF_RegridStore(srcMesh, srcArray, srcPointList, src_pl_used, & srcTermProcessing, pipelineDepth, & routehandle, has_rh, has_iw, & nentries, tweights, & + transposeRoutehandle, has_trh, & has_udl, num_udl, tudl, & has_statusArrayInt, statusArray, & checkFlagInt, & @@ -401,13 +410,20 @@ subroutine ESMF_RegridStore(srcMesh, srcArray, srcPointList, src_pl_used, & endif endif - ! Mark route handle created + ! Mark routeHandle created if (present(routeHandle)) then call ESMF_RouteHandleSetInitCreated(routeHandle, localrc) if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, & ESMF_CONTEXT, rcToReturn=rc)) return endif + ! Mark transpose routeHandle created + if (present(transposeRoutehandle)) then + call ESMF_RouteHandleSetInitCreated(transposeRoutehandle, localrc) + if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, & + ESMF_CONTEXT, rcToReturn=rc)) return + endif + rc = ESMF_SUCCESS end subroutine ESMF_RegridStore From 2f2f919d69e3d7aaa5d96e23aaf39690d538774b Mon Sep 17 00:00:00 2001 From: Robert Oehmke Date: Mon, 21 Oct 2024 12:17:16 -0600 Subject: [PATCH 2/4] Add initial version of transpose routeHandle code to regrid create. --- .../Mesh/src/ESMCI_Mesh_Regrid_Glue.C | 115 ++++++++++++++++-- 1 file changed, 105 insertions(+), 10 deletions(-) diff --git a/src/Infrastructure/Mesh/src/ESMCI_Mesh_Regrid_Glue.C b/src/Infrastructure/Mesh/src/ESMCI_Mesh_Regrid_Glue.C index 44404b1daf..f6a97675f5 100644 --- a/src/Infrastructure/Mesh/src/ESMCI_Mesh_Regrid_Glue.C +++ b/src/Infrastructure/Mesh/src/ESMCI_Mesh_Regrid_Glue.C @@ -613,6 +613,10 @@ void ESMCI_regrid_create( //TODO: also drop PointList objects here if possible to reduce Store() memory footrint #endif + + + //// Creation of routeHandle //// + #ifdef MEMLOG_on VM::logMemInfo(std::string("RegridCreate5.2")); #endif @@ -620,7 +624,7 @@ void ESMCI_regrid_create( ESMCI_REGRID_TRACE_ENTER("NativeMesh ArraySMMStore"); // Build the RouteHandle using ArraySMMStore() - if (*has_rh != 0) { + if (*has_rh) { // Set some flags enum ESMC_TypeKind_Flag tk = ESMC_TYPEKIND_R8; @@ -645,18 +649,99 @@ void ESMCI_regrid_create( ESMC_LogDefault.Write("c_esmc_regrid_create(): Returned from ArraySMMStore().", ESMC_LOGMSG_INFO); #endif + + + //// Creation of transpose routeHandle //// + #ifdef MEMLOG_on VM::logMemInfo(std::string("RegridCreate6.0")); #endif + + + // If requested, build the transpose RouteHandle using ArraySMMStore() + if (*has_trh) { - *nentries = num_entries; - // Clean up. If has_iw, then we will use the arrays to - // fill out the users pointers. These will be deleted following a copy. - if (*has_iw == 0) { - delete [] factors; - delete [] iientries; - *nentries = 0; - } else { + ESMCI_REGRID_TRACE_ENTER("NativeMesh Transpose ArraySMMStore"); + + // Allocate transpose matrix + int *transpose_iientries = new int[iientries_entry_size*num_entries]; + + // Init to beginning entries of factor index lists + int *entry=iientries; + int *t_entry=transpose_iientries; + + // Depending on size of matrix entries, loop constructing transpose matrix + if (iientries_entry_size == 2) { + + // Loop through entries + for (int i=0; i transpose_ii(transpose_iientries, 2, transpose_larg); + + // Call into Array sparse matrix multiply store to create RouteHandle + FTN_X(c_esmc_arraysmmstoreind4)(arraydstpp, arraysrcpp, trh, &tk, factors, + &num_entries, &transpose_ii, &ignoreUnmatched, + srcTermProcessing, pipelineDepth, &localrc); + if (ESMC_LogDefault.MsgFoundError(localrc, ESMCI_ERR_PASSTHRU, + ESMC_CONTEXT, NULL)) throw localrc; // bail out with exception + + // Get rid of transposed factor index list + delete [] transpose_iientries; + + ESMCI_REGRID_TRACE_EXIT("NativeMesh Transpose ArraySMMStore"); + +#ifdef PROGRESSLOG_on + ESMC_LogDefault.Write("c_esmc_regrid_create(): Returned from transpose ArraySMMStore().", ESMC_LOGMSG_INFO); +#endif + +#ifdef MEMLOG_on + VM::logMemInfo(std::string("RegridCreate7.0")); +#endif + } + + + + //// Output of weight matrix //// + + // If user has requested weights, then save them + if (*has_iw) { + // Record the number of entries + *nentries = num_entries; + // Save off the weights so the F90 caller can allocate arrays and // copy the values. if (num_entries>0) { @@ -669,11 +754,21 @@ void ESMCI_regrid_create( // Make sure copying method below takes this into account *tweights = NULL; } + + } else { //...else get rid of them + *nentries = 0; + delete [] factors; + delete [] iientries; } - // Setup structure to transfer unmappedDstList + + //// Handle output of unmapped destination list //// + + // Init info for transferring list *_num_udl=0; *_tudl=NULL; + + // If user wants umappedDstList, then setup structure to transfer it if (*has_udl) { // Get number of unmapped points int num_udl=unmappedDstList.size(); From 6bfd140dcaa4597652729efa5f1feb116f7c2867 Mon Sep 17 00:00:00 2001 From: Robert Oehmke Date: Thu, 31 Oct 2024 13:30:10 -0600 Subject: [PATCH 3/4] Simple test and some adjustments for transpose RH. --- .../Field/src/ESMF_FieldRegrid.F90 | 4 +- .../Field/tests/ESMF_FieldRegridUTest.F90 | 562 +++++++++++++++++- 2 files changed, 562 insertions(+), 4 deletions(-) diff --git a/src/Infrastructure/Field/src/ESMF_FieldRegrid.F90 b/src/Infrastructure/Field/src/ESMF_FieldRegrid.F90 index 5938139998..c40d411a6b 100644 --- a/src/Infrastructure/Field/src/ESMF_FieldRegrid.F90 +++ b/src/Infrastructure/Field/src/ESMF_FieldRegrid.F90 @@ -667,8 +667,8 @@ subroutine ESMF_FieldRegridStoreNX(srcField, dstField, keywordEnforcer, & ! \apiDeprecatedArgWithReplacement{factorList} ! \item [{[indices]}] ! \apiDeprecatedArgWithReplacement{factorIndexList} -! \item [transposeRoutehandle] -! A routeHandle to the transpose of the regrid sparse matrix. The +! \item [transposeRoutehandle] +! A routeHandle for the transpose of the regrid sparse matrix. The ! transposed operation goes from {\tt dstField} to {\tt srcField}. ! \item [{[srcFracField]}] ! The fraction of each source cell participating in the regridding. Only diff --git a/src/Infrastructure/Field/tests/ESMF_FieldRegridUTest.F90 b/src/Infrastructure/Field/tests/ESMF_FieldRegridUTest.F90 index 97007e3763..dbd2aa38ab 100644 --- a/src/Infrastructure/Field/tests/ESMF_FieldRegridUTest.F90 +++ b/src/Infrastructure/Field/tests/ESMF_FieldRegridUTest.F90 @@ -1549,6 +1549,22 @@ program ESMF_FieldRegridUTest ! return result call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE) !------------------------------------------------------------------------ + + !------------------------------------------------------------------------ + !EX_UTest + ! Test regrid matrix + write(failMsg, *) "Test unsuccessful" + write(name, *) "Test regrid transpose." + + ! initialize + rc=ESMF_SUCCESS + + ! do test + call test_regrid_transpose(rc) + + ! return result + call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE) + !------------------------------------------------------------------------ #endif #endif @@ -46924,8 +46940,550 @@ end subroutine test_sph_vec_blnr_csG_to_llG_p + subroutine test_regrid_transpose(rc) + integer, intent(out) :: rc + logical :: correct + integer :: localrc + type(ESMF_Grid) :: srcGrid + type(ESMF_Grid) :: dstGrid + type(ESMF_Field) :: srcField + type(ESMF_Field) :: dstField + type(ESMF_Field) :: tmpField + type(ESMF_Field) :: xdstField + type(ESMF_Field) :: xsrcField + type(ESMF_Array) :: dstArray + type(ESMF_Array) :: srcArray + type(ESMF_Array) :: tmpArray + type(ESMF_RouteHandle) :: routeHandle + type(ESMF_RouteHandle) :: transposeRouteHandle + type(ESMF_ArraySpec) :: arrayspec + type(ESMF_VM) :: vm + real(ESMF_KIND_R8), pointer :: farrayPtrXC(:,:) + real(ESMF_KIND_R8), pointer :: farrayPtrYC(:,:) + real(ESMF_KIND_R8), pointer :: farrayPtr1DXC(:) + real(ESMF_KIND_R8), pointer :: farrayPtr1DYC(:) + real(ESMF_KIND_R8), pointer :: farrayPtr(:,:,:), farrayPtr2(:,:) + real(ESMF_KIND_R8), pointer :: xfarrayPtr(:,:,:) + real(ESMF_KIND_R8), pointer :: tmpfarrayPtr(:,:) + integer :: clbnd(2),cubnd(2) + integer :: fclbnd(3),fcubnd(3) + integer :: i1,i2,i3, index(2) + integer :: lDE, srclocalDECount, dstlocalDECount + real(ESMF_KIND_R8) :: coord(2) + character(len=ESMF_MAXSTR) :: string + integer src_nx, src_ny, dst_nx, dst_ny + integer num_arrays + real(ESMF_KIND_R8) :: dx,dy + + real(ESMF_KIND_R8) :: src_dx, src_dy + real(ESMF_KIND_R8) :: dst_dx, dst_dy + + real(ESMF_KIND_R8) :: lon, lat, theta, phi + real(ESMF_KIND_R8), parameter :: DEG2RAD = 3.141592653589793/180.0_ESMF_KIND_R8 + + integer :: localPet, petCount - - + ! result code + integer :: finalrc + + ! init success flag + correct=.true. + + rc=ESMF_SUCCESS + + ! 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 + + + + ! Src Grid + srcGrid=ESMF_GridCreate1PeriDimUfrm(maxIndex=(/180,180/),& + minCornerCoord=(/0.0_ESMF_KIND_R8,-90.0_ESMF_KIND_R8/), & + maxCornerCoord=(/360.0_ESMF_KIND_R8,90.0_ESMF_KIND_R8/), & + regDecomp=(/1,petCount/), & + staggerLocList=(/ESMF_STAGGERLOC_CENTER/), & + rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + ! Src Field + srcField = ESMF_FieldCreate(srcGrid, typekind=ESMF_TYPEKIND_R8, & + ungriddedLBound=(/1/), ungriddedUBound=(/2/), & ! 2D vector + staggerloc=ESMF_STAGGERLOC_CENTER, name="source", rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + ! Exact src Field + xsrcField = ESMF_FieldCreate(srcGrid, typekind=ESMF_TYPEKIND_R8, & + ungriddedLBound=(/1/), ungriddedUBound=(/2/), & ! 2D vector + staggerloc=ESMF_STAGGERLOC_CENTER, name="source exact", rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + ! Get srcArray from Field + call ESMF_FieldGet(srcField, array=srcArray, rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + ! Get number of local DEs + call ESMF_GridGet(srcGrid, localDECount=srclocalDECount, rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Destination grid + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + ! Setup dest. grid + ! Make a grid that still matches up with identical points, but is + ! only the center, so that matrix is identity, but the src/dst indices aren't + ! the same, this'll let us test the trasponse where the indices will change. + dstGrid=ESMF_GridCreate1PeriDimUfrm(maxIndex=(/180,90/),& + minCornerCoord=(/0.0_ESMF_KIND_R8,-45.0_ESMF_KIND_R8/), & + maxCornerCoord=(/360.0_ESMF_KIND_R8,45.0_ESMF_KIND_R8/), & + regDecomp=(/1,petCount/), & + staggerLocList=(/ESMF_STAGGERLOC_CENTER/), & + rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + ! Create Fields + dstField = ESMF_FieldCreate(dstGrid, typekind=ESMF_TYPEKIND_R8, & + ungriddedLBound=(/1/), ungriddedUBound=(/2/), & ! 2D vector + staggerloc=ESMF_STAGGERLOC_CENTER, name="dest", rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + xdstField = ESMF_FieldCreate(dstGrid, typekind=ESMF_TYPEKIND_R8, & + ungriddedLBound=(/1/), ungriddedUBound=(/2/), & ! 2D vector + staggerloc=ESMF_STAGGERLOC_CENTER, name="xdest", rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + ! Get dstArray from Field + call ESMF_FieldGet(dstField, array=dstArray, rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + ! Get number of local DEs + call ESMF_GridGet(dstGrid, localDECount=dstlocalDECount, rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + +#if 0 + call ESMF_GridWriteVTK(srcGrid,staggerloc=ESMF_STAGGERLOC_CENTER, & + filename="srcGrid", rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif +#endif + + !! Create routehandles + call ESMF_FieldRegridStore( & + srcField, & + dstField=dstField, & + routeHandle=routeHandle, & + transposeRouteHandle=transposeRouteHandle, & + regridmethod=ESMF_REGRIDMETHOD_NEAREST_STOD, & ! Gives a nice matrix full of just 1.0 + rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + + !!!!! Forward direction fill Src and init Dst and check results !!!! + + + ! Fill src data + do lDE=0,srclocalDECount-1 + + !! get coord 1 + call ESMF_GridGetCoord(srcGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=1, & + farrayPtr=farrayPtr1DXC, rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + call ESMF_GridGetCoord(srcGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=2, & + farrayPtr=farrayPtr1DYC, rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + ! get src pointer + call ESMF_FieldGet(srcField, lDE, farrayPtr, & + computationalLBound=fclbnd, computationalUBound=fcubnd, & + rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + !! Set Field value + do i1=fclbnd(1),fcubnd(1) + + ! Get X coord from Grid + lon = farrayPtr1DXC(i1) + theta = DEG2RAD*(lon) + + do i2=fclbnd(2),fcubnd(2) + + ! Get Y coord from Grid + lat = farrayPtr1DYC(i2) + phi = DEG2RAD*(90.-lat) + + ! initialize source field to lon and lat + farrayPtr(i1,i2,1)=lon + farrayPtr(i1,i2,2)=lat + enddo + enddo + + enddo ! lDE + + ! Init Dst and set exact answers + do lDE=0,dstlocalDECount-1 + + !! get coords + call ESMF_GridGetCoord(dstGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=1, & + farrayPtr=farrayPtr1DXC,rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + call ESMF_GridGetCoord(dstGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=2, & + farrayPtr=farrayPtr1DYC, rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + call ESMF_FieldGet(dstField, lDE, farrayPtr, & + computationalLBound=fclbnd, computationalUBound=fcubnd, & + rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + call ESMF_FieldGet(xdstField, lDE, xfarrayPtr, rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + !! Set Field value + do i1=fclbnd(1),fcubnd(1) + + ! Get X coord from Grid + lon = farrayPtr1DXC(i1) + theta = DEG2RAD*(lon) + + do i2=fclbnd(2),fcubnd(2) + + ! Get Y coord from Grid + lat = farrayPtr1DYC(i2) + phi = DEG2RAD*(90.-lat) + + ! Init exact field to lon and lat + xfarrayPtr(i1,i2,1) = lon + xfarrayPtr(i1,i2,2) = lat + + ! initialize destination field + farrayPtr(i1,i2,1)=1000.0 + farrayPtr(i1,i2,2)=1000.0 + enddo + enddo + + enddo ! lDE + + ! Do regrid + call ESMF_FieldRegrid(srcField, dstField, routeHandle, rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + ! Check results + do lDE=0,dstlocalDECount-1 + + call ESMF_FieldGet(dstField, lDE, farrayPtr, & + computationalLBound=fclbnd, computationalUBound=fcubnd, & + rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + call ESMF_FieldGet(xdstField, lDE, xfarrayPtr, & + rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + !! Make sure things look ok + do i1=fclbnd(1),fcubnd(1) + do i2=fclbnd(2),fcubnd(2) + do i3=fclbnd(3),fcubnd(3) + ! if working everything should be close to exact answer + if (abs(farrayPtr(i1,i2,i3)-xfarrayPtr(i1,i2,i3)) .gt. 1.0E-10) then + ! write(*,*) i1,i2,i3," ",farrayPtr(i1,i2,i3),xfarrayPtr(i1,i2,i3) + correct=.false. + endif + enddo + enddo + enddo + + enddo ! lDE + +#if 0 + call ESMF_GridWriteVTK(srcGrid,staggerloc=ESMF_STAGGERLOC_CENTER, & + filename="srcGrid", array1=srcArray, & + rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + call ESMF_GridWriteVTK(dstGrid,staggerloc=ESMF_STAGGERLOC_CENTER, & + filename="dstGrid", array1=dstArray, & + rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif +#endif + + ! Get rid of forward routehandle + call ESMF_FieldRegridRelease(routeHandle, rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + + !!!!! Backward direction fill Dst and init Src and check results !!!! + + + ! Fill dst data + do lDE=0,dstlocalDECount-1 + + !! get coords + call ESMF_GridGetCoord(dstGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=1, & + farrayPtr=farrayPtr1DXC,rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + call ESMF_GridGetCoord(dstGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=2, & + farrayPtr=farrayPtr1DYC, rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + call ESMF_FieldGet(dstField, lDE, farrayPtr, & + computationalLBound=fclbnd, computationalUBound=fcubnd, & + rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + !! Set Field value + do i1=fclbnd(1),fcubnd(1) + + ! Get X coord from Grid + lon = farrayPtr1DXC(i1) + theta = DEG2RAD*(lon) + + do i2=fclbnd(2),fcubnd(2) + + ! Get Y coord from Grid + lat = farrayPtr1DYC(i2) + phi = DEG2RAD*(90.-lat) + + ! initialize destination field + farrayPtr(i1,i2,1)=lon + farrayPtr(i1,i2,2)=lat + enddo + enddo + + enddo ! lDE + + + ! Init Src and set exact answers + do lDE=0,srclocalDECount-1 + + !! get coord 1 + call ESMF_GridGetCoord(srcGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=1, & + farrayPtr=farrayPtr1DXC, rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + call ESMF_GridGetCoord(srcGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=2, & + farrayPtr=farrayPtr1DYC, rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + ! get src pointer + call ESMF_FieldGet(srcField, lDE, farrayPtr, & + computationalLBound=fclbnd, computationalUBound=fcubnd, & + rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + ! Get exact src pointer + call ESMF_FieldGet(xsrcField, lDE, xfarrayPtr, & + rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + !! Set Field value + do i1=fclbnd(1),fcubnd(1) + + ! Get X coord from Grid + lon = farrayPtr1DXC(i1) + theta = DEG2RAD*(lon) + + do i2=fclbnd(2),fcubnd(2) + + ! Get Y coord from Grid + lat = farrayPtr1DYC(i2) + phi = DEG2RAD*(90.-lat) + + ! initialize source field to bad value + farrayPtr(i1,i2,1)=1000.00 + farrayPtr(i1,i2,2)=1000.00 + + ! initialize exact source field to lon,lat + xfarrayPtr(i1,i2,1)=lon + xfarrayPtr(i1,i2,2)=lat + + enddo + enddo + + enddo ! lDE + + + ! Do traspose regrid + call ESMF_FieldRegrid(dstField, srcField, transposeRouteHandle, & + zeroregion=ESMF_REGION_SELECT, rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + + ! TODO: Check results + + + ! Get rid of transpose routehandle + call ESMF_FieldRegridRelease(transposeRouteHandle, rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + 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 + + + ! Free the grids + call ESMF_GridDestroy(srcGrid, rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + call ESMF_GridDestroy(dstGrid, rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + ! return answer based on correct flag + if (correct) then + rc=ESMF_SUCCESS + else + rc=ESMF_FAILURE + endif + + end subroutine test_regrid_transpose end program ESMF_FieldRegridUTest From 3e33d87d3bf562e3014a165e8c88d9fca45d10f7 Mon Sep 17 00:00:00 2001 From: Robert Oehmke Date: Mon, 4 Nov 2024 09:43:21 -0700 Subject: [PATCH 4/4] Further additions to transpose routeHandle and better test. --- .../Field/tests/ESMF_FieldRegridUTest.F90 | 45 +++++++++++++++++-- 1 file changed, 42 insertions(+), 3 deletions(-) diff --git a/src/Infrastructure/Field/tests/ESMF_FieldRegridUTest.F90 b/src/Infrastructure/Field/tests/ESMF_FieldRegridUTest.F90 index dbd2aa38ab..f03838f4b6 100644 --- a/src/Infrastructure/Field/tests/ESMF_FieldRegridUTest.F90 +++ b/src/Infrastructure/Field/tests/ESMF_FieldRegridUTest.F90 @@ -1554,7 +1554,7 @@ program ESMF_FieldRegridUTest !EX_UTest ! Test regrid matrix write(failMsg, *) "Test unsuccessful" - write(name, *) "Test regrid transpose." + write(name, *) "Test transpose regrid routeHandle." ! initialize rc=ESMF_SUCCESS @@ -47061,7 +47061,7 @@ subroutine test_regrid_transpose(rc) ! Setup dest. grid ! Make a grid that still matches up with identical points, but is ! only the center, so that matrix is identity, but the src/dst indices aren't - ! the same, this'll let us test the trasponse where the indices will change. + ! the same, this'll let us test the transponse where the indices will change. dstGrid=ESMF_GridCreate1PeriDimUfrm(maxIndex=(/180,90/),& minCornerCoord=(/0.0_ESMF_KIND_R8,-45.0_ESMF_KIND_R8/), & maxCornerCoord=(/360.0_ESMF_KIND_R8,45.0_ESMF_KIND_R8/), & @@ -47439,8 +47439,47 @@ subroutine test_regrid_transpose(rc) endif - ! TODO: Check results + + ! Check results + do lDE=0,srclocalDECount-1 + + call ESMF_FieldGet(srcField, lDE, farrayPtr, & + computationalLBound=fclbnd, computationalUBound=fcubnd, & + rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + call ESMF_FieldGet(xsrcField, lDE, xfarrayPtr, & + rc=localrc) + if (localrc /=ESMF_SUCCESS) then + rc=ESMF_FAILURE + return + endif + + !! Make sure things look ok + do i1=fclbnd(1),fcubnd(1) + do i2=fclbnd(2),fcubnd(2) + do i3=fclbnd(3),fcubnd(3) + + ! Skip values close to 1000.0 because those are the ones that + ! won't have been regridded to. + if (abs(farrayPtr(i1,i2,i3) - 1000.0) .lt. 1.0E-10) cycle + + ! if working everything should be close to exact answer + if (abs(farrayPtr(i1,i2,i3)-xfarrayPtr(i1,i2,i3)) .gt. 1.0E-10) then +! write(*,*) "T:",i1,i2,i3," ",farrayPtr(i1,i2,i3),xfarrayPtr(i1,i2,i3) + correct=.false. + endif + enddo + enddo + enddo + + enddo ! lDE + + ! Get rid of transpose routehandle call ESMF_FieldRegridRelease(transposeRouteHandle, rc=localrc)