Skip to content

Commit 7635543

Browse files
dhyan1272dhyan1272
committed
CleanUp for PR_v3
1 parent 5856e02 commit 7635543

File tree

6 files changed

+66
-61
lines changed

6 files changed

+66
-61
lines changed

‎src/pmpo_MPMesh.cpp‎

Lines changed: 13 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,10 @@ void MPMesh::calculateStrain(){
1818
auto gnomProjElmCenter = p_mesh->getMeshField<MeshF_ElmCenterGnomProj>();
1919
auto elm2VtxConn = p_mesh->getElm2VtxConn();
2020
auto velField = p_mesh->getMeshField<MeshF_Vel>();
21+
double radius = 1.0;
22+
if(p_mesh->getGeomType() == geom_spherical_surf)
23+
radius=p_mesh->getSphereRadius();
24+
2125
bool isRotated = p_mesh->getRotatedFlag();
2226

2327
auto setMPStrainRate = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
@@ -39,7 +43,7 @@ void MPMesh::calculateStrain(){
3943
double gradBasisByArea[2*maxVtxsPerElm] = {0.0};
4044
initArray(gradBasisByArea,maxVtxsPerElm,0.0);
4145

42-
wachpress_weights_grads_2D(numVtx, gnom_vtx_subview, mpProjX, mpProjY, basisByArea, gradBasisByArea);
46+
wachpress_weights_grads_2D(numVtx, gnom_vtx_subview, mpProjX, mpProjY, radius, basisByArea, gradBasisByArea);
4347

4448
double v11 = 0.0;
4549
double v12 = 0.0;
@@ -70,7 +74,10 @@ void MPMesh::calcBasis() {
7074

7175
auto elm2VtxConn = p_mesh->getElm2VtxConn();
7276
auto vtxCoords = p_mesh->getMeshField<MeshF_VtxCoords>();
73-
double radius = p_mesh->getSphereRadius();
77+
double radius = 1.0;
78+
if(p_mesh->getGeomType() == geom_spherical_surf)
79+
radius=p_mesh->getSphereRadius();
80+
7481
//For Gnomonic Projection
7582
auto gnomProjVtx = p_mesh->getMeshField<polyMPO::MeshF_VtxGnomProj>();
7683
auto gnomProjElmCenter = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterGnomProj>();
@@ -90,14 +97,14 @@ void MPMesh::calcBasis() {
9097
auto gnomProjElmCenter_sub = Kokkos::subview(gnomProjElmCenter, elm, Kokkos::ALL);
9198
computeGnomonicProjectionAtPoint(position3d, gnomProjElmCenter_sub, mpProjX, mpProjY);
9299
auto gnom_vtx_subview = Kokkos::subview(gnomProjVtx, elm, Kokkos::ALL, Kokkos::ALL);
93-
100+
94101
double basisByArea[maxVtxsPerElm] = {0.0};
95102
initArray(basisByArea,maxVtxsPerElm, 0.0);
96103
double gradBasisByArea[2*maxVtxsPerElm] = {0.0};
97104
initArray(gradBasisByArea,maxVtxsPerElm, 0.0);
98-
99-
wachpress_weights_grads_2D(numVtx, gnom_vtx_subview, mpProjX, mpProjY, basisByArea, gradBasisByArea);
100-
105+
106+
wachpress_weights_grads_2D(numVtx, gnom_vtx_subview, mpProjX, mpProjY, radius, basisByArea, gradBasisByArea);
107+
101108
for(int i=0; i<= numVtx; i++){
102109
MPsBasis(mp,i) = basisByArea[i];
103110
}

‎src/pmpo_MPMesh.hpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -29,15 +29,15 @@ class MPMesh{
2929
std::vector<int> numOwnersOnOtherProcs;
3030
std::vector<int> numHalosOnOtherProcs;
3131
std::vector<int>haloOwnerProcs;
32-
std::vector<std::vector<int>> haloOwnerLocalIDs;
32+
std::vector<std::vector<int>> haloOwnerLocalIDs;
3333
std::vector<std::vector<int>> ownerOwnerLocalIDs;
3434
std::vector<std::vector<int>> ownerHaloLocalIDs;
3535

3636
void startCommunication();
3737
void communicateFields(const std::vector<std::vector<double>>& fieldData, const int numEntities, const int numEntries, int mode,
3838
std::vector<std::vector<int>>& recvIDVec, std::vector<std::vector<double>>& recvDataVec);
3939
void communicate_and_take_halo_contributions(const Kokkos::View<double**>& meshField, int nEntities, int numEntries, int mode, int op);
40-
40+
4141
MPMesh(Mesh* inMesh, MaterialPoints* inMPs):
4242
p_mesh(inMesh), p_MPs(inMPs) {
4343
};

‎src/pmpo_materialPoints.hpp‎

Lines changed: 37 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -210,47 +210,47 @@ class MaterialPoints {
210210
ps::parallel_for(MPs, swap, "swap");
211211
}
212212
void updateMPSliceAll(){
213-
updateMPElmID();
214-
updateMPSlice<MPF_Cur_Pos_Rot_Lat_Lon,MPF_Tgt_Pos_Rot_Lat_Lon>();
215-
updateMPSlice<MPF_Cur_Pos_XYZ,MPF_Tgt_Pos_XYZ>();
213+
updateMPElmID();
214+
updateMPSlice<MPF_Cur_Pos_Rot_Lat_Lon,MPF_Tgt_Pos_Rot_Lat_Lon>();
215+
updateMPSlice<MPF_Cur_Pos_XYZ,MPF_Tgt_Pos_XYZ>();
216216
}
217217

218218
void updateRotLatLonAndXYZ2Tgt(const double radius, const bool isRotated){
219-
Kokkos::Timer timer;
220-
auto curPosRotLatLon = MPs->get<MPF_Cur_Pos_Rot_Lat_Lon>();
221-
auto tgtPosRotLatLon = MPs->get<MPF_Tgt_Pos_Rot_Lat_Lon>();
222-
auto tgtPosXYZ = MPs->get<MPF_Tgt_Pos_XYZ>();
223-
auto rotLatLonIncr = MPs->get<MPF_Rot_Lat_Lon_Incr>();
224-
//Velocity
225-
auto velMPs = MPs->get<MPF_Vel>();
226-
auto velIncr = MPs->get<MPF_Vel_Incr>();
219+
Kokkos::Timer timer;
220+
auto curPosRotLatLon = MPs->get<MPF_Cur_Pos_Rot_Lat_Lon>();
221+
auto tgtPosRotLatLon = MPs->get<MPF_Tgt_Pos_Rot_Lat_Lon>();
222+
auto tgtPosXYZ = MPs->get<MPF_Tgt_Pos_XYZ>();
223+
auto rotLatLonIncr = MPs->get<MPF_Rot_Lat_Lon_Incr>();
224+
//Velocity
225+
auto velMPs = MPs->get<MPF_Vel>();
226+
auto velIncr = MPs->get<MPF_Vel_Incr>();
227227

228-
auto mpAppID = MPs->get<MPF_MP_APP_ID>();
229-
230-
auto updateRotLatLon = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
231-
if(mask){
232-
auto rotLat = curPosRotLatLon(mp,0) + rotLatLonIncr(mp,0); // phi
233-
auto rotLon = curPosRotLatLon(mp,1) + rotLatLonIncr(mp,1); // lambda
234-
tgtPosRotLatLon(mp,0) = rotLat;
235-
tgtPosRotLatLon(mp,1) = rotLon;
236-
auto geoLat = rotLat;
237-
auto geoLon = rotLon;
238-
if(isRotated){
239-
auto xyz_rot = xyz_from_lat_lon(rotLat, rotLon, radius);
240-
auto xyz_geo = grid_rotation_backward(xyz_rot);
241-
lat_lon_from_xyz(geoLat, geoLon, xyz_geo, radius);
242-
}
243-
// x=cosLon cosLat, y=sinLon cosLat, z= sinLat
244-
tgtPosXYZ(mp,0) = radius * std::cos(geoLon) * std::cos(geoLat);
245-
tgtPosXYZ(mp,1) = radius * std::sin(geoLon) * std::cos(geoLat);
246-
tgtPosXYZ(mp,2) = radius * std::sin(geoLat);
247-
velMPs(mp,0) = velMPs(mp,0) + velIncr(mp,0);
248-
velMPs(mp,1) = velMPs(mp,1) + velIncr(mp,1);
249-
}
250-
};
251-
ps::parallel_for(MPs, updateRotLatLon,"updateRotationalLatitudeLongitude");
252-
pumipic::RecordTime("PolyMPO_updateRotLatLonAndXYZ2Tgt", timer.seconds());
253-
}
228+
auto mpAppID = MPs->get<MPF_MP_APP_ID>();
229+
230+
auto updateRotLatLon = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
231+
if(mask){
232+
auto rotLat = curPosRotLatLon(mp,0) + rotLatLonIncr(mp,0); // phi
233+
auto rotLon = curPosRotLatLon(mp,1) + rotLatLonIncr(mp,1); // lambda
234+
tgtPosRotLatLon(mp,0) = rotLat;
235+
tgtPosRotLatLon(mp,1) = rotLon;
236+
auto geoLat = rotLat;
237+
auto geoLon = rotLon;
238+
if(isRotated){
239+
auto xyz_rot = xyz_from_lat_lon(rotLat, rotLon, radius);
240+
auto xyz_geo = grid_rotation_backward(xyz_rot);
241+
lat_lon_from_xyz(geoLat, geoLon, xyz_geo, radius);
242+
}
243+
// x=cosLon cosLat, y=sinLon cosLat, z= sinLat
244+
tgtPosXYZ(mp,0) = radius * std::cos(geoLon) * std::cos(geoLat);
245+
tgtPosXYZ(mp,1) = radius * std::sin(geoLon) * std::cos(geoLat);
246+
tgtPosXYZ(mp,2) = radius * std::sin(geoLat);
247+
velMPs(mp,0) = velMPs(mp,0) + velIncr(mp,0);
248+
velMPs(mp,1) = velMPs(mp,1) + velIncr(mp,1);
249+
}
250+
};
251+
ps::parallel_for(MPs, updateRotLatLon,"updateRotationalLatitudeLongitude");
252+
pumipic::RecordTime("PolyMPO_updateRotLatLonAndXYZ2Tgt", timer.seconds());
253+
}
254254

255255
template <int index>
256256
auto getData() {

‎src/pmpo_wachspressBasis.hpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -46,7 +46,7 @@ void sphericalInterpolation(MPMesh& mpMesh){
4646
KOKKOS_INLINE_FUNCTION
4747
void wachpress_weights_grads_2D(int numVtx, const Kokkos::View<double[maxVtxsPerElm][2],
4848
Kokkos::LayoutStride, Kokkos::MemoryTraits<Kokkos::Unmanaged>>& gnom_vtx_subview,
49-
double mpProjX, double mpProjY, double* basis, double* grad_basis){
49+
double mpProjX, double mpProjY, double radius, double* basis, double* grad_basis){
5050

5151
double vertCoords[2][maxVtxsPerElm + 1];
5252
for (int i = 0; i < numVtx; ++i) {
@@ -119,9 +119,9 @@ void wachpress_weights_grads_2D(int numVtx, const Kokkos::View<double[maxVtxsPer
119119

120120
for (int i = 0; i < numVtx; ++i){
121121
grad_basis[i*2 + 0] = derivative[0][i] / denominator - (W[i] / (denominator * denominator)) * derivative_sum[0];
122-
grad_basis[i*2 + 0] = grad_basis[i*2 + 0] / 6371229;
122+
grad_basis[i*2 + 0] = grad_basis[i*2 + 0] / radius;
123123
grad_basis[i*2 + 1] = derivative[1][i] / denominator - (W[i] / (denominator * denominator)) * derivative_sum[1];
124-
grad_basis[i*2 + 1] = grad_basis[i*2 + 1] / 6371229;
124+
grad_basis[i*2 + 1] = grad_basis[i*2 + 1] / radius;
125125
basis[i] = W[i] / denominator;
126126
}
127127
}

‎test/testFortranMPAdvection.f90‎

Lines changed: 5 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -14,9 +14,9 @@ subroutine setProcWedges(mpMesh, nCells, comm_size, nEdgesOnCell, verticesOnCell
1414
integer, dimension(:,:), pointer :: verticesOnCell
1515
real(kind=MPAS_RKIND), dimension(:), pointer :: lonCell
1616
real(kind=MPAS_RKIND) :: normalizedLong, min, max
17-
17+
1818
allocate(owningProc(nCells))
19-
19+
2020
min = 1000
2121
max = -1000
2222
do i = 1, nCells
@@ -74,10 +74,8 @@ subroutine runAdvectionTest2(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell
7474

7575
call calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius, -numPush)
7676
call polympo_push(mpMesh)
77-
78-
end subroutine
79-
8077

78+
end subroutine
8179

8280
subroutine runReconstructionTest(mpMesh, numMPs, numPush, nCells, nVertices, mp2Elm, &
8381
latVertex, lonVertex, nEdgesOnCell, verticesOnCell, sphereRadius)
@@ -274,7 +272,7 @@ program main
274272
latVertex, &
275273
xCell, yCell, zCell, &
276274
verticesOnCell, cellsOnCell)
277-
275+
278276
call polympo_setGnomonicProjection(mpMesh)
279277

280278
call polympo_setMPICommunicator(mpMesh, mpi_comm_handle);
@@ -368,7 +366,7 @@ program main
368366
call polympo_createMPs(mpMesh,nCells,numMPs,c_loc(mpsPerElm),c_loc(mp2Elm),c_loc(isMPActive))
369367
call polympo_setMPRotLatLon(mpMesh,2,numMPs,c_loc(mpLatLon))
370368
call polympo_setMPPositions(mpMesh,3,numMPs,c_loc(mpPosition))
371-
369+
372370
if (testType == "API") then
373371
call runApiTest(mpMesh, numMPs, nVertices, nCells, numPush, mpLatLon, mpPosition, xVertex, yVertex, zVertex, latVertex)
374372
else if (testType == "MIGRATION") then

‎test/testFortranMPReconstruction.f90‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ program main
2929
real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex
3030
real(kind=MPAS_RKIND), dimension(:), pointer :: xCell, yCell, zCell
3131
real(kind=MPAS_RKIND), dimension(:), pointer :: areaTriangle
32-
32+
3333
integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell
3434
integer :: numMPs, vID
3535
integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive, globalElms
@@ -79,7 +79,7 @@ program main
7979
latVertex, &
8080
xCell, yCell, zCell, &
8181
verticesOnCell, cellsOnCell, areaTriangle)
82-
82+
8383
call polympo_setGnomonicProjection(mpMesh)
8484
nCompsDisp = 2
8585
allocate(dispIncr(nCompsDisp,nVertices))
@@ -112,7 +112,7 @@ program main
112112
mpPosition(2,i) = yCell(i)
113113
mpPosition(3,i) = zCell(i)
114114
end do
115-
115+
116116
call polympo_setElmGlobal(mpMesh, nCells, c_loc(globalElms))
117117
call polympo_createMPs(mpMesh,nCells,numMPs,c_loc(mpsPerElm),c_loc(mp2Elm),c_loc(isMPActive))
118118
call polympo_setMPICommunicator(mpMesh, mpi_comm_handle)
@@ -145,7 +145,7 @@ program main
145145
do i = 1, nVertices
146146
call assert(meshVtxMass(i) < TEST_VAL+TOLERANCE .and. meshVtxMass(i) > TEST_VAL-TOLERANCE, "Error: wrong vtx mass")
147147
end do
148-
148+
149149

150150
! Test push reconstruction
151151
do j = 1, 5
@@ -165,8 +165,8 @@ program main
165165
.and. meshElmMass(mp2Elm(i)) > TEST_VAL-TOLERANCE, "Error: wrong elm mass")
166166
end do
167167
end do
168-
169-
168+
169+
170170
call polympo_deleteMPMesh(mpMesh)
171171
call polympo_finalize()
172172

0 commit comments

Comments
 (0)