Skip to content

Commit 36837b3

Browse files
dhyan1272dhyan1272
committed
RotatedFlag in p_mesh
1 parent a5ea31a commit 36837b3

File tree

7 files changed

+45
-40
lines changed

7 files changed

+45
-40
lines changed

‎src/pmpo_MPMesh.cpp‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,8 @@ void MPMesh::calcBasis(bool use3DArea) {
1919
auto gnomProjVtx = p_mesh->getMeshField<polyMPO::MeshF_VtxGnomProj>();
2020
auto gnomProjElmCenter = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterGnomProj>();
2121

22-
bool isRotated = true;
22+
bool isRotated = p_mesh->getRotatedFlag();
23+
2324
auto calcbasis = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
2425
if(mask) { //if material point is 'active'/'enabled'
2526
Vec3d position3d(MPsPosition(mp,0),MPsPosition(mp,1),MPsPosition(mp,2));
@@ -351,7 +352,7 @@ void MPMesh::push_ahead(){
351352
sphericalInterpolationDispVelIncr(*this);
352353

353354
//Push the MPs
354-
p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius());
355+
p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius(), p_mesh->getRotatedFlag());
355356
pumipic::RecordTime("PolyMPO_interpolateAndPush", timer.seconds());
356357
}
357358

@@ -387,7 +388,7 @@ void MPMesh::push(){
387388

388389
sphericalInterpolation<MeshF_RotLatLonIncr>(*this);
389390

390-
p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius()); // set Tgt_XYZ
391+
p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius(), p_mesh->getRotatedFlag());
391392

392393
auto elm2Process = p_mesh->getElm2Process();
393394

‎src/pmpo_MPMesh_assembly.hpp‎

