Skip to content

Commit 036fa3d

Browse files
dhyan1272dhyan1272
committed
Stress Divergence Verified
1 parent f525ddb commit 036fa3d

3 files changed

Lines changed: 78 additions & 19 deletions

File tree

‎src/pmpo_MPMesh.cpp‎

Lines changed: 19 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -54,7 +54,7 @@ void MPMesh::calculateStrain(){
5454
vTanOverR = vTanOverR + MPsBasis(mp, i) * tanLatVertexRotatedOverRadius(iVertex, 0) * velField(iVertex, 1);
5555
//Debugging
5656
if(MPsAppID(mp)==0){
57-
printf("Strain Calc: iVertex %d vel field %.15e %.15e \n", iVertex, velField(iVertex, 0), velField(iVertex, 1));
57+
printf("Strain Calc: iVertex %d vel field %.15e %.15e \n", iVertex, velField(iVertex, 0), velField(iVertex, 1));
5858
}
5959
}
6060
MPsStrainRate(mp, 0) = v11 - vTanOverR;
@@ -66,7 +66,10 @@ void MPMesh::calculateStrain(){
6666
}
6767

6868
void MPMesh::calculateStress(){
69-
//MeshFields
69+
static int int_xx =0;
70+
//std::cout<<"Counting stress: "<< int_xx << std::endl;
71+
72+
//MeshFields
7073
auto solveStress = p_mesh->getMeshField<polyMPO::MeshF_SolveStress>();
7174
auto elasticTimeStep = p_mesh->getElasticTimeStep();
7275
auto dynamicTimeStep = p_mesh->getDynamicTimeStep();
@@ -89,12 +92,13 @@ void MPMesh::calculateStress(){
8992
MPsStress(mp, m) = stress[m];
9093
//Debugging
9194
if(MPsAppID(mp)==0){
92-
printf("Strain in GPU: %.15e %.15e %.15e\n", MPsStrainRate(mp, 0), MPsStrainRate(mp, 1), MPsStrainRate(mp, 2));
93-
printf("Stress in GPU: %.15e %.15e %.15e\n", MPsStress(mp, 0), MPsStress(mp, 1), MPsStress(mp, 2));
95+
//printf("Strain in GPU: %.15e %.15e %.15e\n", MPsStrainRate(mp, 0), MPsStrainRate(mp, 1), MPsStrainRate(mp, 2));
96+
//printf("Stress in GPU: %.15e %.15e %.15e\n", MPsStress(mp, 0), MPsStress(mp, 1), MPsStress(mp, 2));
9497
}
9598
}
9699
};
97100
p_MPs->parallel_for(setMPStress, "setMPStress");
101+
int_xx++;
98102
}
99103

