From 3d8bdd63f64a31b368576dfae8c6c35722ce0165 Mon Sep 17 00:00:00 2001 From: dhyan1272 Date: Thu, 29 Aug 2024 13:52:19 -0700 Subject: [PATCH] Merge yg/elem_centre to ac/reconstrcuction --- src/pmpo_MPMesh.cpp | 22 ++++++++++- src/pmpo_c.cpp | 38 ++++++++++++++++++ src/pmpo_c.h | 2 + src/pmpo_fortran.f90 | 26 ++++++++++++ src/pmpo_mesh.cpp | 5 ++- src/pmpo_mesh.hpp | 11 +++++- test/readMPAS.f90 | 59 ++++++++++++++++++++++++++++ test/testFortran.f90 | 20 ++++++++++ test/testFortranCreateRebuildMPs.f90 | 6 +++ test/testFortranMPAdvection.f90 | 6 +++ test/testFortranMPReconstruction.f90 | 3 ++ test/testFortranReadMPAS.f90 | 6 +++ 12 files changed, 199 insertions(+), 5 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index 81c686f0..4d459438 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -117,6 +117,8 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ auto numMPs = p_MPs->getCount(); const auto vtxCoords = p_mesh->getMeshField(); + const auto elCenters = p_mesh->getMeshField(); + auto elm2VtxConn = p_mesh->getElm2VtxConn(); auto elm2ElmConn = p_mesh->getElm2ElmConn(); @@ -156,18 +158,36 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ while(true){ int numConnElms = elm2ElmConn(iElm,0); Vec3d delta = MPnew - elmCenter(iElm); - double minDistSq = delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2]; + //New delta + Vec3d center(elCenters(iElm, 0), elCenters(iElm, 1), elCenters(iElm, 2)); + delta = MPnew - center; + + double minDistSq = delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2]; int closestElm = -1; //go through all the connected elm, calc distance for(int i=1; i<=numConnElms; i++){ int elmID = elm2ElmConn(iElm,i)-1; delta = MPnew - elmCenter(elmID); + //New delta + Vec3d center(elCenters(elmID, 0), elCenters(elmID, 1), elCenters(elmID, 2)); + delta = MPnew - center; + double neighborDistSq = delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2]; if(neighborDistSq < minDistSq){ closestElm = elmID; minDistSq = neighborDistSq; } + /* + if(elm==2095){ + printf("elm %d Loop %d Adj elm %d Mins %.15e %.15e=>%.15e %.15e %.15e C new: %.15e %.15e %.15e\n", + elm, i, elmID, neighborDistSq, minDistSq, mpTgtPos(mp,0), mpTgtPos(mp,1), mpTgtPos(mp,2), + elCenters(elmID, 0), elCenters(elmID, 1), elCenters(elmID, 2)); + } + */ } + + //if(elm==2095) printf("Starting %d Ending %d %d \n", elm, iElm, closestElm); + if(closestElm<0){ MPs2Elm(mp) = iElm; break; diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index dddb16b5..a67c6d1c 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -696,6 +696,44 @@ void polympo_getMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, double } } +void polympo_setMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nElements, const double* xArray, const double* yArray, const double* zArray){ + //chech validity + checkMPMeshValid(p_mpmesh); + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + + //check the size + PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nElements); + + //copy the host array to the device + auto coordsArray = p_mesh->getMeshField(); + auto h_coordsArray = Kokkos::create_mirror_view(coordsArray); + for(int i=0; ip_mesh; + + //check the size + PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nElements); + + //copy the device to host + auto coordsArray = p_mesh->getMeshField(); + auto h_coordsArray = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), + coordsArray); + for(int i=0; i @brief set the polympo mesh elements/cells center + !> @param mpmesh(in/out) MPMesh object + !> @param nElements(in) length of array in, use for assertion + !> @param x/y/zArray(in) the 1D arrays of element centers coords + !--------------------------------------------------------------------------- + subroutine polympo_setMeshElmCenter(mpMesh, nElements, xArray, yArray, zArray) & + bind(C, NAME='polympo_setMeshElmCenter_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: nElements + type(c_ptr), intent(in), value :: xArray, yArray, zArray + end subroutine + !--------------------------------------------------------------------------- + !> @brief get the polympo mesh elements/cells center + !> @param mpmesh(in/out) MPMesh object + !> @param nElements(in) length of array in, use for assertion + !> @param x/y/zArray(in/out) the 1D arrays of element centers coords + !--------------------------------------------------------------------------- + subroutine polympo_getMeshElmCenter(mpMesh, nElements, xArray, yArray, zArray) & + bind(C, NAME='polympo_getMeshElmCenter_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: nElements + type(c_ptr), value :: xArray, yArray, zArray + end subroutine + !--------------------------------------------------------------------------- !> @brief set the vertices velocity from a host array !> @param mpmesh(in/out) MPMesh object !> @param nVertices(in) numVertices diff --git a/src/pmpo_mesh.cpp b/src/pmpo_mesh.cpp index 11f4b1a6..9f0fce50 100644 --- a/src/pmpo_mesh.cpp +++ b/src/pmpo_mesh.cpp @@ -48,9 +48,10 @@ namespace polyMPO{ auto elmMassMapEntry = meshFields2TypeAndString.at(MeshF_ElmMass); PMT_ALWAYS_ASSERT(elmMassMapEntry.first == MeshFType_ElmBased); - elmMass_ = MeshFView(elmMassMapEntry.second,numElms_); + elmMass_ = MeshFView(elmMassMapEntry.second,numElms_); + elmCenterXYZ_ = MeshFView(meshFields2TypeAndString.at(MeshF_ElmCenterXYZ).second, numElms_); } - + void Mesh::computeRotLatLonIncr(){ Kokkos::Timer timer; PMT_ALWAYS_ASSERT(geomType_ == geom_spherical_surf); diff --git a/src/pmpo_mesh.hpp b/src/pmpo_mesh.hpp index 28634928..64538e9d 100644 --- a/src/pmpo_mesh.hpp +++ b/src/pmpo_mesh.hpp @@ -24,7 +24,8 @@ enum MeshFieldIndex{ MeshF_ElmMass, MeshF_OnSurfVeloIncr, MeshF_OnSurfDispIncr, - MeshF_RotLatLonIncr + MeshF_RotLatLonIncr, + MeshF_ElmCenterXYZ }; enum MeshFieldType{ MeshFType_Invalid = -2, @@ -42,6 +43,7 @@ template <> struct meshFieldToType < MeshF_ElmMass > { using type = Ko template <> struct meshFieldToType < MeshF_OnSurfVeloIncr > { using type = Kokkos::View; }; template <> struct meshFieldToType < MeshF_OnSurfDispIncr > { using type = Kokkos::View; }; template <> struct meshFieldToType < MeshF_RotLatLonIncr > { using type = Kokkos::View; }; +template <> struct meshFieldToType < MeshF_ElmCenterXYZ > { using type = Kokkos::View; }; template using MeshFView = typename meshFieldToType::type; @@ -57,7 +59,8 @@ const std::map vtxOnSurfVeloIncr_; MeshFView vtxOnSurfDispIncr_; MeshFView vtxRotLatLonIncr_; + MeshFViewelmCenterXYZ_; //DoubleMat2DView vtxStress_; public: @@ -187,6 +191,9 @@ auto Mesh::getMeshField(){ else if constexpr (index==MeshF_RotLatLonIncr){ return vtxRotLatLonIncr_; } + else if constexpr (index==MeshF_ElmCenterXYZ){ + return elmCenterXYZ_; + } fprintf(stderr,"Mesh Field Index error!\n"); exit(1); } diff --git a/test/readMPAS.f90 b/test/readMPAS.f90 index 0e3cf2c3..ac7549d1 100644 --- a/test/readMPAS.f90 +++ b/test/readMPAS.f90 @@ -26,6 +26,7 @@ subroutine loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, & onSphere, sphereRadius, & xVertex, yVertex, zVertex, & latVertex, & + xCell, yCell, zCell, & verticesOnCell, cellsOnCell) use :: netcdf use :: iso_c_binding @@ -39,6 +40,7 @@ subroutine loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, & integer, dimension(:), pointer :: nEdgesOnCell real(kind=MPAS_RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex + real(kind=MPAS_RKIND), dimension(:), pointer :: xCell, yCell, zCell integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell call polympo_checkPrecisionForRealKind(MPAS_RKIND) @@ -70,6 +72,9 @@ subroutine loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, & !set vtxCoords which is a mesh field call polympo_setMeshVtxCoords(mpMesh,nVertices,c_loc(xVertex),c_loc(yVertex),c_loc(zVertex)) call polympo_setMeshVtxRotLat(mpMesh,nVertices,c_loc(latVertex)) + + !set mesh element center + call polympo_setMeshElmCenter(mpMesh,nCells,c_loc(xCell),c_loc(yCell),c_loc(zCell)) end subroutine subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, & @@ -77,6 +82,7 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, & onSphere, sphereRadius, & xVertex, yVertex, zVertex, & latVertex, lonVertex, & + xCell, yCell, zCell, & verticesOnCell, cellsOnCell) use :: netcdf use :: iso_c_binding @@ -91,11 +97,13 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, & integer, dimension(:), pointer :: nEdgesOnCell real(kind=MPAS_RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex + real(kind=MPAS_RKIND), dimension(:), pointer :: xCell, yCell, zCell integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell integer :: ncid, status, nCellsID, nVerticesID, maxEdgesID, vertexDegreeID, & nEdgesOnCellID, xVertexID, yVertexID, zVertexID, & latVertexID, lonVertexID, & + xCellID, yCellID, zCellID, & verticesOnCellID, cellsOnCellID status = nf90_open(path=trim(filename), mode=nf90_nowrite, ncid=ncid) @@ -166,6 +174,9 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, & allocate(zVertex(nVertices)) allocate(latVertex(nVertices)) allocate(lonVertex(nVertices)) + allocate(xCell(nCells)) + allocate(yCell(nCells)) + allocate(zCell(nCells)) allocate(nEdgesOnCell(nCells)) allocate(verticesOnCell(maxEdges, nCells)) allocate(cellsOnCell(maxEdges, nCells)) @@ -222,6 +233,27 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, & stop end if + status = nf90_inq_varid(ncid, 'xCell', xCellID) + if (status /= nf90_noerr) then + write(0, *) "readMPASMeshFromNCFile: Error on inquire varid of 'xCell'" + write(0, *) trim(nf90_strerror(status)) + stop + end if + + status = nf90_inq_varid(ncid, 'yCell', yCellID) + if (status /= nf90_noerr) then + write(0, *) "readMPASMeshFromNCFile: Error on inquire varid of 'yCell'" + write(0, *) trim(nf90_strerror(status)) + stop + end if + + status = nf90_inq_varid(ncid, 'zCell', zCellID) + if (status /= nf90_noerr) then + write(0, *) "readMPASMeshFromNCFile: Error on inquire varid of 'zCell'" + write(0, *) trim(nf90_strerror(status)) + stop + end if + status = nf90_inq_varid(ncid, 'nEdgesOnCell', nEdgesOnCellID) if (status /= nf90_noerr) then write(0, *) "readMPASMeshFromNCFile: Error on inquire varid of 'nEdgesOnCell'" @@ -278,6 +310,27 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, & stop end if + status = nf90_get_var(ncid, xCellID, xCell) + if (status /= nf90_noerr) then + write(0, *) "readMPASMeshFromNCFile: Error on get var of 'xCell'" + write(0, *) trim(nf90_strerror(status)) + stop + end if + + status = nf90_get_var(ncid, yCellID, yCell) + if (status /= nf90_noerr) then + write(0, *) "readMPASMeshFromNCFile: Error on get var of 'yCell'" + write(0, *) trim(nf90_strerror(status)) + stop + end if + + status = nf90_get_var(ncid, zCellID, zCell) + if (status /= nf90_noerr) then + write(0, *) "readMPASMeshFromNCFile: Error on get var of 'zCell'" + write(0, *) trim(nf90_strerror(status)) + stop + end if + status = nf90_get_var(ncid, nEdgesOnCellID, nEdgesOnCell) if (status /= nf90_noerr) then write(0, *) "readMPASMeshFromNCFile: Error on get var of 'nEdgesOnCell'" @@ -315,6 +368,7 @@ subroutine setWithMPASMeshByFortran(mpMesh, fileName, n) bind(C, name="setWithMP integer, dimension(:), pointer :: nEdgesOnCell real(kind=MPAS_RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex + real(kind=MPAS_RKIND), dimension(:), pointer :: xCell, yCell, zCell integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell fileNameFortran = transfer(fileName(1:n), fileNameFortran) @@ -325,18 +379,23 @@ subroutine setWithMPASMeshByFortran(mpMesh, fileName, n) bind(C, name="setWithMP onSphere, sphereRadius, & xVertex, yVertex, zVertex, & latVertex, lonVertex, & + xCell, yCell, zCell, & verticesOnCell, cellsOnCell) call loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, & nCells, nVertices, nEdgesOnCell, & onSphere, sphereRadius, & xVertex, yVertex, zVertex, & latVertex, & + xCell, yCell, zCell, & verticesOnCell, cellsOnCell) deallocate(nEdgesOnCell) deallocate(xVertex) deallocate(yVertex) deallocate(zVertex) + deallocate(xCell) + deallocate(yCell) + deallocate(zCell) deallocate(verticesOnCell) deallocate(cellsOnCell) end subroutine diff --git a/test/testFortran.f90 b/test/testFortran.f90 index d8a08a48..74316c09 100644 --- a/test/testFortran.f90 +++ b/test/testFortran.f90 @@ -168,6 +168,26 @@ program main deallocate(yArray) deallocate(zArray) + ! test elmcenter + allocate(xArray(numElms)) + allocate(yArray(numElms)) + allocate(zArray(numElms)) + xArray = value1 + yArray = value2 + zArray = value1 + value2 + call polympo_setMeshElmCenter(mpMesh, numElms, c_loc(xArray), c_loc(yArray), c_loc(zArray)) + xArray = -1 + yArray = -1 + zArray = -1 + call polympo_getMeshElmCenter(mpMesh, numElms, c_loc(xArray), c_loc(yArray), c_loc(zArray)) + call assert(all(xArray .eq. value1), "Assert xArray == value1 Failed!") + call assert(all(yArray .eq. value2), "Assert yArray == value2 Failed!") + call assert(all(zArray .eq. value1 + value2), "Assert zArray == value1 + value2 Failed!") + + deallocate(xArray) + deallocate(yArray) + deallocate(zArray) + call polympo_deleteMPMesh(mpMesh) call polympo_finalize() diff --git a/test/testFortranCreateRebuildMPs.f90 b/test/testFortranCreateRebuildMPs.f90 index aa34c47b..110359b6 100644 --- a/test/testFortranCreateRebuildMPs.f90 +++ b/test/testFortranCreateRebuildMPs.f90 @@ -254,6 +254,7 @@ program main integer, dimension(:), pointer :: nEdgesOnCell real(kind=MPAS_RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex + real(kind=MPAS_RKIND), dimension(:), pointer :: xCell, yCell, zCell integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell real(kind=MPAS_RKIND), dimension(:,:), pointer :: mpPosition @@ -279,12 +280,14 @@ program main onSphere, sphereRadius, & xVertex, yVertex, zVertex, & latVertex, lonVertex, & + xCell, yCell, zCell, & verticesOnCell, cellsOnCell) call loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, & nCells, nVertices, nEdgesOnCell, & onSphere, sphereRadius, & xVertex, yVertex, zVertex, & latVertex, & + xCell, yCell, zCell, & verticesOnCell, cellsOnCell) !check for allocation @@ -311,6 +314,9 @@ program main deallocate(xVertex) deallocate(yVertex) deallocate(zVertex) + deallocate(xCell) + deallocate(yCell) + deallocate(zCell) deallocate(verticesOnCell) deallocate(cellsOnCell) diff --git a/test/testFortranMPAdvection.f90 b/test/testFortranMPAdvection.f90 index 72909622..2e7573e6 100644 --- a/test/testFortranMPAdvection.f90 +++ b/test/testFortranMPAdvection.f90 @@ -169,6 +169,7 @@ program main integer, dimension(:), pointer :: nEdgesOnCell real(kind=MPAS_RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex + real(kind=MPAS_RKIND), dimension(:), pointer :: xCell, yCell, zCell integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell integer :: numMPs, numMPsCount, numPush integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive @@ -201,6 +202,7 @@ program main onSphere, sphereRadius, & xVertex, yVertex, zVertex, & latVertex, lonVertex, & + xCell, yCell, zCell, & verticesOnCell, cellsOnCell) if (onSphere .ne. 'YES') then write (*,*) "The mesh is not spherical!" @@ -215,6 +217,7 @@ program main onSphere, sphereRadius, & xVertex, yVertex, zVertex, & latVertex, & + xCell, yCell, zCell, & verticesOnCell, cellsOnCell) !createMPs @@ -315,6 +318,9 @@ program main deallocate(zVertex) deallocate(latVertex) deallocate(lonVertex) + deallocate(xCell) + deallocate(yCell) + deallocate(zCell) deallocate(verticesOnCell) deallocate(cellsOnCell) deallocate(mpsPerElm) diff --git a/test/testFortranMPReconstruction.f90 b/test/testFortranMPReconstruction.f90 index c1ab18ec..58312917 100644 --- a/test/testFortranMPReconstruction.f90 +++ b/test/testFortranMPReconstruction.f90 @@ -26,6 +26,7 @@ program main integer, dimension(:), pointer :: nEdgesOnCell real(kind=MPAS_RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex + real(kind=MPAS_RKIND), dimension(:), pointer :: xCell, yCell, zCell integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell integer :: numMPs, vID integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive @@ -56,6 +57,7 @@ program main onSphere, sphereRadius, & xVertex, yVertex, zVertex, & latVertex, lonVertex, & + xCell, yCell, zCell, & verticesOnCell, cellsOnCell) if (onSphere .ne. 'YES') then write (*,*) "The mesh is not spherical!" @@ -71,6 +73,7 @@ program main onSphere, sphereRadius, & xVertex, yVertex, zVertex, & latVertex, & + xCell, yCell, zCell, & verticesOnCell, cellsOnCell) nCompsDisp = 2 diff --git a/test/testFortranReadMPAS.f90 b/test/testFortranReadMPAS.f90 index f32537ab..b51c9d37 100644 --- a/test/testFortranReadMPAS.f90 +++ b/test/testFortranReadMPAS.f90 @@ -22,6 +22,7 @@ program main integer, dimension(:), pointer :: nEdgesOnCell real(kind=MPAS_RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex + real(kind=MPAS_RKIND), dimension(:), pointer :: xCell, yCell, zCell integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell call mpi_init(ierr) @@ -46,12 +47,14 @@ program main onSphere, sphereRadius, & xVertex, yVertex, zVertex, & latVertex, lonVertex, & + xCell, yCell, zCell, & verticesOnCell, cellsOnCell) call loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, & nCells, nVertices, nEdgesOnCell, & onSphere, sphereRadius, & xVertex, yVertex, zVertex, & latVertex, & + xCell, yCell, zCell, & verticesOnCell, cellsOnCell) !todo check the value using get functions. @@ -65,6 +68,9 @@ program main deallocate(xVertex) deallocate(yVertex) deallocate(zVertex) + deallocate(xCell) + deallocate(yCell) + deallocate(zCell) deallocate(verticesOnCell) deallocate(cellsOnCell)