Skip to content

Commit cf5afec

Browse files
dhyan1272dhyan1272
committed
Stress Divergence calculation updated
1 parent 8c917ff commit cf5afec

9 files changed

Lines changed: 89 additions & 28 deletions

‎src/pmpo_MPMesh.cpp‎

Lines changed: 23 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -92,6 +92,7 @@ void MPMesh::calculateStress(){
9292
MPsStress(mp, m) = stress[m];
9393
//Debugging
9494
if(MPsAppID(mp)==0){
95+
printf("Strain in GPU: %.15e %.15e %.15e\n", MPsStrainRate(mp, 0), MPsStrainRate(mp, 1), MPsStrainRate(mp, 2));
9596
printf("Stress in GPU: %.15e %.15e %.15e\n", MPsStress(mp, 0), MPsStress(mp, 1), MPsStress(mp, 2));
9697
}
9798
}
@@ -107,7 +108,8 @@ void MPMesh::calculateStressDivergence(){
107108
auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
108109
int numVertices = p_mesh->getNumVertices();
109110
auto tanLatVertexRotatedOverRadius = p_mesh->getMeshField<MeshF_TanLatVertexRotatedOverRadius>();
110-
111+
auto interiorVertex = p_mesh->getMeshField<MeshF_InteriorVertex>();
112+
111113
//Material Points
112114
auto MPsAppID = p_MPs->getData<MPF_MP_APP_ID>();
113115
auto weight = p_MPs->getData<MPF_Basis_Vals>();
@@ -116,20 +118,22 @@ void MPMesh::calculateStressDivergence(){
116118
auto MPsStress = p_MPs->getData<MPF_Stress>();
117119
auto mpPositions = p_MPs->getData<MPF_Cur_Pos_XYZ>();
118120

119-
auto VtxCoeffs_new = this->precomputedVtxCoeffs_new;
121+
auto VtxCoeffs_new = this->precomputedVtxCoeffs_new;
122+
auto vtxMatrixMass_l = this->vtxMatrixMass;
123+
auto nearAnEdge_l = this->nearAnEdge;
120124

121125
//Earth Radius
122126
double radius = 1.0;
123127
if(p_mesh->getGeomType() == geom_spherical_surf)
124128
radius=p_mesh->getSphereRadius();
125129

126130
//Reconstructed the stress
127-
Kokkos::View<double*[3]> stress_rec("stress_rec", p_mesh->getNumVertices());
128131
Kokkos::View<double*> stress_divU("stress_divu", p_mesh->getNumVertices());
129132
Kokkos::View<double*> stress_divV("stress_divv", p_mesh->getNumVertices());
130-
Kokkos::View<double*> dSdX("dSdX", p_mesh->getNumVertices());
131-
Kokkos::View<double*> dSdY("dSdY", p_mesh->getNumVertices());
132-
133+
134+
Kokkos::View<double*> divU_edge("divUedge", p_mesh->getNumVertices());
135+
Kokkos::View<double*> divV_edge("divVedge", p_mesh->getNumVertices());
136+
133137
//Assemble fields for Stress Divergence
134138
auto stress_div = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
135139
if(mask) { //if material point is 'active'/'enabled'
@@ -152,16 +156,17 @@ void MPMesh::calculateStressDivergence(){
152156
VtxCoeffs_new(vID,2, 2)*CoordDiffs[2] +
153157
VtxCoeffs_new(vID,2, 3)*CoordDiffs[3]);
154158

155-
for (int k=0; k<3; k++){
156-
auto val = factor*MPsStress(mp,k);
157-
Kokkos::atomic_add(&stress_rec(vID,k), val);
158-
}
159159
Kokkos::atomic_add(&stress_divU(vID), factor1 * MPsStress(mp, 0) + factor2 * MPsStress(mp, 2) -
160160
2 * tanLatVertexRotatedOverRadius(vID, 0) * factor * MPsStress(mp, 2));
161161
Kokkos::atomic_add(&stress_divV(vID), factor2 * MPsStress(mp, 1) + factor1*MPsStress(mp, 2) +
162162
factor * tanLatVertexRotatedOverRadius(vID, 0) * (MPsStress(mp, 0)-MPsStress(mp, 1)));
163-
Kokkos::atomic_add(&dSdX(vID), -factor * weight_grads(mp, i*2 + 0));
164-
Kokkos::atomic_add(&dSdY(vID), -factor * weight_grads(mp, i*2 + 1));
163+
164+
165+
Kokkos::atomic_add(&divU_edge(vID), - weight_grads(mp, i*2 + 0) * MPsStress(mp, 0) - weight_grads(mp, i*2 + 1) * MPsStress(mp, 2) -
166+
2 * tanLatVertexRotatedOverRadius(vID, 0) * w_vtx * MPsStress(mp, 2));
167+
Kokkos::atomic_add(&divV_edge(vID), - weight_grads(mp, i*2 + 1) * MPsStress(mp, 1) - weight_grads(mp, i*2 + 0) * MPsStress(mp, 2) +
168+
w_vtx * tanLatVertexRotatedOverRadius(vID, 0) * (MPsStress(mp, 0)- MPsStress(mp, 1)));
169+
165170
}
166171
}
167172
};
@@ -171,22 +176,17 @@ void MPMesh::calculateStressDivergence(){
171176
//COMMUNICATE THE VERTEX FIELDS
172177

173178
//TODO put as mesh field
174-
Kokkos::View<doubleSclr_t*> stressDivergence("stressDivergence", p_mesh->getNumVertices());
175-
auto areaVertex = p_mesh->getMeshField<MeshF_DualTriangleArea>();
179+
Kokkos::View<vec2d_t*> stressDivergence("stressDivergence", p_mesh->getNumVertices());
176180

177181
Kokkos::parallel_for("calculate_divergence", numVtx, KOKKOS_LAMBDA(const int vtx){
178-
double threshold = 0.15 / Kokkos::sqrt(areaVertex(vtx, 0));
179-
double valX = Kokkos::max(Kokkos::abs(dSdX(vtx)) - threshold, 0.0);
180-
double dSdX_filtered = Kokkos::copysign(valX, dSdX(vtx));
181-
double valY = Kokkos::max(Kokkos::abs(dSdY(vtx)) - threshold, 0.0);
182-
double dSdY_filtered = Kokkos::copysign(valY, dSdY(vtx));
183182

184-
stressDivergence(vtx, 0) = stress_divU(vtx) + dSdX_filtered * stress_rec(vtx, 0) + dSdY_filtered * stress_rec(vtx, 2);
185-
stressDivergence(vtx, 1) = stress_divV(vtx) + dSdY_filtered * stress_rec(vtx, 1) + dSdX_filtered * stress_rec(vtx, 2);
183+
stressDivergence(vtx, 0) = stress_divU(vtx);
184+
stressDivergence(vtx, 1) = stress_divV(vtx);
186185
//Debugging
187186
if (vtx >= 10 && vtx <= 11) {
188-
printf("Vtx %d Area %.15e ds: %.15e %.15e \n", vtx, areaVertex(vtx, 0), dSdX_filtered, dSdY_filtered);
189-
printf("Vtx %d Divergence %.15e %.15e \n", vtx, stressDivergence(vtx, 0), stressDivergence(vtx, 1));
187+
printf("Vtx %d Divergence %.15e %.15e %.15e %.15e %.15e %.15e \n", vtx, nearAnEdge_l(vtx), vtxMatrixMass_l(vtx),
188+
stress_divU(vtx), stress_divV(vtx),
189+
divU_edge(vtx), divV_edge(vtx));
190190
}
191191
});
192192

‎src/pmpo_MPMesh.hpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -72,6 +72,8 @@ class MPMesh{
7272
void reconstruct_coeff_full();
7373
void invertMatrix(const Kokkos::View<double**>& vtxMatrices, const double& radius);
7474
Kokkos::View<double*[vec3d_nEntries][vec4d_nEntries]> precomputedVtxCoeffs_new;
75+
Kokkos::View<double*> nearAnEdge;
76+
Kokkos::View<double*> vtxMatrixMass;
7577

7678
//Not used currently
7779
std::map<MeshFieldIndex, std::function<void()>> reconstructSlice = std::map<MeshFieldIndex, std::function<void()>>();

‎src/pmpo_MPMesh_assembly.hpp‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -165,22 +165,32 @@ void MPMesh::reconstruct_coeff_full(){
165165
communicate_and_take_halo_contributions(vtxMatrices, numVertices, numEntriesMatrix, mode, op);
166166
}
167167
pumipic::RecordTime("Communicate Matrix Values" + std::to_string(self), timer.seconds());
168+
169+
//Stroe the 1st matrix element
170+
Kokkos::View<double*>vtxMatrixMass_l("vtxMass", numVertices);
171+
Kokkos::parallel_for("storeMatrixMass", numVertices, KOKKOS_LAMBDA(const int vtx){
172+
vtxMatrixMass_l(vtx) = vtxMatrices(vtx, 0);
173+
});
174+
this->vtxMatrixMass = vtxMatrixMass_l;
168175

169176
invertMatrix(vtxMatrices, radius);
170177
}
171178

172179
void MPMesh::invertMatrix(const Kokkos::View<double**>& vtxMatrices, const double& radius){
173180
std::cout<<__FUNCTION__<<std::endl;
181+
174182
int nVertices = p_mesh->getNumVertices();
175183
auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
176184
auto dual_triangle_area = p_mesh->getMeshField<MeshF_DualTriangleArea>();
185+
auto interiorVertex = p_mesh->getMeshField<MeshF_InteriorVertex>();
177186
bool isRotated = p_mesh->getRotatedFlag();
178187

179188
double eps = 1e-7;
180189
double truncateFactor = 0.05;
181190

182191
Kokkos::View<double*[3][vec4d_nEntries]> VtxCoeffs("VtxCoeffs", nVertices);
183192
Kokkos::deep_copy(VtxCoeffs, 0.0);
193+
Kokkos::View<double*> nearAnEdge_l("nearAnEdge_l", nVertices);
184194

185195
Kokkos::parallel_for("invertMatrix", nVertices, KOKKOS_LAMBDA(const int vtx){
186196
if(vtxMatrices(vtx, 0) < eps)
@@ -287,8 +297,25 @@ void MPMesh::invertMatrix(const Kokkos::View<double**>& vtxMatrices, const doubl
287297
VtxCoeffs(vtx, 2, 1) = invM2D[1] * rotateScaleM(0, 0) + invM2D[3] * rotateScaleM(0, 1) + invM2D[4] * rotateScaleM(0, 2);
288298
VtxCoeffs(vtx, 2, 2) = invM2D[1] * rotateScaleM(1, 0) + invM2D[3] * rotateScaleM(1, 1) + invM2D[4] * rotateScaleM(1, 2);
289299
VtxCoeffs(vtx, 2, 3) = invM2D[1] * rotateScaleM(2, 0) + invM2D[3] * rotateScaleM(2, 1) + invM2D[4] * rotateScaleM(2, 2);
300+
301+
//Calculate a smmooth flag for if we are near a material edge (from MPAS)
302+
double pMassGradNorm = vtxMatrices(vtx, 1) * vtxMatrices(vtx, 1) + vtxMatrices(vtx, 2) * vtxMatrices(vtx, 2) +
303+
vtxMatrices(vtx, 3) * vtxMatrices(vtx, 3);
304+
pMassGradNorm = (sqrt(pMassGradNorm) / vtx_area_sqrt) / Kokkos::max(vtxMatrices(vtx, 0), 1e-4);
305+
306+
double ramp = 2.2 - 10.0 * pMassGradNorm;
307+
ramp = ramp < 0.0 ? 0.0 : ramp;
308+
ramp = ramp > 1.0 ? 1.0 : ramp;
309+
310+
double massRamp = 4.0 * (vtxMatrices(vtx, 0) - 1.0);
311+
massRamp = massRamp < 0.0 ? 0.0 : massRamp;
312+
massRamp = massRamp > 1.0 ? 1.0 : massRamp;
313+
314+
ramp *= massRamp;
315+
nearAnEdge_l(vtx) = (interiorVertex(vtx) == 1) ? ramp : 0;
290316
});
291317
this->precomputedVtxCoeffs_new = VtxCoeffs;
318+
this->nearAnEdge = nearAnEdge_l;
292319
}
293320

294321
template <MeshFieldIndex meshFieldIndex>

‎src/pmpo_c.cpp‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -927,6 +927,20 @@ void polympo_setVtxGlobal_f(MPMesh_ptr p_mpmesh, const int nVertices, const int*
927927
p_mesh->setVtxGlobal(vtxGlobal);
928928
}
929929

930+
void polympo_setInteriorVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array){
931+
checkMPMeshValid(p_mpmesh);
932+
auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
933+
PMT_ALWAYS_ASSERT(nVertices==p_mesh->getNumVertices());
934+
Kokkos::View<int*, Kokkos::HostSpace> arrayHost("arrayHost", nVertices);
935+
for (int i = 0; i < nVertices; i++) {
936+
arrayHost(i) = array[i];
937+
}
938+
939+
auto vtxInterior = p_mesh->getMeshField<polyMPO::MeshF_InteriorVertex>();
940+
Kokkos::deep_copy(vtxInterior, arrayHost);
941+
}
942+
943+
930944
//Mesh Fields
931945
int polympo_getMeshFVtxType_f() {
932946
return polyMPO::MeshFType_VtxBased;

‎src/pmpo_c.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -74,6 +74,7 @@ void polympo_setOwningProc_f(MPMesh_ptr p_mpmesh, const int nCells, const int* a
7474
void polympo_setOwningProcVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array);
7575
void polympo_setElmGlobal_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array);
7676
void polympo_setVtxGlobal_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array);
77+
void polympo_setInteriorVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array);
7778

