Skip to content

Commit 59bb2ca

Browse files
dhyan1272dhyan1272
committed
Stress Divergence Consistent with MPAS
1 parent 639b1fb commit 59bb2ca

1 file changed

Lines changed: 31 additions & 50 deletions

File tree

‎src/pmpo_MPMesh.cpp‎

Lines changed: 31 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -16,14 +16,10 @@ void MPMesh::calculateStrain(){
1616
auto MPsStrainRate = p_MPs->getData<MPF_Strain_Rate>();
1717
//Mesh Fields
1818
auto tanLatVertexRotatedOverRadius = p_mesh->getMeshField<MeshF_TanLatVertexRotatedOverRadius>();
19-
auto gnomProjVtx = p_mesh->getMeshField<MeshF_VtxGnomProj>();
20-
auto gnomProjElmCenter = p_mesh->getMeshField<MeshF_ElmCenterGnomProj>();
2119
auto elm2VtxConn = p_mesh->getElm2VtxConn();
2220
auto velField = p_mesh->getMeshField<MeshF_Vel>();
2321
auto solveStress = p_mesh->getMeshField<polyMPO::MeshF_SolveStress>();
2422

25-
bool isRotated = p_mesh->getRotatedFlag();
26-
2723
auto setMPStrainRate = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
2824
if(mask){
2925

@@ -36,16 +32,6 @@ void MPMesh::calculateStrain(){
3632

3733
int numVtx = elm2VtxConn(elm,0);
3834

39-
Vec3d position3d(MPsPosition(mp,0),MPsPosition(mp,1),MPsPosition(mp,2));
40-
if(isRotated){
41-
position3d[0] = -MPsPosition(mp, 2);
42-
position3d[2] = MPsPosition(mp, 0);
43-
}
44-
auto gnomProjElmCenter_sub = Kokkos::subview(gnomProjElmCenter, elm, Kokkos::ALL);
45-
double mpProjX, mpProjY;
46-
computeGnomonicProjectionAtPoint(position3d, gnomProjElmCenter_sub, mpProjX, mpProjY);
47-
auto gnom_vtx_subview = Kokkos::subview(gnomProjVtx, elm, Kokkos::ALL, Kokkos::ALL);
48-
4935
double v11 = 0.0;
5036
double v12 = 0.0;
5137
double v21 = 0.0;
@@ -134,65 +120,60 @@ void MPMesh::calculateStressDivergence(){
134120
if(p_mesh->getGeomType() == geom_spherical_surf)
135121
radius=p_mesh->getSphereRadius();
136122

137-
//Reconstructed the stress
138-
Kokkos::View<vec2d_t*> stress_divUV("stress_divUV", p_mesh->getNumVertices());
139-
Kokkos::View<vec2d_t*> divUV_edge("divUVedge", p_mesh->getNumVertices());
123+
auto stress_divUV = p_mesh->getMeshField<MeshF_StressDivergence>();
124+
Kokkos::deep_copy(stress_divUV, 0.0);
140125

141126
//Assemble fields for Stress Divergence
142127
auto stress_div = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
143128
if(mask) { //if material point is 'active'/'enabled'
144129
int nVtxE = elm2VtxConn(elm,0); //number of vertices bounding the element
145130
for(int i=0; i<nVtxE; i++){
146131
int vID = elm2VtxConn(elm,i+1)-1;
132+
133+
double ramp = nearAnEdge_l(vID);
134+
double invM = 1.0/vtxMatrixMass_l(vID);
135+
invM = vtxMatrixMass_l(vID) >1e-4 ? invM : 0;
136+
147137
double w_vtx=weight(mp,i);
148138
double CoordDiffs[vec4d_nEntries] = {1, (-vtxCoords(vID,0) + mpPositions(mp,0))/radius,
149139
(-vtxCoords(vID,1) + mpPositions(mp,1))/radius,
150140
(-vtxCoords(vID,2) + mpPositions(mp,2))/radius};
151141

152-
auto factor = w_vtx*(VtxCoeffs_new(vID,0, 0) + VtxCoeffs_new(vID,0, 1)*CoordDiffs[1] +
153-
VtxCoeffs_new(vID,0, 2)*CoordDiffs[2] +
154-
VtxCoeffs_new(vID,0, 3)*CoordDiffs[3]);
142+
auto factor = ramp * w_vtx * (VtxCoeffs_new(vID,0, 0) + VtxCoeffs_new(vID,0, 1)*CoordDiffs[1] +
143+
VtxCoeffs_new(vID,0, 2)*CoordDiffs[2] +
144+
VtxCoeffs_new(vID,0, 3)*CoordDiffs[3]) +
145+
(1.0 - ramp) * invM * w_vtx;
146+
147+
factor = factor * tanLatVertexRotatedOverRadius(vID, 0);
148+
149+
auto factor1 = ramp * (w_vtx/radius) * (VtxCoeffs_new(vID, 1, 0) + VtxCoeffs_new(vID, 1, 1)*CoordDiffs[1] +
150+
VtxCoeffs_new(vID, 1, 2)*CoordDiffs[2] +
151+
VtxCoeffs_new(vID, 1, 3)*CoordDiffs[3]) -
152+
(1.0 - ramp) * invM * weight_grads(mp, i*2+0);
155153

156-
auto factor1 = (w_vtx/radius)*(VtxCoeffs_new(vID,1, 0) + VtxCoeffs_new(vID,1, 1)*CoordDiffs[1] +
157-
VtxCoeffs_new(vID,1, 2)*CoordDiffs[2] +
158-
VtxCoeffs_new(vID,1, 3)*CoordDiffs[3]);
159-
auto factor2 = (w_vtx/radius)*(VtxCoeffs_new(vID,2, 0) + VtxCoeffs_new(vID,2, 1)*CoordDiffs[1] +
160-
VtxCoeffs_new(vID,2, 2)*CoordDiffs[2] +
161-
VtxCoeffs_new(vID,2, 3)*CoordDiffs[3]);
154+
auto factor2 = ramp * (w_vtx/radius) * (VtxCoeffs_new(vID, 2, 0) + VtxCoeffs_new(vID, 2, 1)*CoordDiffs[1] +
155+
VtxCoeffs_new(vID, 2, 2)*CoordDiffs[2] +
156+
VtxCoeffs_new(vID, 2, 3)*CoordDiffs[3]) -
157+
(1.0 - ramp) * invM * weight_grads(mp, i*2+1);
162158

163159
Kokkos::atomic_add(&stress_divUV(vID, 0), factor1 * MPsStress(mp, 0) + factor2 * MPsStress(mp, 2) -
164-
2 * tanLatVertexRotatedOverRadius(vID, 0) * factor * MPsStress(mp, 2));
165-
Kokkos::atomic_add(&stress_divUV(vID, 1), factor2 * MPsStress(mp, 1) + factor1*MPsStress(mp, 2) +
166-
factor * tanLatVertexRotatedOverRadius(vID, 0) * (MPsStress(mp, 0)-MPsStress(mp, 1)));
167-
168-
Kokkos::atomic_add(&divUV_edge(vID, 0), - weight_grads(mp, i*2 + 0) * MPsStress(mp, 0) - weight_grads(mp, i*2 + 1) * MPsStress(mp, 2) -
169-
2 * tanLatVertexRotatedOverRadius(vID, 0) * w_vtx * MPsStress(mp, 2));
170-
Kokkos::atomic_add(&divUV_edge(vID, 1), - weight_grads(mp, i*2 + 1) * MPsStress(mp, 1) - weight_grads(mp, i*2 + 0) * MPsStress(mp, 2) +
171-
w_vtx * tanLatVertexRotatedOverRadius(vID, 0) * (MPsStress(mp, 0)- MPsStress(mp, 1)));
160+
2 * factor * MPsStress(mp, 2));
161+
162+
Kokkos::atomic_add(&stress_divUV(vID, 1), factor2 * MPsStress(mp, 1) + factor1 * MPsStress(mp, 2) +
163+
factor * (MPsStress(mp, 0)-MPsStress(mp, 1)));
172164
}
173165
}
174166
};
175167
p_MPs->parallel_for(stress_div, " stress_div_assembly");
176168

169+
/*
177170
if(numProcsTot>1){
178171
//Takes contribution of halo vertices and adds it in owner procs
179172
communicate_and_take_halo_contributions(stress_divUV, numVertices, 2, 0, 0);
180-
communicate_and_take_halo_contributions(divUV_edge, numVertices, 2, 0, 0);
181173
//Transfer the correct values at owned vertices to halo vertices
182174
communicate_and_take_halo_contributions(stress_divUV, numVertices, 2, 1, 1);
183-
communicate_and_take_halo_contributions(divUV_edge, numVertices, 2, 1, 1);
184175
}
185-
186-
auto stressDivergence = p_mesh->getMeshField<MeshF_StressDivergence>();
187-
Kokkos::parallel_for("calculate_divergence", numVtxOwned, KOKKOS_LAMBDA(const int vtx){
188-
double ramp = nearAnEdge_l(vtx);
189-
double invM = 1.0/vtxMatrixMass_l(vtx);
190-
invM = vtxMatrixMass_l(vtx) >1e-4 ? invM : 0;
191-
192-
stressDivergence(vtx, 0) = ramp * stress_divUV(vtx, 0) + (1 - ramp) * divUV_edge(vtx, 0) * invM ;
193-
stressDivergence(vtx, 1) = ramp * stress_divUV(vtx, 1) + (1 - ramp) * divUV_edge(vtx, 1) * invM ;
194-
});
195-
176+
*/
196177
pumipic::RecordTime("Stress_Divergence_Reconstruction" + std::to_string(self), timer.seconds());
197178
}
198179

@@ -281,7 +262,7 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
281262
auto elm2global = p_mesh->getElmGlobal();
282263

283264
if(printVTPIndex>=0) {
284-
printVTP_mesh(printVTPIndex);
265+
printVTP_mesh(printVTPIndex);
285266
}
286267

287268
Vec3dView history("positionHistory",numMPs);
@@ -383,11 +364,11 @@ bool getAnyIsMigrating(MaterialPoints* p_MPs, bool isMigrating) {
383364

384365
void MPMesh::push_ahead(){
385366
Kokkos::Timer timer;
386-
//Latitude Longitude increment at mesh vertices and interpolate to particle position
367+
//Latitude Longitude increment at mesh vertices
387368
p_mesh->computeRotLatLonIncr();
388369

389370
//Interpolates latitude longitude, mesh velocity increments to MPs
390-
calcBasis();
371+
//calcBasis();
391372
sphericalInterpolation<MeshF_RotLatLonIncr>(*this);
392373
sphericalInterpolation<MeshF_OnSurfVeloIncr>(*this);
393374
sphericalInterpolation2Fields(*this);

0 commit comments

Comments
 (0)