diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index f4d57aa1..236cfa1a 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -449,6 +449,7 @@ void polympo_setMeshNumElms_f(MPMesh_ptr p_mpmesh, const int numElms){ p_mesh->setNumElms(numElms); p_mesh->setElm2VtxConn(elm2Vtx); p_mesh->setElm2ElmConn(elm2Elm); + p_mesh->setMeshElmBasedFieldSize(); } void polympo_getMeshNumElms_f(MPMesh_ptr p_mpmesh, int & numElms) { @@ -588,6 +589,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 78b3a7be..8d6d723a 100644 --- a/src/pmpo_mesh.cpp +++ b/src/pmpo_mesh.cpp @@ -38,6 +38,15 @@ namespace polyMPO{ PMT_ALWAYS_ASSERT(vtxRotLatLonIncrMapEntry.first == MeshFType_VtxBased); vtxRotLatLonIncr_ = DoubleVec2dView(vtxRotLatLonIncrMapEntry.second,numVtxs_); } + + void Mesh::setMeshElmBasedFieldSize(){ + PMT_ALWAYS_ASSERT(meshEdit_); + + auto elmCoordsMapEntry = meshFields2TypeAndString.at(MeshF_ElmCenterXYZ); + PMT_ALWAYS_ASSERT(elmCoordsMapEntry.first == MeshFType_ElmBased); + elmCenterXYZ_ = DoubleVec3dView(elmCoordsMapEntry.second, numElms_); + + } void Mesh::computeRotLatLonIncr(){ PMT_ALWAYS_ASSERT(geomType_ == geom_spherical_surf); diff --git a/src/pmpo_mesh.hpp b/src/pmpo_mesh.hpp index 3f94256a..e4c9e18e 100644 --- a/src/pmpo_mesh.hpp +++ b/src/pmpo_mesh.hpp @@ -25,7 +25,8 @@ enum MeshFieldIndex{ MeshF_Vel, MeshF_OnSurfVeloIncr, MeshF_OnSurfDispIncr, - MeshF_RotLatLonIncr + MeshF_RotLatLonIncr, + MeshF_ElmCenterXYZ }; enum MeshFieldType{ MeshFType_Invalid = -2, @@ -42,7 +43,8 @@ const std::map auto getMeshField(); void setMeshVtxBasedFieldSize(); + void setMeshElmBasedFieldSize(); void setMeshEdit(bool meshEdit) { meshEdit_ = meshEdit; } //onec MeshType/GeomType is set to valid types, we can't change them anymore @@ -161,6 +166,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 8bf0f367..d826206e 100644 --- a/test/testFortran.f90 +++ b/test/testFortran.f90 @@ -146,6 +146,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 a1b3fb88..b7fb8f8b 100644 --- a/test/testFortranMPAdvection.f90 +++ b/test/testFortranMPAdvection.f90 @@ -25,6 +25,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 integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive @@ -53,6 +54,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!" @@ -67,6 +69,7 @@ program main onSphere, sphereRadius, & xVertex, yVertex, zVertex, & latVertex, & + xCell, yCell, zCell, & verticesOnCell, cellsOnCell) nCompsDisp = 2 @@ -177,6 +180,9 @@ program main deallocate(zVertex) deallocate(latVertex) deallocate(lonVertex) + deallocate(xCell) + deallocate(yCell) + deallocate(zCell) deallocate(verticesOnCell) deallocate(cellsOnCell) deallocate(dispIncr) 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)