Lines changed: 15 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -281,7 +281,7 @@ void MPMesh::startCommunication(){
281281

282282
pumipic::RecordTime("Start Communication" + std::to_string(self), timer.seconds());
283283

284-
bool isRotated=false;
284+
bool isRotated = p_mesh->getRotatedFlag();
285285
p_mesh->setGnomonicProjection(isRotated);
286286

287287
/*
@@ -416,7 +416,7 @@ void MPMesh::invertMatrix(const Kokkos::View<double**>& vtxMatrices, const doubl
416416
int nVertices = p_mesh->getNumVertices();
417417
auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
418418
auto dual_triangle_area = p_mesh->getMeshField<MeshF_DualTriangleArea>();
419-
bool isRotated=false;
419+
bool isRotated = p_mesh->getRotatedFlag();
420420

421421
double eps = 1e-7;
422422
double truncateFactor = 0.05;
@@ -435,7 +435,8 @@ void MPMesh::invertMatrix(const Kokkos::View<double**>& vtxMatrices, const doubl
435435
double Y = vtxCoords(vtx, 1)/radius;
436436
double Z = vtxCoords(vtx, 2)/radius;
437437
if(isRotated){
438-
//Change X and Z
438+
X = -vtxCoords(vtx, 2)/radius;
439+
Z = vtxCoords(vtx, 0)/radius;
439440
}
440441

441442
auto cosLat = sqrt(pow(X, 2) + pow(Y, 2));
@@ -446,7 +447,9 @@ void MPMesh::invertMatrix(const Kokkos::View<double**>& vtxMatrices, const doubl
446447
Vec3d v1 = { X * invCosLat, -Y * Z * invCosLat, Y / vtx_area_sqrt };
447448
Vec3d v2 = { 0.0, cosLat, Z /vtx_area_sqrt };
448449
if(isRotated){
449-
450+
v0 = {0.0, cosLat, Z / vtx_area_sqrt};
451+
v1 = {X * invCosLat, -Y * Z * invCosLat, Y / vtx_area_sqrt};
452+
v2 = {Y * invCosLat, X * Z *invCosLat, -X / vtx_area_sqrt};
450453
}
451454
Matrix3d rotateScaleM = {v0, v1, v2};
452455

@@ -518,7 +521,7 @@ void MPMesh::invertMatrix(const Kokkos::View<double**>& vtxMatrices, const doubl
518521
VtxCoeffs(vtx, 0, 3) = temp[2];
519522

520523
//Debugging
521-
if(vtx == 33821){
524+
if(vtx == 2630){
522525
printf("Matrices %d vtx \n", vtx);
523526
printf("[ %.15e %.15e %.15e %.15e ]\n", vtxMatrices(vtx,0), vtxMatrices(vtx,1), vtxMatrices(vtx,2), vtxMatrices(vtx,3));
524527
printf("[ %.15e %.15e %.15e %.15e ]\n", vtxMatrices(vtx,1), vtxMatrices(vtx,4), vtxMatrices(vtx,5), vtxMatrices(vtx,6));
@@ -607,11 +610,15 @@ void MPMesh::reconstruct_full() {
607610
communicate_and_take_halo_contributions(meshField, numVertices, numEntries, 0, 0);
608611
pumipic::RecordTime("Communicate Field Values" + std::to_string(self), timer.seconds());
609612

613+
Kokkos::fence();
614+
assert(cudaDeviceSynchronize() == cudaSuccess);
610615
//Debug Delete
611616
Kokkos::parallel_for("printSymmetricBlock", numVertices, KOKKOS_LAMBDA(const int vtx){
612-
if (vtx == 33821) {
613-
printf("IceArea in %d: ", vtx);
614-
printf("%25.15e \n", meshField(vtx, 0));
617+
if (vtx == 2630) {
618+
printf("Field in %d: ", vtx);
619+
for (int k=0; k<numEntries; k++)
620+
printf(" %.15e ", meshField(vtx, k));
621+
printf("\n");
615622
}
616623
});
617624

‎src/pmpo_c.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -303,7 +303,7 @@ void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh,
303303
void polympo_setMPLatLonRotatedFlag_f(MPMesh_ptr p_mpmesh, const int isRotateFlag){
304304
//chech validity
305305
checkMPMeshValid(p_mpmesh);
306-
((polyMPO::MPMesh*)p_mpmesh)->p_MPs->setRotatedFlag(isRotateFlag>0);
306+
((polyMPO::MPMesh*)p_mpmesh)->p_mesh->setRotatedFlag(isRotateFlag>0);
307307

308308
}
309309

‎src/pmpo_createTestMPMesh.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -217,7 +217,7 @@ MaterialPoints* initTestMPs(Mesh* mesh, int testMPOption){
217217
}
218218
if(geomType == geom_spherical_surf){
219219
auto p_MPs = new MaterialPoints(numElms,numMPs,positions,numMPsPerElement,MPToElement);
220-
p_MPs->setRotatedFlag(false);
220+
mesh->setRotatedFlag(false);
221221
auto mpRotLatLonField = p_MPs->getData<MPF_Cur_Pos_Rot_Lat_Lon>();
222222
auto setRotLatLon = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
223223
mpRotLatLonField(mp,0) = latLonPositions(mp,0);

‎src/pmpo_materialPoints.hpp‎

Lines changed: 6 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -125,7 +125,6 @@ class MaterialPoints {
125125
PS* MPs;
126126
int elmIDoffset = -1;
127127
int maxAppID = -1;
128-
bool isRotatedFlag = false;
129128
Operating_Mode operating_mode;
130129
RebuildHelper rebuildFields;
131130
IntFunc getAppID;
@@ -215,7 +214,8 @@ class MaterialPoints {
215214
updateMPSlice<MPF_Cur_Pos_Rot_Lat_Lon,MPF_Tgt_Pos_Rot_Lat_Lon>();
216215
updateMPSlice<MPF_Cur_Pos_XYZ,MPF_Tgt_Pos_XYZ>();
217216
}
218-
void updateRotLatLonAndXYZ2Tgt(const double radius){
217+
218+
void updateRotLatLonAndXYZ2Tgt(const double radius, const bool isRotated){
219219
Kokkos::Timer timer;
220220
auto curPosRotLatLon = MPs->get<MPF_Cur_Pos_Rot_Lat_Lon>();
221221
auto tgtPosRotLatLon = MPs->get<MPF_Tgt_Pos_Rot_Lat_Lon>();
@@ -224,8 +224,9 @@ class MaterialPoints {
224224
//Velocity
225225
auto velMPs = MPs->get<MPF_Vel>();
226226
auto velIncr = MPs->get<MPF_Vel_Incr>();
227-
228-
auto is_rotated = getRotatedFlag();
227+
228+
auto mpAppID = MPs->get<MPF_MP_APP_ID>();
229+
229230
auto updateRotLatLon = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
230231
if(mask){
231232
auto rotLat = curPosRotLatLon(mp,0) + rotLatLonIncr(mp,0); // phi
@@ -234,7 +235,7 @@ class MaterialPoints {
234235
tgtPosRotLatLon(mp,1) = rotLon;
235236
auto geoLat = rotLat;
236237
auto geoLon = rotLon;
237-
if(is_rotated){
238+
if(isRotated){
238239
auto xyz_rot = xyz_from_lat_lon(rotLat, rotLon, radius);
239240
auto xyz_geo = grid_rotation_backward(xyz_rot);
240241
lat_lon_from_xyz(geoLat, geoLon, xyz_geo, radius);
@@ -243,7 +244,6 @@ class MaterialPoints {
243244
tgtPosXYZ(mp,0) = radius * std::cos(geoLon) * std::cos(geoLat);
244245
tgtPosXYZ(mp,1) = radius * std::sin(geoLon) * std::cos(geoLat);
245246
tgtPosXYZ(mp,2) = radius * std::sin(geoLat);
246-
247247
velMPs(mp,0) = velMPs(mp,0) + velIncr(mp,0);
248248
velMPs(mp,1) = velMPs(mp,1) + velIncr(mp,1);
249249
}
@@ -281,13 +281,6 @@ class MaterialPoints {
281281
PMT_ALWAYS_ASSERT(maxAppID != -1);
282282
return maxAppID;
283283
}
284-
bool getRotatedFlag() {
285-
return isRotatedFlag;
286-
}
287-
void setRotatedFlag(bool flagSet) {
288-
isRotatedFlag = flagSet;
289-
}
290-
291284
// MUTATOR
292285
template <MaterialPointSlice index> void fillData(double value);//use PS_LAMBDA fill up to 1
293286
};// End MaterialPoints

‎src/pmpo_mesh.hpp‎

Lines changed: 17 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -112,6 +112,7 @@ class Mesh {
112112
MeshFView<MeshF_VtxGnomProj> vtxGnomProj_;
113113
MeshFView<MeshF_ElmCenterGnomProj> elmCenterGnomProj_;
114114
//DoubleMat2DView vtxStress_;
115+
bool isRotatedFlag = false;
115116

116117
public:
117118
Mesh(){};
@@ -189,6 +190,13 @@ class Mesh {
189190
void setGnomonicProjection(bool isRotated);
190191

191192
void computeRotLatLonIncr();
193+
194+
bool getRotatedFlag() {
195+
return isRotatedFlag;
196+
}
197+
void setRotatedFlag(bool flagSet) {
198+
isRotatedFlag = flagSet;
199+
}
192200
};
193201

194202
template<MeshFieldIndex index>
@@ -253,19 +261,15 @@ void Mesh::fillMeshField(int size, int numEntries, double val){
253261
KOKKOS_INLINE_FUNCTION
254262
void computeGnomonicProjectionAtPoint(const Vec3d& Coord,
255263
const Kokkos::View<double[4], Kokkos::LayoutStride, Kokkos::MemoryTraits<Kokkos::Unmanaged>>& gnomProjElmCenter_sub,
256-
double& outX, double& outY){
257-
/*
258-
printf("Inline %.15e %.15e %.15e %.15e \n", gnomProjElmCenter_sub(0), gnomProjElmCenter_sub(1), gnomProjElmCenter_sub(2),
259-
gnomProjElmCenter_sub(3));
260-
*/
261-
const double iDen = 1.0 / (gnomProjElmCenter_sub(1) * gnomProjElmCenter_sub(3) * Coord[0] +
262-
gnomProjElmCenter_sub(0) * gnomProjElmCenter_sub(3) * Coord[1] +
263-
gnomProjElmCenter_sub(2) * Coord[2]);
264-
outX = iDen * (Coord[1] * gnomProjElmCenter_sub(1) -
265-
Coord[0] * gnomProjElmCenter_sub(0));
266-
outY = iDen * (Coord[2] * gnomProjElmCenter_sub(3) - Coord[1] * gnomProjElmCenter_sub(2) * gnomProjElmCenter_sub(0) -
267-
Coord[0] * gnomProjElmCenter_sub(1) * gnomProjElmCenter_sub(2));
268-
}
264+
double& outX, double& outY){
265+
const double iDen = 1.0 / (gnomProjElmCenter_sub(1) * gnomProjElmCenter_sub(3) * Coord[0] +
266+
gnomProjElmCenter_sub(0) * gnomProjElmCenter_sub(3) * Coord[1] +
267+
gnomProjElmCenter_sub(2) * Coord[2]);
268+
outX = iDen * (Coord[1] * gnomProjElmCenter_sub(1) -
269+
Coord[0] * gnomProjElmCenter_sub(0));
270+
outY = iDen * (Coord[2] * gnomProjElmCenter_sub(3) - Coord[1] * gnomProjElmCenter_sub(2) * gnomProjElmCenter_sub(0) -
271+
Coord[0] * gnomProjElmCenter_sub(1) * gnomProjElmCenter_sub(2));
272+
}
269273

270274
}
271275

‎src/pmpo_wachspressBasis.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -451,7 +451,7 @@ inline void sphericalInterpolationDispVelIncr(MPMesh& mpMesh){
451451

452452
// Field required for calculting gnomonic projection of MPs
453453
bool use3DArea = false;
454-
bool isRotated = true;
454+
bool isRotated = p_mesh->getRotatedFlag();
455455
auto gnomProjVtx = p_mesh->getMeshField<polyMPO::MeshF_VtxGnomProj>();
456456
auto gnomProjElmCenter = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterGnomProj>();
457457

0 commit comments

Comments
 (0)