Skip to content

Port particle-cell mapping to GPU #91

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

Open
Shahrear2000 wants to merge 17 commits into
base: sj/gpu-aware-mpi-halo-exchange
Choose a base branch
from
Open

Port particle-cell mapping to GPU #91

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
b98b69a
First commit to port
Jun 29, 2026
12417d3
Conditional for HALO exchange of nProcsTot=1 in polyMPO
Jun 30, 2026
229c8cb
SetMPArea and modify in strain rate
Jul 13, 2026
f0c25f7
Passing oceanDragCoeff and using same vel corrections as rebase 29 June
Jul 14, 2026
3b6fccd
Reconstrucntion of velocity test commented out as now in src involves…
Jul 14, 2026
215764a
More formatting
Jul 15, 2026
10c3a39
Particle to cell mapping for a scalar field
Aug 17, 2026
72f5916
Cell to MP mapping also
Sep 3, 2026
8df665e
Removed un-necessay files
Sep 3, 2026
30d1ba7
Clean up of code
Sep 7, 2026
71e7aea
Get the mapped fields back to MPAS
Sep 8, 2026
6b4eb2a
Intent of c_ptr
Sep 8, 2026
28b7fbe
Merge remote-tracking branch 'origin/dn/Particle-cell-Mapping' into s…
Shahrear2000 Sep 26, 2026
044c955
Particle-cell-map port to GPU/polyMPO - category fields
Shahrear2000 Sep 26, 2026
d782960
Add particle-cell mapping unit test
Shahrear2000 Sep 29, 2026
7c3c088
Run particle-cell mapping unit test on the spherical SCVT mesh
Shahrear2000 Sep 29, 2026
6418bf7
Add unit testcase with SCVT mesh
Shahrear2000 Sep 30, 2026
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
27 changes: 17 additions & 10 deletions src/pmpo_MPMesh.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,11 +14,13 @@ void MPMesh::calculateStrain(){
auto MPsBasisGrads = p_MPs->getData<MPF_Basis_Grad_Vals>();
auto MPsAppID = p_MPs->getData<MPF_MP_APP_ID>();
auto MPsStrainRate = p_MPs->getData<MPF_Strain_Rate>();
auto MPsArea = p_MPs->getData<polyMPO::MPF_Area>();
//Mesh Fields
auto tanLatVertexRotatedOverRadius = p_mesh->getMeshField<MeshF_TanLatVertexRotatedOverRadius>();
auto elm2VtxConn = p_mesh->getElm2VtxConn();
auto velField = p_mesh->getMeshField<MeshF_Vel>();
auto solveStress = p_mesh->getMeshField<polyMPO::MeshF_SolveStress>();
auto elasticTimeStep = p_mesh->getElasticTimeStep();

auto setMPStrainRate = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
if(mask){
Expand All @@ -38,7 +40,7 @@ void MPMesh::calculateStrain(){
double v22 = 0.0;
double uTanOverR = 0.0;
double vTanOverR = 0.0;

for (int i = 0; i < numVtx; i++){
int iVertex = elm2VtxConn(elm, i+1)-1;
v11 = v11 + MPsBasisGrads(mp, i*2 + 0) * velField(iVertex, 0);
Expand All @@ -52,6 +54,8 @@ void MPMesh::calculateStrain(){
MPsStrainRate(mp, 0) = v11 - vTanOverR;
MPsStrainRate(mp, 1) = v22;
MPsStrainRate(mp, 2) = 0.5*(v12 + v21 + uTanOverR);

MPsArea(mp, 0) = MPsArea(mp, 0) * exp((v11+v22-vTanOverR)*elasticTimeStep);
}
};
p_MPs->parallel_for(setMPStrainRate, "setMPStrainRate");
Expand All @@ -75,13 +79,16 @@ void MPMesh::calculateStress(const int constitutive_relation){
if(mask){
Vec3d strain_rate (MPsStrainRate(mp, 0), MPsStrainRate(mp, 1), MPsStrainRate(mp, 2));
Vec3d stress(MPsStress(mp, 0), MPsStress(mp, 1), MPsStress(mp, 2));
double rep_pressure=MPsRepPressure(mp,0);

if (constitutive_relation == 1)
constitutive_evp(strain_rate, stress, MPsIcePressure(mp,0), MPsRepPressure(mp,0), MPsArea(mp,0), elasticTimeStep, dampingTimescale);
constitutive_evp(strain_rate, stress, MPsIcePressure(mp,0), rep_pressure, MPsArea(mp,0), elasticTimeStep, dampingTimescale);
else if(constitutive_relation == 3)
constitutive_linear(strain_rate, stress);

for (int m=0 ; m<3; m++)
MPsStress(mp, m) = stress[m]*solveStress(elm);
MPsRepPressure(mp,0)=rep_pressure;
}
};
p_MPs->parallel_for(setMPStress, "setMPStress");
Expand Down Expand Up @@ -145,7 +152,7 @@ void MPMesh::calculateStressDivergence(){
(1.0 - ramp) * invM * w_vtx;

factor = factor * tanLatVertexRotatedOverRadius(vID, 0);

auto factor1 = ramp * (w_vtx/radius) * (VtxCoeffs_new(vID, 1, 0) + VtxCoeffs_new(vID, 1, 1)*CoordDiffs[1] +
VtxCoeffs_new(vID, 1, 2)*CoordDiffs[2] +
VtxCoeffs_new(vID, 1, 3)*CoordDiffs[3]) -
Expand Down Expand Up @@ -280,7 +287,7 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
Vec3d dx = MPnew-MP;
while(true){
int numConnElms = elm2ElmConn(iElm,0);

Vec3d center(elmCenter(iElm, 0), elmCenter(iElm, 1), elmCenter(iElm, 2));
Vec3d delta = MPnew - center;

Expand Down Expand Up @@ -742,8 +749,8 @@ void MPMesh::T2LTracking(Vec2dView dx){
Vec2d MP(mpPositions(mp,0),mpPositions(mp,1));//XXX:the input is XYZ, but we only support 2d vector
if(mask){
int iElm = elm;
Vec2d MPnew = MP + dx(mp);
Vec2d MPnew = MP + dx(mp);

while(true){
int numVtx = elm2VtxConn(iElm,0);
bool goToNeighbour = false;
Expand All @@ -753,25 +760,25 @@ void MPMesh::T2LTracking(Vec2dView dx){
v[i] = elm2VtxConn(iElm,i+1)-1;
//get edges and perpendiculardx
Vec2d e[maxVtxsPerElm];
double pdx[maxVtxsPerElm];
double pdx[maxVtxsPerElm];
for(int i=0; i< numVtx; i++){
int idx_ip1 = (i+1)%numVtx;
Vec2d v_i(vtxCoords(v[i],0),vtxCoords(v[i],1));
Vec2d v_ip1(vtxCoords(v[idx_ip1],0),vtxCoords(v[idx_ip1],1));
e[i] = v_ip1 - v_i;
pdx[i] = (v_i - MP).cross(dx(mp));
}

for(int i=0; i<numVtx; i++){
int ip1 = (i+1)%numVtx;
//pdx*pdx<0 and edge is acrossed
//pdx*pdx<0 and edge is acrossed
if(pdx[i]*pdx[ip1] <0 && e[i].cross(Vec2d(MPnew[0]-vtxCoords(v[i],0),
MPnew[1]-vtxCoords(v[i],1)))<0){
//go to the next elm
iElm = elm2ElmConn(iElm,i+1);
goToNeighbour = true;
if(iElm <0){
mpStatus(mp) = 0;
mpStatus(mp) = 0;
MPs2Elm(mp) = -1;
goToNeighbour = false;
}
Expand Down
224 changes: 216 additions & 8 deletions src/pmpo_MPMesh.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,19 @@ template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_VtxMass
template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_ElmMass > = MPF_Mass;
template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_RotLatLonIncr > = MPF_Rot_Lat_Lon_Incr;
template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_OnSurfVeloIncr > = MPF_Vel_Incr;
template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_OceanStressCell> = MPF_OceanStress;

template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_IceAreaCategory > = MPF_IceAreaCategory;
template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_IceVolumeCategory > = MPF_IceVolumeCategory;
template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_SnowVolumeCategory > = MPF_SnowVolumeCategory;

template <MaterialPointSlice>
const MeshFieldIndex MPSliceToMeshFieldIndex;
template <> const MeshFieldIndex MPSliceToMeshFieldIndex < MPF_OpenWaterArea > = MeshF_OpenWaterArea;

template <> const MeshFieldIndex MPSliceToMeshFieldIndex < MPF_IceAreaCategory > = MeshF_IceAreaCategory;
template <> const MeshFieldIndex MPSliceToMeshFieldIndex < MPF_IceVolumeCategory > = MeshF_IceVolumeCategory;
template <> const MeshFieldIndex MPSliceToMeshFieldIndex < MPF_SnowVolumeCategory > = MeshF_SnowVolumeCategory;

#define maxMPsPerElm 8

Expand All @@ -43,9 +56,6 @@ class MPMesh{

void startCommunication();

void communicate_and_take_halo_contributions(const Kokkos::View<double**>& meshField, int nEntities, int numEntries, int mode, int op);

//Now Kokkos views are made 1D
template <typename ViewType>
void communicate_and_take_halo_contributions_staged(
const ViewType& meshField,
Expand Down Expand Up @@ -134,10 +144,6 @@ class MPMesh{
pumipic::RecordTime("SD: Contribution" + std::to_string(self), timer.seconds());
}


void communicateFields(const std::vector<std::vector<double>>& fieldData, const int numEntities, const int numEntries, int mode,
std::vector<std::vector<int>>& recvIDVec, std::vector<std::vector<double>>& recvDataVec);

template <class ViewType>
void communicateFieldsFromHostView(
const ViewType& fieldData,
Expand Down Expand Up @@ -272,8 +278,9 @@ class MPMesh{
template <MeshFieldIndex meshFieldIndex>
void assemblyVtx1();
void reconstruct_coeff_full();
void invertMatrix(const Kokkos::View<double**>& vtxMatrices, const double& radius);
void invertMatrix(const Kokkos::View<double**>& vtxMatrices, const double& radius, bool vtxCoeffCalc);
Kokkos::View<double*[vec3d_nEntries][vec4d_nEntries]> precomputedVtxCoeffs_new;
Kokkos::View<double*[vec3d_nEntries][vec4d_nEntries]> precomputedElmCoeffs_new;
Kokkos::View<double*> nearAnEdge;
Kokkos::View<double*> vtxMatrixMass;

Expand Down Expand Up @@ -689,6 +696,207 @@ class MPMesh{

#endif

template<MaterialPointSlice MPSlice>
void mapMPsToCells(){
//MP field mesh field mapping
constexpr MeshFieldIndex meshFieldIndex = MPSliceToMeshFieldIndex<MPSlice>;
const int numEntries = mpSliceToNumEntries<MPSlice>();
auto mpField = p_MPs->getData<MPSlice>();
auto meshField = p_mesh->getMeshField<meshFieldIndex>();
Kokkos::deep_copy(meshField, 0.0);

//MP fields
auto mpArea = p_MPs->getData<MPF_Area>();
auto mpPos = p_MPs->getData<MPF_Cur_Pos_XYZ>();

//Mesh fields and finding number of MPs per cell and areaMP per Cell
auto nElms = p_mesh->getNumElements();
const auto elmCenter = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterXYZ>();
double radius = 1.0;
if(p_mesh->getGeomType() == geom_spherical_surf)
radius=p_mesh->getSphereRadius();

Kokkos::View<int*> nMPsPerCell("nMPsPerCell", nElms);
Kokkos::View<double*> sumAreaMP("sumAreaMP", nElms);
Kokkos::deep_copy(nMPsPerCell, 0);
Kokkos::deep_copy(sumAreaMP, 0.0);
auto calcMPsCell = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
if(mask){
Kokkos::atomic_increment(&nMPsPerCell(elm));
Kokkos::atomic_add(&sumAreaMP(elm), mpArea(mp, 0));
}
};
p_MPs->parallel_for(calcMPsCell, "calcN_MPsCell");

//TODO put them as inputs
bool higherOrderRemap = true;

//Calculate the matrix
constexpr int numEntriesMatrix=10;
Kokkos::View<double*[numEntriesMatrix]> lsMatrices("VtxMatrices", p_mesh->getNumElements());
Kokkos::deep_copy(lsMatrices, 0);

auto assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
if(mask) { //if material point is 'active'/'enabled'
if(nMPsPerCell(elm) <= 0) return;
double wp = mpArea(mp,0)/sumAreaMP(elm);
Kokkos::atomic_add(&lsMatrices(elm,0), wp);
Kokkos::atomic_add(&lsMatrices(elm,1), wp*(-elmCenter(elm,0)+mpPos(mp,0))/radius);
Kokkos::atomic_add(&lsMatrices(elm,2), wp*(-elmCenter(elm,1)+mpPos(mp,1))/radius);
Kokkos::atomic_add(&lsMatrices(elm,3), wp*(-elmCenter(elm,2)+mpPos(mp,2))/radius);
Kokkos::atomic_add(&lsMatrices(elm,4), wp*(-elmCenter(elm,0)+mpPos(mp,0)) *
(-elmCenter(elm,0)+mpPos(mp,0))/(radius*radius));
Kokkos::atomic_add(&lsMatrices(elm,5), wp*(-elmCenter(elm,0)+mpPos(mp,0)) * (
-elmCenter(elm,1)+mpPos(mp,1))/(radius*radius));
Kokkos::atomic_add(&lsMatrices(elm,6), wp*(-elmCenter(elm,0)+mpPos(mp,0)) *
(-elmCenter(elm,2)+mpPos(mp,2))/(radius*radius));
Kokkos::atomic_add(&lsMatrices(elm,7), wp*(-elmCenter(elm,1)+mpPos(mp,1)) *
(-elmCenter(elm,1)+mpPos(mp,1))/(radius*radius));
Kokkos::atomic_add(&lsMatrices(elm,8), wp*(-elmCenter(elm,1)+mpPos(mp,1)) *
(-elmCenter(elm,2)+mpPos(mp,2))/(radius*radius));
Kokkos::atomic_add(&lsMatrices(elm,9), wp*(-elmCenter(elm,2)+mpPos(mp,2)) *
(-elmCenter(elm,2)+mpPos(mp,2))/(radius*radius));
}
};
p_MPs->parallel_for(assemble, "assembly");

invertMatrix(lsMatrices, radius, false);
auto precomputedElmCoeffs_l = this->precomputedElmCoeffs_new;

Kokkos::View<double**> minVals("minVals", nElms, numEntries);
Kokkos::View<double**> maxVals("maxVals", nElms, numEntries);
Kokkos::deep_copy(minVals, std::numeric_limits<double>::max());
Kokkos::deep_copy(maxVals, std::numeric_limits<double>::lowest());

auto calcCellField= PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
if(mask) { //if material point is 'active'/'enabled
if(nMPsPerCell(elm) <= 0) return;
double wp = mpArea(mp,0)/sumAreaMP(elm);
Vec3d pXYZ = {mpPos(mp,0) - elmCenter(elm,0), mpPos(mp,1) - elmCenter(elm,1), mpPos(mp,2) - elmCenter(elm,2)};
double psi = wp * (precomputedElmCoeffs_l(elm, 0, 0) + (precomputedElmCoeffs_l(elm, 0, 1) * pXYZ[0] +
precomputedElmCoeffs_l(elm, 0, 2) * pXYZ[1] +
precomputedElmCoeffs_l(elm, 0, 3) * pXYZ[2])/radius);


for (int comp=0; comp<numEntries; comp++)
{
Kokkos::atomic_add(&meshField(elm, comp), mpField(mp,comp)*psi);
Kokkos::atomic_fetch_min(&minVals(elm,comp), mpField(mp,comp));
Kokkos::atomic_fetch_max(&maxVals(elm,comp), mpField(mp,comp));
}
}
};
p_MPs->parallel_for(calcCellField, "calcCellField");

Kokkos::fence();

Kokkos::parallel_for("clampCellField", nElms, KOKKOS_LAMBDA(const int elm){
if(nMPsPerCell(elm) <= 0) return;
for (int comp=0; comp<numEntries; comp++) {
double v = meshField(elm, comp);
if (v < minVals(elm,comp)) v = minVals(elm,comp);
if (v > maxVals(elm,comp)) v = maxVals(elm,comp);
meshField(elm, comp) = v;
}
});
}

template<MeshFieldIndex mfIndex>
void mapCellsToMPs(){

//MP Field
constexpr MaterialPointSlice mpSlice = meshFieldIndexToMPSlice<mfIndex>;
auto mpField = p_MPs->getData<mpSlice>();
auto mpPos = p_MPs->getData<MPF_Cur_Pos_XYZ>();

//Mesh Fields
auto nElms = p_mesh->getNumElements();
auto meshField = p_mesh->getMeshField<mfIndex>();
auto elm2ElmConn = p_mesh->getElm2ElmConn();
auto gnomProjVtx = p_mesh->getMeshField<MeshF_VtxGnomProj>();
auto gnomProjCell = p_mesh->getMeshField<MeshF_CellGnomProj>();
auto gnomProjCellOnCell = p_mesh->getMeshField<MeshF_CellOnCellGnomProj>();
auto gnomProjElmCenter = p_mesh->getMeshField<MeshF_ElmCenterGnomProj>();
bool isRotated = p_mesh->getRotatedFlag();

//Mps Per Cell
Kokkos::View<int*> nMPsPerCell("nMPsPerCell", nElms);
Kokkos::deep_copy(nMPsPerCell, 0);
auto calcMPsCell = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
if(mask){
Kokkos::atomic_increment(&nMPsPerCell(elm));
}
};
p_MPs->parallel_for(calcMPsCell, "calcN_MPsCell");

const int numEntries = mpSliceToNumEntries<mpSlice>();
for (int comp=0; comp<numEntries; comp++){

//Calculate mapped Value
Kokkos::View<double*> centerValues("centerValues", nElms);
Kokkos::View<double*> grads("grads", 2*nElms);
Kokkos::View<double*> cellMinValues("cellMinValues", nElms);
Kokkos::View<double*> cellMaxValues("cellMaxValues", nElms);

Kokkos::parallel_for("calcMatrix", nElms, KOKKOS_LAMBDA(const int elm){
if(nMPsPerCell(elm) <= 0) return;
double lsMatrix[3]={0.0};
double rsMatrix[2]={0.0};
double cellMin = meshField(elm, comp);
double cellMax = meshField(elm, comp);

int numConnElms = elm2ElmConn(elm, 0);
for(int i=1; i<=numConnElms; i++){
int elmID = elm2ElmConn(elm,i)-1;
if(elmID < 0 || elmID >= nElms)
continue;
lsMatrix[0] += (gnomProjCellOnCell(elm, i-1, 0) - gnomProjCell(elm, 0)) * (gnomProjCellOnCell(elm, i-1, 0) - gnomProjCell(elm, 0));
lsMatrix[1] += (gnomProjCellOnCell(elm, i-1, 0) - gnomProjCell(elm, 0)) * (gnomProjCellOnCell(elm, i-1, 1) - gnomProjCell(elm, 1));
lsMatrix[2] += (gnomProjCellOnCell(elm, i-1, 1) - gnomProjCell(elm, 1)) * (gnomProjCellOnCell(elm, i-1, 1) - gnomProjCell(elm, 1));
rsMatrix[0] += (meshField(elmID, comp)- meshField(elm, comp)) * (gnomProjCellOnCell(elm, i-1, 0) - gnomProjCell(elm, 0));
rsMatrix[1] += (meshField(elmID, comp)- meshField(elm, comp)) * (gnomProjCellOnCell(elm, i-1, 1) - gnomProjCell(elm, 1));
if (meshField(elmID, comp) < cellMin) cellMin = meshField(elmID, comp);
if (meshField(elmID, comp) > cellMax) cellMax = meshField(elmID, comp);
}
double det = lsMatrix[0] * lsMatrix[2] - lsMatrix[1] * lsMatrix[1];
double grad[2] = {0.0, 0.0};
if (std::abs(det) > 1.0e-12) {
grad[0] = (1.0 / det) * ( lsMatrix[2] * rsMatrix[0] - lsMatrix[1] * rsMatrix[1]);
grad[1] = (1.0 / det) * (-lsMatrix[1] * rsMatrix[0] + lsMatrix[0] * rsMatrix[1]);
}
double centerVal = meshField(elm, comp) - grad[0] * gnomProjCell(elm, 0) - grad[1] * gnomProjCell(elm, 1);

grads(2*elm + 0) = grad[0];
grads(2*elm + 1) = grad[1];
centerValues(elm) = centerVal;
cellMinValues(elm) = cellMin;
cellMaxValues(elm) = cellMax;
});

Kokkos::fence();

auto calcMPValue = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
if(mask){
if(nMPsPerCell(elm) <= 0) return;
Vec3d position3d(mpPos(mp, 0), mpPos(mp, 1), mpPos(mp, 2));
if(isRotated){
position3d[0] = -mpPos(mp, 2);
position3d[2] = mpPos(mp, 0);
}
double mpProjX, mpProjY;
auto gnomProjElmCenter_sub = Kokkos::subview(gnomProjElmCenter, elm, Kokkos::ALL);
computeGnomonicProjectionAtPoint(position3d, gnomProjElmCenter_sub, mpProjX, mpProjY);

double mpVal = centerValues(elm) + grads(2*elm + 0) * mpProjX + grads(2*elm + 1) * mpProjY;
if (mpVal < cellMinValues(elm)) mpVal = cellMinValues(elm);
if (mpVal > cellMaxValues(elm)) mpVal = cellMaxValues(elm);
mpField(mp, comp) = mpVal;
}
};
p_MPs->parallel_for(calcMPValue, "calcMPValue");
}
}

};

}//namespace polyMPO end
Expand Down
Loading
Loading