Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add element center to Mesh fields and corresponding APIs and subroutines #55

Open
wants to merge 8 commits into
base: cws/pumipicDps
Choose a base branch
from
Open

Add element center to Mesh fields and corresponding APIs and subroutines #55

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
8bb7fb8
add elmCenterXZY in mesh field
RPIFisherman Mar 1, 2024
d280ff5
APIs for set/get elmCenters
RPIFisherman Mar 1, 2024
8059d6f
fix typo
RPIFisherman Mar 12, 2024
0bda51a
add elmcenter set/get test
RPIFisherman Mar 12, 2024
367c7ea
readMPAS with cell center
RPIFisherman Mar 13, 2024
92b3271
Merge branch 'cws/pumipicDps' into yg/elmCenterAPIs
RPIFisherman Mar 13, 2024
cc25d9e
Merge branch 'cws/pumipicDps' into yg/elmCenterAPIs
RPIFisherman Mar 29, 2024
370dd51
Merge branch 'cws/pumipicDps' into yg/elmCenterAPIs
RPIFisherman Mar 29, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
39 changes: 39 additions & 0 deletions src/pmpo_c.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down Expand Up @@ -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<polyMPO::MeshF_ElmCenterXYZ>();
auto h_coordsArray = Kokkos::create_mirror_view(coordsArray);
for(int i=0; i<nElements; i++){
h_coordsArray(i,0) = xArray[i];
h_coordsArray(i,1) = yArray[i];
h_coordsArray(i,2) = zArray[i];
}
Kokkos::deep_copy(coordsArray, h_coordsArray);
}

void polympo_getMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nElements, double* xArray, double* yArray, 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 device to host
auto coordsArray = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterXYZ>();
auto h_coordsArray = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),
coordsArray);
for(int i=0; i<nElements; i++){
xArray[i] = h_coordsArray(i,0);
yArray[i] = h_coordsArray(i,1);
zArray[i] = h_coordsArray(i,2);
}
}

void polympo_setMeshVel_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* uVelIn, const double* vVelIn){
//check mpMesh is valid
checkMPMeshValid(p_mpmesh);
Expand Down
2 changes: 2 additions & 0 deletions src/pmpo_c.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,8 @@ void polympo_setMeshVtxCoords_f(MPMesh_ptr p_mpmesh, const int nVertices, const
void polympo_getMeshVtxCoords_f(MPMesh_ptr p_mpmesh, const int nVertices, double* xArray, double* yArray, double* zArray);
void polympo_setMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* latitude);
void polympo_getMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, double* latitude);
void polympo_setMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nElements, const double* xArray, const double* yArray, const double* zArray);
void polympo_getMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nElements, double* xArray, double* yArray, double* zArray);
void polympo_setMeshVel_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* uVelocity, const double* vVelocity);
void polympo_getMeshVel_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uVelocity, double* vVelocity);
void polympo_setMeshOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array);//vec2d
Expand Down
26 changes: 26 additions & 0 deletions src/pmpo_fortran.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -445,6 +445,32 @@ subroutine polympo_getMeshVtxRotLat(mpMesh, nVertices, latitude) &
type(c_ptr), value :: latitude
end subroutine
!---------------------------------------------------------------------------
!> @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
Expand Down
9 changes: 9 additions & 0 deletions src/pmpo_mesh.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
12 changes: 10 additions & 2 deletions src/pmpo_mesh.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,8 @@ enum MeshFieldIndex{
MeshF_Vel,
MeshF_OnSurfVeloIncr,
MeshF_OnSurfDispIncr,
MeshF_RotLatLonIncr
MeshF_RotLatLonIncr,
MeshF_ElmCenterXYZ
};
enum MeshFieldType{
MeshFType_Invalid = -2,
Expand All @@ -42,7 +43,8 @@ const std::map<MeshFieldIndex, std::pair<MeshFieldType,
{MeshF_Vel, {MeshFType_VtxBased,"MeshField_Velocity"}},
{MeshF_OnSurfVeloIncr, {MeshFType_VtxBased,"MeshField_OnSurfaceVelocityIncrement"}},
{MeshF_OnSurfDispIncr, {MeshFType_VtxBased,"MeshField_OnSurfaceDisplacementIncrement"}},
{MeshF_RotLatLonIncr, {MeshFType_VtxBased,"MeshField_RotationalLatitudeLongitudeIncreasement"}}};
{MeshF_RotLatLonIncr, {MeshFType_VtxBased,"MeshField_RotationalLatitudeLongitudeIncreasement"}},
{MeshF_ElmCenterXYZ, {MeshFType_ElmBased,"MeshField_ElementCenterXYZ"}}};

enum mesh_type {mesh_unrecognized_lower = -1,
mesh_general_polygonal, //other meshes
Expand Down Expand Up @@ -70,6 +72,7 @@ class Mesh {
//start of meshFields
DoubleVec3dView vtxCoords_;
DoubleSclrView vtxRotLat_;
DoubleVec3dView elmCenterXYZ_;
DoubleVec2dView vtxVel_;
DoubleVec2dView vtxOnSurfVeloIncr_;
DoubleVec2dView vtxOnSurfDispIncr_;
Expand All @@ -95,6 +98,7 @@ class Mesh {
elm2ElmConn_(elm2ElmConn){
meshEdit_ = true;
setMeshVtxBasedFieldSize();
setMeshElmBasedFieldSize();
meshEdit_ = false;
vtxCoords_ = vtxCoords;
}
Expand All @@ -112,6 +116,7 @@ class Mesh {
IntElm2ElmView getElm2ElmConn() { return elm2ElmConn_; }
template<MeshFieldIndex index> 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
Expand Down Expand Up @@ -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);
}
Expand Down
59 changes: 59 additions & 0 deletions test/readMPAS.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -70,13 +72,17 @@ 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, &
nCells, nVertices, nEdgesOnCell, &
onSphere, sphereRadius, &
xVertex, yVertex, zVertex, &
latVertex, lonVertex, &
xCell, yCell, zCell, &
verticesOnCell, cellsOnCell)
use :: netcdf
use :: iso_c_binding
Expand All @@ -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)
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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'"
Expand Down Expand Up @@ -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'"
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand Down
20 changes: 20 additions & 0 deletions test/testFortran.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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()

Expand Down
6 changes: 6 additions & 0 deletions test/testFortranCreateRebuildMPs.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
Expand All @@ -311,6 +314,9 @@ program main
deallocate(xVertex)
deallocate(yVertex)
deallocate(zVertex)
deallocate(xCell)
deallocate(yCell)
deallocate(zCell)
deallocate(verticesOnCell)
deallocate(cellsOnCell)

Expand Down
Loading
Loading