100104
void MPMesh::calculateStressDivergence(){
@@ -170,19 +174,20 @@ void MPMesh::calculateStressDivergence(){
170174
p_MPs->parallel_for(stress_div, "assembly");
171175

172176
//TODO COMMUNICATE THE VERTEX FIELDS
173-
174177

175178
auto stressDivergence = p_mesh->getMeshField<MeshF_StressDivergence>();
176179

177-
Kokkos::parallel_for("calculate_divergence", numVtx, KOKKOS_LAMBDA(const int vtx){
178-
stressDivergence(vtx, 0) = stress_divU(vtx);
179-
stressDivergence(vtx, 1) = stress_divV(vtx);
180+
Kokkos::parallel_for("calculate_divergence", numVtx, KOKKOS_LAMBDA(const int vtx){
181+
182+
double ramp = nearAnEdge_l(vtx);
183+
double invM = 1.0/vtxMatrixMass_l(vtx);
184+
invM = vtxMatrixMass_l(vtx) >1e-4 ? invM : 0;
185+
186+
stressDivergence(vtx, 0) = ramp * stress_divU(vtx) + (1 - ramp) * divU_edge(vtx) * invM ;
187+
stressDivergence(vtx, 1) = ramp * stress_divV(vtx) + (1 - ramp) * divV_edge(vtx) * invM ;
180188
//Debugging
181-
if (vtx >= 10 && vtx <= 11) {
182-
//printf("Vtx %d Divergence %.15e %.15e %.15e %.15e %.15e %.15e \n", vtx, nearAnEdge_l(vtx), vtxMatrixMass_l(vtx),
183-
// stress_divU(vtx), stress_divV(vtx),
184-
// divU_edge(vtx), divV_edge(vtx));
185-
printf("Vtx %d Divergence %.15e %.15e \n", vtx, stressDivergence(vtx, 0), stressDivergence(vtx, 1));
189+
if (vtx == 10 || vtx == 11 || vtx == 1476 || vtx == 1481) {
190+
//printf("Vtx %d Divergence %.15e %.15e \n", vtx, stressDivergence(vtx, 0), stressDivergence(vtx, 1));
186191
}
187192
});
188193

@@ -544,7 +549,7 @@ void MPMesh::push_ahead(){
544549
calcBasis();
545550
sphericalInterpolation<MeshF_RotLatLonIncr>(*this);
546551
sphericalInterpolation<MeshF_OnSurfVeloIncr>(*this);
547-
552+
sphericalInterpolation2Fields(*this);
548553
//Push the MPs
549554
p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius(), p_mesh->getRotatedFlag());
550555
pumipic::RecordTime("PolyMPO_interpolateAndPush", timer.seconds());

‎src/pmpo_materialPoints.hpp‎

Lines changed: 19 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -43,7 +43,8 @@ enum MaterialPointSlice {
4343
MPF_Tgt_Proc_ID,
4444
MPF_Area,
4545
MPF_IcePressure,
46-
MPF_ReplacementPressure
46+
MPF_ReplacementPressure,
47+
MPF_Vel_IncrTimesTanLatVertexOverRadius
4748
};
4849

4950
enum Operating_Mode{
@@ -75,6 +76,7 @@ template <> struct mpSliceToMeshField < MPF_Tgt_Proc_ID > { using type =
7576
template <> struct mpSliceToMeshField < MPF_Area > { using type = doubleSclr_t; };
7677
template <> struct mpSliceToMeshField < MPF_IcePressure > { using type = doubleSclr_t; };
7778
template <> struct mpSliceToMeshField < MPF_ReplacementPressure > { using type = doubleSclr_t; };
79+
template <> struct mpSliceToMeshField < MPF_Vel_IncrTimesTanLatVertexOverRadius> { using type = vec2d_t; };
7880

7981
template <MaterialPointSlice slice>
8082
static constexpr int mpSliceToNumEntries() {
@@ -113,7 +115,8 @@ typedef MemberTypes<mpSliceToMeshField < MPF_Status >::type,
113115
mpSliceToMeshField < MPF_Tgt_Proc_ID >::type,
114116
mpSliceToMeshField < MPF_Area >::type,
115117
mpSliceToMeshField < MPF_IcePressure >::type,
116-
mpSliceToMeshField < MPF_ReplacementPressure >::type
118+
mpSliceToMeshField < MPF_ReplacementPressure >::type,
119+
mpSliceToMeshField < MPF_Vel_IncrTimesTanLatVertexOverRadius >::type
117120
>MaterialPointTypes;
118121
typedef ps::ParticleStructure<MaterialPointTypes> PS;
119122

@@ -230,7 +233,10 @@ class MaterialPoints {
230233
//Velocity
231234
auto velMPs = MPs->get<MPF_Vel>();
232235
auto velIncr = MPs->get<MPF_Vel_Incr>();
233-
236+
//Stress
237+
auto MPsStress = MPs->get<MPF_Stress>();
238+
auto MPsStressMetric = MPs->get<MPF_Vel_IncrTimesTanLatVertexOverRadius>();
239+
234240
auto mpAppID = MPs->get<MPF_MP_APP_ID>();
235241

236242
auto updateRotLatLon = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
@@ -252,6 +258,16 @@ class MaterialPoints {
252258
tgtPosXYZ(mp,2) = radius * std::sin(geoLat);
253259
velMPs(mp,0) = velMPs(mp,0) + velIncr(mp,0);
254260
velMPs(mp,1) = velMPs(mp,1) + velIncr(mp,1);
261+
//Stress MP term
262+
Vec3d stress_prev(MPsStress(mp,0), MPsStress(mp,1), MPsStress(mp,2));
263+
MPsStress(mp,0) = stress_prev[0] * (1 - MPsStressMetric(mp, 1)) + 2 * MPsStressMetric(mp, 1) * stress_prev[2];
264+
MPsStress(mp,1) = stress_prev[1] + 2*MPsStressMetric(mp, 0) * stress_prev[2];
265+
MPsStress(mp,2) = stress_prev[2] + MPsStressMetric(mp, 0) * (stress_prev[0] - stress_prev[1]);
266+
267+
if(mpAppID(mp)==0){
268+
// printf("Correction: %.15e %.15e\n", MPsStressMetric(mp, 0), MPsStressMetric(mp, 1));
269+
// printf("Stress in GPU After Correction: %.15e %.15e %.15e\n", MPsStress(mp, 0), MPsStress(mp, 1), MPsStress(mp, 2));
270+
}
255271
}
256272
};
257273
ps::parallel_for(MPs, updateRotLatLon,"updateRotationalLatitudeLongitude");

‎src/pmpo_wachspressBasis.hpp‎

Lines changed: 40 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -15,8 +15,6 @@ void sphericalInterpolation(MPMesh& mpMesh){
1515
auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
1616
int numVtxs = p_mesh->getNumVertices();
1717
auto elm2VtxConn = p_mesh->getElm2VtxConn();
18-
double radius = p_mesh->getSphereRadius();
19-
PMT_ALWAYS_ASSERT(radius >0);
2018

2119
auto p_MPs = mpMesh.p_MPs;
2220
auto MPsPosition = p_MPs->getPositions();
@@ -43,6 +41,46 @@ void sphericalInterpolation(MPMesh& mpMesh){
4341
pumipic::RecordTime("PolyMPO_sphericalInterpolation", timer.seconds());
4442
}
4543

44+
inline void sphericalInterpolation2Fields(MPMesh& mpMesh){
45+
46+
Kokkos::Timer timer;
47+
48+
auto p_mesh = mpMesh.p_mesh;
49+
auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
50+
int numVtxs = p_mesh->getNumVertices();
51+
auto elm2VtxConn = p_mesh->getElm2VtxConn();
52+
53+
//Material Points Data
54+
auto p_MPs = mpMesh.p_MPs;
55+
auto MPsPosition = p_MPs->getPositions();
56+
auto MPsBasis = p_MPs->getData<MPF_Basis_Vals>();
57+
58+
constexpr MaterialPointSlice mpfIndex = MPF_Vel_IncrTimesTanLatVertexOverRadius;
59+
auto mpField = p_MPs->getData<mpfIndex>();
60+
const int numEntries = mpSliceToNumEntries<mpfIndex>();
61+
62+
constexpr MeshFieldIndex mfIndex1 = MeshF_OnSurfVeloIncr;
63+
auto meshField1 = p_mesh->getMeshField<mfIndex1>();
64+
constexpr MeshFieldIndex mfIndex2 = MeshF_TanLatVertexRotatedOverRadius;
65+
auto meshField2 = p_mesh->getMeshField<mfIndex2>();
66+
67+
auto interpolation2 = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
68+
if(mask) { //if material point is 'active'/'enabled'
69+
int numVtx = elm2VtxConn(elm,0);
70+
for(int entry=0; entry<numEntries; entry++){
71+
double mpValue = 0.0;
72+
for(int i=1; i<= numVtx; i++)
73+
mpValue += meshField1(elm2VtxConn(elm,i)-1, entry) * meshField2(elm2VtxConn(elm,i)-1, 0) * MPsBasis(mp,i-1);
74+
mpField(mp,entry) = mpValue;
75+
}
76+
}
77+
};
78+
p_MPs->parallel_for(interpolation2, "interpolation");
79+
80+
pumipic::RecordTime("PolyMPO_sphericalInterpolation2Fields", timer.seconds());
81+
}
82+
83+
4684
KOKKOS_INLINE_FUNCTION
4785
void wachpress_weights_grads_2D(int numVtx, const Kokkos::View<double[maxVtxsPerElm][2],
4886
Kokkos::LayoutStride, Kokkos::MemoryTraits<Kokkos::Unmanaged>>& gnom_vtx_subview,

0 commit comments

Comments
 (0)