7879
//Mesh fields
7980
int polympo_getMeshFVtxType_f();

‎src/pmpo_fortran.f90‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -632,6 +632,15 @@ subroutine polympo_setVtxGlobal(mpMesh, nVertices, array) &
632632
type(c_ptr), intent(in), value :: array
633633
end subroutine
634634

635+
subroutine polympo_setInteriorVertex(mpMesh, nVertices, array) &
636+
bind(C, NAME='polympo_setInteriorVertex_f')
637+
use :: iso_c_binding
638+
type(c_ptr), value :: mpMesh
639+
integer(c_int), value :: nVertices
640+
type(c_ptr), intent(in), value :: array
641+
end subroutine
642+
643+
635644
!---------------------------------------------------------------------------
636645
!> @brief get enum for vertex mesh fields
637646
!---------------------------------------------------------------------------
@@ -914,6 +923,7 @@ subroutine polympo_getMeshDualTriangleArea(mpMesh, nVertices, areaTriangle) &
914923
type(c_ptr), value :: areaTriangle
915924
end subroutine
916925

926+
917927
subroutine polympo_setGnomonicProjection(mpMesh) &
918928
bind(C, NAME='polympo_setGnomonicProjection_f')
919929
use :: iso_c_binding

‎src/pmpo_materialPoints.hpp‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -37,7 +37,6 @@ enum MaterialPointSlice {
3737
MPF_Vel_Incr,
3838
MPF_Strain_Rate,
3939
MPF_Stress,
40-
MPF_Stress_Div,
4140
MPF_Shear_Traction,
4241
MPF_Constv_Mdl_Param,
4342
MPF_MP_APP_ID,
@@ -69,7 +68,6 @@ template <> struct mpSliceToMeshField < MPF_Rot_Lat_Lon_Incr > { using type =
6968
template <> struct mpSliceToMeshField < MPF_Vel_Incr > { using type = vec2d_t; };
7069
template <> struct mpSliceToMeshField < MPF_Strain_Rate > { using type = double[3]; };
7170
template <> struct mpSliceToMeshField < MPF_Stress > { using type = double[3]; };
72-
template <> struct mpSliceToMeshField < MPF_Stress_Div > { using type = vec3d_t; };
7371
template <> struct mpSliceToMeshField < MPF_Shear_Traction > { using type = vec3d_t; };
7472
template <> struct mpSliceToMeshField < MPF_Constv_Mdl_Param > { using type = double[12]; };
7573
template <> struct mpSliceToMeshField < MPF_MP_APP_ID > { using type = int; };
@@ -109,7 +107,6 @@ typedef MemberTypes<mpSliceToMeshField < MPF_Status >::type,
109107
mpSliceToMeshField < MPF_Vel_Incr >::type,
110108
mpSliceToMeshField < MPF_Strain_Rate >::type,
111109
mpSliceToMeshField < MPF_Stress >::type,
112-
mpSliceToMeshField < MPF_Stress_Div >::type,
113110
mpSliceToMeshField < MPF_Shear_Traction >::type,
114111
mpSliceToMeshField < MPF_Constv_Mdl_Param >::type,
115112
mpSliceToMeshField < MPF_MP_APP_ID >::type,

‎src/pmpo_mesh.cpp‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -47,7 +47,10 @@ namespace polyMPO{
4747
auto tanLatVertexRotOverRadiusEntry = meshFields2TypeAndString.at(MeshF_TanLatVertexRotatedOverRadius);
4848
PMT_ALWAYS_ASSERT(tanLatVertexRotOverRadiusEntry.first == MeshFType_VtxBased);
4949
tanLatVertexRotatedOverRadius_ = MeshFView<MeshF_TanLatVertexRotatedOverRadius>(tanLatVertexRotOverRadiusEntry.second, numVtxs_);
50-
50+
51+
auto interiorVertexEntry = meshFields2TypeAndString.at(MeshF_InteriorVertex);
52+
PMT_ALWAYS_ASSERT(interiorVertexEntry.first == MeshFType_VtxBased);
53+
interiorVertex_ = MeshFView<MeshF_InteriorVertex>(interiorVertexEntry.second, numVtxs_);
5154
}
5255

5356
void Mesh::setMeshElmBasedFieldSize(){

‎src/pmpo_mesh.hpp‎

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -30,7 +30,8 @@ enum MeshFieldIndex{
3030
MeshF_VtxGnomProj,
3131
MeshF_ElmCenterGnomProj,
3232
MeshF_TanLatVertexRotatedOverRadius,
33-
MeshF_SolveStress
33+
MeshF_SolveStress,
34+
MeshF_InteriorVertex
3435
};
3536
enum MeshFieldType{
3637
MeshFType_Invalid = -2,
@@ -54,6 +55,7 @@ template <> struct meshFieldToType < MeshF_VtxGnomProj > { using type = Ko
5455
template <> struct meshFieldToType < MeshF_ElmCenterGnomProj > { using type = Kokkos::View<double*[4]>; };
5556
template <> struct meshFieldToType < MeshF_TanLatVertexRotatedOverRadius > { using type = Kokkos::View<doubleSclr_t*>; };
5657
template <> struct meshFieldToType < MeshF_SolveStress > { using type = IntView; };
58+
template <> struct meshFieldToType < MeshF_InteriorVertex > { using type = IntView; };
5759

5860
template <MeshFieldIndex index>
5961
using MeshFView = typename meshFieldToType<index>::type;
@@ -75,6 +77,7 @@ const std::map<MeshFieldIndex, std::pair<MeshFieldType, std::string>> meshFields
7577
{MeshF_ElmCenterGnomProj,{MeshFType_ElmBased,"MeshField_ElementCenterGnomonicprojection"}},
7678
{MeshF_TanLatVertexRotatedOverRadius, {MeshFType_VtxBased,"MeshField_TanLatVertexRotatedOverRadius"}},
7779
{MeshF_SolveStress, {MeshFType_ElmBased,"MeshField_SolveStress"}},
80+
{MeshF_InteriorVertex, {MeshFType_VtxBased,"MeshField_InteriorVertex"}}
7881
};
7982

8083
enum mesh_type {mesh_unrecognized_lower = -1,
@@ -120,6 +123,7 @@ class Mesh {
120123
//DoubleMat2DView vtxStress_;
121124
MeshFView<MeshF_TanLatVertexRotatedOverRadius> tanLatVertexRotatedOverRadius_;
122125
MeshFView<MeshF_SolveStress> solveStress_;
126+
MeshFView<MeshF_InteriorVertex> interiorVertex_;
123127
bool isRotatedFlag = false;
124128
double elasticTimeStep_;
125129
double dynamicTimeStep_;
@@ -274,6 +278,9 @@ auto Mesh::getMeshField(){
274278
else if constexpr (index==MeshF_SolveStress){
275279
return solveStress_;
276280
}
281+
else if constexpr (index==MeshF_InteriorVertex){
282+
return interiorVertex_;
283+
}
277284
fprintf(stderr,"Mesh Field Index error!\n");
278285
exit(1);
279286
}

0 commit comments

Comments
 (0)