Skip to content

Commit 4d9e866

Browse files
AngelyrRPIFishermandhyan1272onkarsahni
authored
Reconstruction and Migration (#68)
* migration test * copied old code * owning proc * initialize material points * migration test working * moved out divide procs * divide in wedges * rename procWedges * use mpi_allreduce * do while * update comment * migration timer * only run test if mpi * change MPF_Constv_Mdl_Param to MeshF_Invalid * skeleton of reconstruction APIs * setReconstructSlice * assembly * more reconstruct slice types * API changes * constexpr mpSliceToMeshField * mpSliceToMeshFieldSize * switch statement * fixed spacing * renaming * moved assembly to MPMesh * moved reconstruct slice * move reconstruct down * template type conversion * reconsOption * reconstruct error * rename pmpo_assembly.hpp * seperate template index size * fixed comments * slice to type * combine size and type * fixed sizes * MPTypeHelper * compile time getSize * removed type alias * using typedef * variable renaming, notation fix * using util types * reuse vec2d * first scratch * remove mesh types * add two options for recons * add reconstruct API * test Recons * removed extra types * reuse doubleview * rename mpSliceToNumEntries * remove setReconstructionOption, recons APIs name changes * swap mesh field and material point field * add vtx/elm Mass mesh fields * set/get APIs for vtx/elm Mass mesh fields * order change, put vtxMass and elmMass together * mesh ENt TYpe * doubleScrl * assembly order * test for set/get vtx/elm Mass APIs * change order of MeshFieldIndex enums * split into two functions * fixed getmeshenttype * added elm or vtx to mesh fields * assembly refactor * elm based assembly * fortran get mesh type * removed basis and mass weight * rename mesh field type * moved mesh field type fortran * moved mpMass mpVel * test vtx mass * error msg * average num added * rename num components * 1D num components * test elm mass * average order 0 * removed for testing * calculate basis * multiply weights * removed recon option * assembly refactor * fixed set elm mass * calculate basis * remove unused * fix recon test bug * tests passing * fixed unit test * test elm push reconstruction * small reconstruct test fixes * moved tolerance * moved calculate surface displacement * moved calculate displacement * mps per edge * remove inbound * find new mp position * rename mps scale factor per vtx * count num mps starting at 0 * fixed num mps * fixed num mps * removed unused * fixed radius * num mps count * numMPs count * fixed mpPosition error * calcSurfDispIncr * bring back xc * record time * summrize time * more timers * print scale factor num mps * skip reconstruct if size == 0 * read input from command line * fixed perlmutter error * fixed elm push test * fixed vtx test * rename tolerance * timing verbosity * more timers * moved advection test * rename advection test * combined reconstruction and advection test * num push * more timers * modify test * api test * function timers * added new files * calcsurfDispIncr * advection passing * pumipic timing * set Proc Wedges * fixed mpi test perlmutter * set owning proc * fixed owning proc * timers * fixed error * cap normalize * clip normalized * remake fill mesh field * Debug statement, run for just 1 step * Debug statemensts in reconstruction * elem of MP as before * Num_elems_with_mps_before and after tracking * Bug of distance, element to element connectivity, and indexing -1 * Remove log file created erroneously * Used offset in adajcency calculated before * Reverting to use of -1 instead of offset for setting elment to element connection * E2E connetion not set a pplication level, -1 at tracking for now * fixed migration bug * removed comments whitespace * removed whitespace * remove unsued warnings * removed debug code * Merge yg/elem_centre to ac/reconstrcuction * Not calculating old detla averaged centroid * Changes requested in PR * Minor formats * PR changes 2 * Option to reverse lat lon increment * Check difference initial and final pos * Fixed spherical interp bug, quantiatative push ctest * Assert statement added * Linear Reconstruction * Linear Construction working, cleaner * enable c * divide advection tests * mpi fortran migration test * typo * Cleaned, Linear Reconstruction with regularization * Some more cleaning * Added rotated cordinates, added linear reconstruction test * Some more cleaning * Stricter tolerance * Linear Reconstruction test commented * PR changes 1 * Linear rec with MPs as Cell Centres * PR Changes 2 * PR Change 3 * renaming * set mpi communicator * assert spherical * linear typo * unimplemented code fix * rename variabel * change version --------- Co-authored-by: RPIFisherman <yuyanggong.rpi@gmail.com> Co-authored-by: dhyan1272 <nathd@rpi.edu> Co-authored-by: onkarsahni <sahni@rpi.edu>
1 parent efa1062 commit 4d9e866

34 files changed

Lines changed: 2254 additions & 495 deletions

‎.github/workflows/cmake.yml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -235,6 +235,7 @@ jobs:
235235
-DCMAKE_BUILD_TYPE="Debug"
236236
-DKokkos_DIR=${{ runner.temp }}/build-kokkos/install/lib64/cmake/Kokkos
237237
-DCMAKE_PREFIX_PATH=${{ runner.temp }}/build-pumi-pic/install
238+
-DIS_TESTING=on
238239
-DCMAKE_CXX_COMPILER=mpicxx
239240
-DCMAKE_INSTALL_PREFIX=${{ runner.temp }}/build-polyMPO/install
240241

‎CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
cmake_minimum_required(VERSION 3.19)
2-
project(polyMPO VERSION 0.0.4 LANGUAGES CXX Fortran)
2+
project(polyMPO VERSION 0.1.0 LANGUAGES CXX Fortran)
33

44
include(cmake/bob.cmake)
55
bob_begin_package()

‎src/CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@ set(HEADERS
44
pmpo_wachspressBasis.hpp
55
pmpo_MPMesh.hpp
66
pmpo_utils.hpp
7-
pmpo_assembly.hpp
7+
pmpo_MPMesh_assembly.hpp
88
pmpo_c.h
99
pmpo_createTestMPMesh.hpp
1010
)

‎src/pmpo_MPMesh.cpp‎

Lines changed: 105 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,44 @@ namespace polyMPO{
77

88
void printVTP_mesh(MPMesh& mpMesh, int printVTPIndex=-1);
99

10+
void MPMesh::calcBasis() {
11+
assert(p_mesh->getGeomType() == geom_spherical_surf);
12+
auto MPsPosition = p_MPs->getPositions();
13+
auto mp_basis_field = p_MPs->getData<MPF_Basis_Vals>(); // we can implement MPs->getBasisVals() like MPs->getPositions()
14+
auto elm2VtxConn = p_mesh->getElm2VtxConn();
15+
auto vtxCoords = p_mesh->getMeshField<MeshF_VtxCoords>();
16+
double radius = p_mesh->getSphereRadius();
17+
18+
auto calcbasis = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
19+
if(mask) { //if material point is 'active'/'enabled'
20+
Vec3d position3d(MPsPosition(mp,0),MPsPosition(mp,1),MPsPosition(mp,2));
21+
// formating
22+
Vec3d v3d[maxVtxsPerElm+1];
23+
int numVtx = elm2VtxConn(elm,0);
24+
for(int i = 1; i<=numVtx; i++){
25+
v3d[i-1][0] = vtxCoords(elm2VtxConn(elm,i)-1,0);
26+
v3d[i-1][1] = vtxCoords(elm2VtxConn(elm,i)-1,1);
27+
v3d[i-1][2] = vtxCoords(elm2VtxConn(elm,i)-1,2);
28+
}
29+
v3d[numVtx][0] = vtxCoords(elm2VtxConn(elm,1)-1,0);
30+
v3d[numVtx][1] = vtxCoords(elm2VtxConn(elm,1)-1,1);
31+
v3d[numVtx][2] = vtxCoords(elm2VtxConn(elm,1)-1,2);
32+
33+
double basisByArea3d[maxVtxsPerElm] = {0.0};
34+
initArray(basisByArea3d,maxVtxsPerElm,0.0);
35+
36+
// calc basis
37+
getBasisByAreaGblFormSpherical2(position3d, numVtx, v3d, radius, basisByArea3d);
38+
39+
// fill step
40+
for(int i=0; i<= numVtx; i++){
41+
mp_basis_field(mp,i) = basisByArea3d[i];
42+
}
43+
}
44+
};
45+
p_MPs->parallel_for(calcbasis, "calcbasis");
46+
}
47+
1048
void MPMesh::CVTTrackingEdgeCenterBased(Vec2dView dx){
1149
int numElms = p_mesh->getNumElements();
1250

@@ -74,36 +112,26 @@ void MPMesh::CVTTrackingEdgeCenterBased(Vec2dView dx){
74112

75113

76114
void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
115+
Kokkos::Timer timer;
77116
int numVtxs = p_mesh->getNumVertices();
78117
int numElms = p_mesh->getNumElements();
79118
auto numMPs = p_MPs->getCount();
80119

81-
const auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
120+
const auto elmCenter = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterXYZ>();
121+
82122
auto elm2VtxConn = p_mesh->getElm2VtxConn();
83123
auto elm2ElmConn = p_mesh->getElm2ElmConn();
84124

85125
auto mpPositions = p_MPs->getData<MPF_Cur_Pos_XYZ>();
86126
auto mpTgtPos = p_MPs->getData<MPF_Tgt_Pos_XYZ>();
87-
auto MPs2Elm = p_MPs->getData<MPF_Tgt_Elm_ID>();;
127+
auto MPs2Elm = p_MPs->getData<MPF_Tgt_Elm_ID>();
128+
auto MPs2Proc = p_MPs->getData<MPF_Tgt_Proc_ID>();
129+
auto elm2Process = p_mesh->getElm2Process();
88130

89131
if(printVTPIndex>=0) {
90132
printVTP_mesh(printVTPIndex);
91133
}
92134

93-
Vec3dView elmCenter("elementCenter",numElms);
94-
Kokkos::parallel_for("calcElementCenter", numElms, KOKKOS_LAMBDA(const int elm){
95-
int numVtx = elm2VtxConn(elm,0);
96-
double sum_x = 0.0, sum_y = 0.0, sum_z = 0.0;
97-
for(int i=1; i<= numVtx; i++){
98-
sum_x += vtxCoords(elm2VtxConn(elm,i)-1,0);
99-
sum_y += vtxCoords(elm2VtxConn(elm,i)-1,1);
100-
sum_z += vtxCoords(elm2VtxConn(elm,i)-1,2);
101-
}
102-
elmCenter(elm)[0] = sum_x/numVtx;
103-
elmCenter(elm)[1] = sum_y/numVtx;
104-
elmCenter(elm)[2] = sum_z/numVtx;
105-
});
106-
107135
Vec3dView history("positionHistory",numMPs);
108136
Vec3dView resultLeft("positionResult",numMPs);
109137
Vec3dView resultRight("positionResult",numMPs);
@@ -116,22 +144,32 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
116144
Vec3d MPnew(mpTgtPos(mp,0),mpTgtPos(mp,1),mpTgtPos(mp,2));
117145
Vec3d dx = MPnew-MP;
118146
while(true){
119-
int numVtx = elm2VtxConn(iElm,0);
120-
Vec3d delta = MPnew - elmCenter(iElm);
121-
double minDistSq = delta[0]*delta[0] + delta[1]*delta[1];
147+
int numConnElms = elm2ElmConn(iElm,0);
148+
149+
Vec3d center(elmCenter(iElm, 0), elmCenter(iElm, 1), elmCenter(iElm, 2));
150+
Vec3d delta = MPnew - center;
151+
152+
double minDistSq = delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2];
122153
int closestElm = -1;
123154
//go through all the connected elm, calc distance
124-
for(int i=1; i<=numVtx; i++){
125-
int elmID = elm2ElmConn(iElm,i);
126-
delta = MPnew - elmCenter(elmID);
127-
double neighborDistSq = delta[0]*delta[0] + delta[1]*delta[1];
155+
for(int i=1; i<=numConnElms; i++){
156+
int elmID = elm2ElmConn(iElm,i)-1;
157+
158+
//New delta
159+
Vec3d center(elmCenter(elmID, 0), elmCenter(elmID, 1), elmCenter(elmID, 2));
160+
delta = MPnew - center;
161+
162+
double neighborDistSq = delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2];
128163
if(neighborDistSq < minDistSq){
129164
closestElm = elmID;
130165
minDistSq = neighborDistSq;
131166
}
132167
}
168+
133169
if(closestElm<0){
134170
MPs2Elm(mp) = iElm;
171+
if (elm2Process.size() > 0)
172+
MPs2Proc(mp) = elm2Process(iElm);
135173
break;
136174
}else{
137175
iElm = closestElm;
@@ -193,6 +231,7 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
193231
fprintf(pFile," </DataArray>\n </Lines>\n </Piece>\n </PolyData>\n</VTKFile>\n");
194232
fclose(pFile);
195233
}
234+
pumipic::RecordTime("PolyMPO_CVTTrackingElmCenterBased", timer.seconds());
196235
}
197236

198237
void MPMesh::T2LTracking(Vec2dView dx){
@@ -258,17 +297,53 @@ void MPMesh::T2LTracking(Vec2dView dx){
258297
p_MPs->parallel_for(T2LCalc,"T2lTrackingCalc");
259298
}
260299

300+
void MPMesh::reconstructSlices() {
301+
if (reconstructSlice.size() == 0) return;
302+
Kokkos::Timer timer;
303+
calcBasis();
304+
for (auto const& [index, reconstruct] : reconstructSlice) {
305+
if (reconstruct) reconstruct();
306+
}
307+
reconstructSlice.clear();
308+
pumipic::RecordTime("PolyMPO_Reconstruct", timer.seconds());
309+
}
310+
311+
bool getAnyIsMigrating(MaterialPoints* p_MPs, bool isMigrating) {
312+
Kokkos::Timer timer;
313+
MPI_Comm comm = p_MPs->getMPIComm();
314+
int comm_rank;
315+
MPI_Comm_rank(comm, &comm_rank);
316+
int comm_size;
317+
MPI_Comm_size(comm, &comm_size);
318+
319+
bool anyIsMigrating = false;
320+
MPI_Allreduce(&isMigrating, &anyIsMigrating, 1, MPI_C_BOOL, MPI_LOR, comm);
321+
pumipic::RecordTime("PolyMPO_getAnyIsMigrating", timer.seconds());
322+
return anyIsMigrating;
323+
}
324+
261325
void MPMesh::push(){
326+
Kokkos::Timer timer;
262327
p_mesh->computeRotLatLonIncr();
263-
sphericalInterpolation<MeshF_RotLatLonIncr, MPF_Rot_Lat_Lon_Incr>(*this);
328+
sphericalInterpolation<MeshF_RotLatLonIncr>(*this);
264329
p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius()); // set Tgt_XYZ
330+
auto elm2Process = p_mesh->getElm2Process();
265331

266-
CVTTrackingElmCenterBased(); // move to Tgt_XYZ
332+
bool anyIsMigrating = false;
333+
do {
334+
CVTTrackingElmCenterBased(); // move to Tgt_XYZ
335+
p_MPs->updateMPSlice<MPF_Cur_Pos_XYZ, MPF_Tgt_Pos_XYZ>(); // Tgt_XYZ becomes Cur_XYZ
336+
p_MPs->updateMPSlice<MPF_Cur_Pos_Rot_Lat_Lon, MPF_Tgt_Pos_Rot_Lat_Lon>(); // Tgt becomes Cur
337+
if (elm2Process.size() > 0)
338+
anyIsMigrating = getAnyIsMigrating(p_MPs, p_MPs->migrate());
339+
else
340+
p_MPs->rebuild(); //rebuild pumi-pic
341+
p_MPs->updateMPElmID(); //update mpElm IDs slices
342+
reconstructSlices();
343+
}
344+
while (anyIsMigrating);
267345

268-
p_MPs->updateMPSlice<MPF_Cur_Pos_XYZ, MPF_Tgt_Pos_XYZ>(); // Tgt_XYZ becomes Cur_XYZ
269-
p_MPs->updateMPSlice<MPF_Cur_Pos_Rot_Lat_Lon, MPF_Tgt_Pos_Rot_Lat_Lon>(); // Tgt becomes Cur
270-
p_MPs->rebuild(); //rebuild pumi-pic
271-
p_MPs->updateMPElmID(); //update mpElm IDs slices
346+
pumipic::RecordTime("PolyMPO_push", timer.seconds());
272347
}
273348

274349
void MPMesh::printVTP_mesh(int printVTPIndex){
@@ -282,7 +357,7 @@ void MPMesh::printVTP_mesh(int printVTPIndex){
282357
FILE * pFile = fopen(fileOutput,"w");
283358
free(fileOutput);
284359

285-
DoubleVec3dView::HostMirror h_vtxCoords = Kokkos::create_mirror_view(vtxCoords);
360+
auto h_vtxCoords = Kokkos::create_mirror_view(vtxCoords);
286361
IntVtx2ElmView::HostMirror h_elm2VtxConn = Kokkos::create_mirror_view(elm2VtxConn);
287362
const int nCells = p_mesh->getNumElements();
288363
const int nVertices = p_mesh->getNumVertices();

‎src/pmpo_MPMesh.hpp‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,12 +7,20 @@
77

88
namespace polyMPO{
99

10+
template <MeshFieldIndex> const MaterialPointSlice meshFieldIndexToMPSlice;
11+
template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_Vel > = MPF_Vel;
12+
template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_VtxMass > = MPF_Mass;
13+
template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_ElmMass > = MPF_Mass;
14+
template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_RotLatLonIncr > = MPF_Rot_Lat_Lon_Incr;
15+
1016
#define maxMPsPerElm 8
1117

1218
class MPMesh{
1319
public:
1420
Mesh* p_mesh;
1521
MaterialPoints* p_MPs;
22+
23+
std::map<MeshFieldIndex, std::function<void()>> reconstructSlice = std::map<MeshFieldIndex, std::function<void()>>();
1624

1725
MPMesh(Mesh* inMesh, MaterialPoints* inMPs):
1826
p_mesh(inMesh), p_MPs(inMPs) {
@@ -26,6 +34,25 @@ class MPMesh{
2634
void CVTTrackingElmCenterBased(const int printVTPIndex = -1);
2735
void T2LTracking(Vec2dView dx);
2836
void push();
37+
void calcBasis();
38+
39+
DoubleView assemblyV0();
40+
template <MaterialPointSlice index>
41+
DoubleView wtScaAssembly();
42+
template <MaterialPointSlice index>
43+
Vec2dView wtVec2Assembly();
44+
template <MeshFieldIndex meshFieldIndex>
45+
void assembly(int order, MeshFieldType type, bool basisWeightFlag, bool massWeightFlag);
46+
template <MeshFieldIndex meshFieldIndex>
47+
void assemblyVtx0();
48+
template <MeshFieldIndex meshFieldIndex>
49+
void assemblyElm0();
50+
template <MeshFieldIndex meshFieldIndex>
51+
void assemblyVtx1();
52+
53+
template<MeshFieldIndex meshFieldIndex>
54+
void setReconstructSlice(int order, MeshFieldType type);
55+
void reconstructSlices();
2956

3057
void printVTP_mesh(int printVTPIndex);
3158
};

0 commit comments

Comments
 (0)