Skip to content

Commit f58ddb6

Browse files
dhyan1272dhyan1272
committed
CleanUp_v2
1 parent 9cdbde5 commit f58ddb6

File tree

3 files changed

+40
-141
lines changed

3 files changed

+40
-141
lines changed

‎src/pmpo_MPMesh.cpp‎

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -388,7 +388,6 @@ void MPMesh::T2LTracking(Vec2dView dx){
388388
void MPMesh::reconstructSlices() {
389389
if (reconstructSlice.size() == 0) return;
390390
Kokkos::Timer timer;
391-
calcBasis(true);
392391
for (auto const& [index, reconstruct] : reconstructSlice) {
393392
if (reconstruct) reconstruct();
394393
}
@@ -415,9 +414,12 @@ void MPMesh::push_ahead(){
415414
//Latitude Longitude increment at mesh vertices and interpolate to particle position
416415
p_mesh->computeRotLatLonIncr();
417416

418-
//Interpolates latitude longitude increments and mesh velocity increments to
419-
//MP positions
420-
sphericalInterpolationDispVelIncr(*this);
417+
//Interpolates latitude longitude, mesh velocity increments to MPs
418+
bool use3DArea=false;
419+
calcBasis(use3DArea);
420+
sphericalInterpolation<MeshF_RotLatLonIncr>(*this);
421+
sphericalInterpolation<MeshF_OnSurfVeloIncr>(*this);
422+
//sphericalInterpolationDispVelIncr(*this);
421423

422424
//Push the MPs
423425
p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius(), p_mesh->getRotatedFlag());

‎src/pmpo_MPMesh.hpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -60,6 +60,8 @@ class MPMesh{
6060
void push_swap_pos();
6161
void push();
6262

63+
//Used before advection to interpolate fields from mesh to MPs
64+
//And also before reconstruction
6365
void calcBasis(bool use3DArea);
6466

6567
//Reconstruction

‎src/pmpo_wachspressBasis.hpp‎

Lines changed: 32 additions & 137 deletions
Original file line numberDiff line numberDiff line change
@@ -302,62 +302,43 @@ void getBasisByAreaGblForm_1(Vec2d MP, int numVtxs, Vec2d* vtxCoords, double* ba
302302
}
303303
*/
304304

305-
// spherical interpolation of values from mesh vertices to MPsi
305+
// spherical interpolation of values from mesh vertices to MPs
306306
template <MeshFieldIndex meshFieldIndex>
307307
void sphericalInterpolation(MPMesh& mpMesh){
308-
Kokkos::Timer timer;
309-
auto p_mesh = mpMesh.p_mesh;
310-
auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
311-
int numVtxs = p_mesh->getNumVertices();
312-
auto elm2VtxConn = p_mesh->getElm2VtxConn();
313-
314-
auto p_MPs = mpMesh.p_MPs;
315-
auto MPsPosition = p_MPs->getPositions();
316-
double radius = p_mesh->getSphereRadius();
317-
PMT_ALWAYS_ASSERT(radius >0);
318-
constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice<meshFieldIndex>;
319-
auto mpField = p_MPs->getData<mpfIndex>();
308+
Kokkos::Timer timer;
309+
310+
auto p_mesh = mpMesh.p_mesh;
311+
auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
312+
int numVtxs = p_mesh->getNumVertices();
313+
auto elm2VtxConn = p_mesh->getElm2VtxConn();
314+
double radius = p_mesh->getSphereRadius();
315+
PMT_ALWAYS_ASSERT(radius >0);
316+
317+
auto p_MPs = mpMesh.p_MPs;
318+
auto MPsPosition = p_MPs->getPositions();
319+
auto MPsBasis = p_MPs->getData<MPF_Basis_Vals>();
320+
321+
constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice<meshFieldIndex>;
322+
auto mpField = p_MPs->getData<mpfIndex>();
320323

321-
const int numEntries = mpSliceToNumEntries<mpfIndex>();
322-
//check field correspondence
323-
auto meshField = p_mesh->getMeshField<meshFieldIndex>();
324-
325-
auto interpolation = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
326-
if(mask) { //if material point is 'active'/'enabled'
327-
Vec3d position3d(MPsPosition(mp,0),MPsPosition(mp,1),MPsPosition(mp,2));
328-
// formating
329-
Vec3d v3d[maxVtxsPerElm+1];
330-
int numVtx = elm2VtxConn(elm,0);
331-
for(int i = 1; i<=numVtx; i++){
332-
v3d[i-1][0] = vtxCoords(elm2VtxConn(elm,i)-1,0);
333-
v3d[i-1][1] = vtxCoords(elm2VtxConn(elm,i)-1,1);
334-
v3d[i-1][2] = vtxCoords(elm2VtxConn(elm,i)-1,2);
335-
}
336-
v3d[numVtx][0] = vtxCoords(elm2VtxConn(elm,1)-1,0);
337-
v3d[numVtx][1] = vtxCoords(elm2VtxConn(elm,1)-1,1);
338-
v3d[numVtx][2] = vtxCoords(elm2VtxConn(elm,1)-1,2);
339-
340-
double basisByArea3d[maxVtxsPerElm] = {0.0};
341-
initArray(basisByArea3d,maxVtxsPerElm,0.0);
342-
343-
// calc basis
344-
getBasisByAreaGblFormSpherical(position3d, numVtx, v3d, radius, basisByArea3d);
345-
346-
// interpolation step
347-
for(int entry=0; entry<numEntries; entry++){
348-
double mpValue = 0.0;
349-
for(int i=1; i<= numVtx; i++){
350-
mpValue += meshField(elm2VtxConn(elm,i)-1,entry)*basisByArea3d[i-1];
351-
}
352-
mpField(mp,entry) = mpValue;
353-
}
354-
}
355-
};
356-
p_MPs->parallel_for(interpolation, "interpolation");
357-
pumipic::RecordTime("PolyMPO_sphericalInterpolation", timer.seconds());
324+
const int numEntries = mpSliceToNumEntries<mpfIndex>();
325+
auto meshField = p_mesh->getMeshField<meshFieldIndex>();
326+
327+
auto interpolation = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
328+
if(mask) { //if material point is 'active'/'enabled'
329+
int numVtx = elm2VtxConn(elm,0);
330+
for(int entry=0; entry<numEntries; entry++){
331+
double mpValue = 0.0;
332+
for(int i=1; i<= numVtx; i++)
333+
mpValue += meshField(elm2VtxConn(elm,i)-1,entry)*MPsBasis(mp,i-1);
334+
mpField(mp,entry) = mpValue;
335+
}
336+
}
337+
};
338+
p_MPs->parallel_for(interpolation, "interpolation");
339+
pumipic::RecordTime("PolyMPO_sphericalInterpolation", timer.seconds());
358340
}
359341

360-
361342
KOKKOS_INLINE_FUNCTION
362343
void compute2DplanarTriangleArea(int numVtx,
363344
const Kokkos::View<double[maxVtxsPerElm][2], Kokkos::LayoutStride, Kokkos::MemoryTraits<Kokkos::Unmanaged>>& gnom_vtx_subview,
@@ -422,7 +403,6 @@ void compute2DplanarTriangleArea(int numVtx,
422403
}
423404
}
424405

425-
426406
KOKKOS_INLINE_FUNCTION
427407
void wachpress_weights_grads_2D(int numVtx,
428408
const Kokkos::View<double[maxVtxsPerElm][2], Kokkos::LayoutStride, Kokkos::MemoryTraits<Kokkos::Unmanaged>>& gnom_vtx_subview,
@@ -515,90 +495,5 @@ void wachpress_weights_grads_2D(int numVtx,
515495
}
516496
}
517497

518-
519-
inline void sphericalInterpolationDispVelIncr(MPMesh& mpMesh){
520-
Kokkos::Timer timer;
521-
auto p_mesh = mpMesh.p_mesh;
522-
auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
523-
int numVtxs = p_mesh->getNumVertices();
524-
auto elm2VtxConn = p_mesh->getElm2VtxConn();
525-
526-
auto p_MPs = mpMesh.p_MPs;
527-
auto MPsPosition = p_MPs->getPositions();
528-
double radius = p_mesh->getSphereRadius();
529-
PMT_ALWAYS_ASSERT(radius > 0);
530-
531-
constexpr MeshFieldIndex meshFieldIndex1 = polyMPO::MeshF_RotLatLonIncr;
532-
constexpr MeshFieldIndex meshFieldIndex2 = polyMPO::MeshF_OnSurfVeloIncr;
533-
534-
auto meshField1 = p_mesh->getMeshField<meshFieldIndex1>();
535-
auto meshField2 = p_mesh->getMeshField<meshFieldIndex2>();
536-
537-
constexpr MaterialPointSlice mpfIndex1 = meshFieldIndexToMPSlice<meshFieldIndex1>;
538-
constexpr MaterialPointSlice mpfIndex2 = meshFieldIndexToMPSlice<meshFieldIndex2>;
539-
540-
const int numEntries1 = mpSliceToNumEntries<mpfIndex1>();
541-
const int numEntries2 = mpSliceToNumEntries<mpfIndex2>();
542-
543-
auto mpField1 = p_MPs->getData<mpfIndex1>();
544-
auto mpField2 = p_MPs->getData<mpfIndex2>();
545-
546-
// Field required for calculting gnomonic projection of MPs
547-
bool use3DArea = false;
548-
bool isRotated = p_mesh->getRotatedFlag();
549-
auto gnomProjVtx = p_mesh->getMeshField<polyMPO::MeshF_VtxGnomProj>();
550-
auto gnomProjElmCenter = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterGnomProj>();
551-
552-
auto interpolation = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
553-
if(mask) {
554-
555-
double basisByArea[maxVtxsPerElm] = {0.0};
556-
initArray(basisByArea, maxVtxsPerElm, 0.0);
557-
int numVtx = elm2VtxConn(elm, 0);
558-
Vec3d position3d(MPsPosition(mp, 0), MPsPosition(mp, 1), MPsPosition(mp, 2));
559-
560-
if(use3DArea){
561-
Vec3d v3d[maxVtxsPerElm + 1];
562-
for (int i = 1; i <= numVtx; i++) {
563-
v3d[i-1][0] = vtxCoords(elm2VtxConn(elm, i) - 1, 0);
564-
v3d[i-1][1] = vtxCoords(elm2VtxConn(elm, i) - 1, 1);
565-
v3d[i-1][2] = vtxCoords(elm2VtxConn(elm, i) - 1, 2);
566-
}
567-
v3d[numVtx][0] = vtxCoords(elm2VtxConn(elm,1)-1,0);
568-
v3d[numVtx][1] = vtxCoords(elm2VtxConn(elm,1)-1,1);
569-
v3d[numVtx][2] = vtxCoords(elm2VtxConn(elm,1)-1,2);
570-
getBasisByAreaGblFormSpherical(position3d, numVtx, v3d, radius, basisByArea);
571-
}
572-
else{ //if using gnomonic Projection for weights
573-
double mpProjX, mpProjY;
574-
auto gnomProjElmCenter_sub = Kokkos::subview(gnomProjElmCenter, elm, Kokkos::ALL);
575-
if(isRotated){
576-
position3d[0] = -MPsPosition(mp, 2);
577-
position3d[2] = MPsPosition(mp, 0);
578-
}
579-
computeGnomonicProjectionAtPoint(position3d, gnomProjElmCenter_sub, mpProjX, mpProjY);
580-
auto gnom_vtx_subview = Kokkos::subview(gnomProjVtx, elm, Kokkos::ALL, Kokkos::ALL);
581-
compute2DplanarTriangleArea(numVtx, gnom_vtx_subview, mpProjX, mpProjY, basisByArea);
582-
}
583-
584-
for(int entry=0; entry<numEntries1; entry++){
585-
double mpValue = 0.0;
586-
for(int i=1; i<= numVtx; i++)
587-
mpValue += meshField1(elm2VtxConn(elm,i)-1,entry)*basisByArea[i-1];
588-
mpField1(mp,entry) = mpValue;
589-
}
590-
591-
for(int entry=0; entry<numEntries2; entry++){
592-
double mpValue = 0.0;
593-
for(int i=1; i<= numVtx; i++)
594-
mpValue += meshField2(elm2VtxConn(elm,i)-1,entry)*basisByArea[i-1];
595-
mpField2(mp,entry) = mpValue;
596-
}
597-
}//mask
598-
};//lambda
599-
p_MPs->parallel_for(interpolation, "sphericalInterpolationMultiField");
600-
pumipic::RecordTime("PolyMPO_sphericalInterpolationDispVelIncr", timer.seconds());
601-
}
602-
603498
} //namespace polyMPO end
604499
#endif

0 commit comments

Comments
 (0)