Skip to content

Commit 5ff1728

Browse files
dhyan1272dhyan1272
committed
StrainRate Stress
1 parent 0ca0073 commit 5ff1728

11 files changed

Lines changed: 311 additions & 19 deletions

‎src/CMakeLists.txt‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@ set(HEADERS
77
pmpo_MPMesh_assembly.hpp
88
pmpo_c.h
99
pmpo_createTestMPMesh.hpp
10+
pmpo_const_relation.hpp
1011
)
1112

1213
set(SOURCES
@@ -53,4 +54,4 @@ add_library(polyMPO INTERFACE)
5354
target_link_libraries(polyMPO INTERFACE ${polyMPO_EXPORTED_TARGETS})
5455
bob_export_target(polyMPO)
5556

56-
bob_end_subdir()
57+
bob_end_subdir()

‎src/pmpo_MPMesh.cpp‎

Lines changed: 36 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22
#include "pmpo_utils.hpp"
33
#include "pmpo_MPMesh.hpp"
44
#include "pmpo_wachspressBasis.hpp"
5+
#include "pmpo_const_relation.hpp"
56

67
namespace polyMPO{
78

@@ -11,7 +12,7 @@ void MPMesh::calculateStrain(){
1112
auto MPsPosition = p_MPs->getPositions();
1213
auto MPsBasis = p_MPs->getData<MPF_Basis_Vals>();
1314
auto MPsAppID = p_MPs->getData<MPF_MP_APP_ID>();
14-
15+
auto MPsStrainRate = p_MPs->getData<MPF_Strain_Rate>();
1516
//Mesh Fields
1617
auto tanLatVertexRotatedOverRadius = p_mesh->getMeshField<MeshF_TanLatVertexRotatedOverRadius>();
1718
auto gnomProjVtx = p_mesh->getMeshField<MeshF_VtxGnomProj>();
@@ -59,12 +60,46 @@ void MPMesh::calculateStrain(){
5960
v22 = v22 + gradBasisByArea[i*2 + 1] * velField(iVertex, 1);
6061
uTanOverR = uTanOverR + basisByArea[i] * tanLatVertexRotatedOverRadius(iVertex, 0) * velField(iVertex, 0);
6162
vTanOverR = vTanOverR + basisByArea[i] * tanLatVertexRotatedOverRadius(iVertex, 0) * velField(iVertex, 1);
63+
MPsStrainRate(mp, 0) = v11 - vTanOverR;
64+
MPsStrainRate(mp, 1) = v22;
65+
MPsStrainRate(mp, 2) = 0.5*(v12 + v21 + uTanOverR);
6266
}
67+
MPsStrainRate(mp, 0) = v11 - vTanOverR;
68+
MPsStrainRate(mp, 1) = v22;
69+
MPsStrainRate(mp, 2) = 0.5*(v12 + v21 + uTanOverR);
6370
}
6471
};
6572
p_MPs->parallel_for(setMPStrainRate, "setMPStrainRate");
6673
}
6774

75+
void MPMesh::calculateStress(){
76+
//MeshFields
77+
auto solveStress = p_mesh->getMeshField<polyMPO::MeshF_SolveStress>();
78+
auto elasticTimeStep = p_mesh->getElasticTimeStep();
79+
auto dynamicTimeStep = p_mesh->getDynamicTimeStep();
80+
auto dampingTimescale = polyMPO::dampingTimescaleParameter * dynamicTimeStep;
81+
//MPFields
82+
auto MPsAppID = p_MPs->getData<MPF_MP_APP_ID>();
83+
auto MPsStrainRate = p_MPs->getData<MPF_Strain_Rate>();
84+
auto MPsArea = p_MPs->getData<polyMPO::MPF_Area>();
85+
auto MPsIcePressure = p_MPs->getData<polyMPO::MPF_IcePressure>();
86+
auto MPsRepPressure = p_MPs->getData<polyMPO::MPF_ReplacementPressure>();
87+
88+
auto setMPStress = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
89+
if(mask){
90+
//Debugging
91+
if(MPsAppID(mp)==0){
92+
printf("Strain %.15e %.15e %.15e\n", MPsStrainRate(mp, 0), MPsStrainRate(mp, 1), MPsStrainRate(mp, 2));
93+
printf("Mesh:%.15e %d, MP:%.15e %.15e %.15e\n",elasticTimeStep,solveStress(elm),MPsArea(mp,0),MPsIcePressure(mp,0),MPsRepPressure(mp,0));
94+
Vec3d strain_rate (MPsStrainRate(mp, 0), MPsStrainRate(mp, 1), MPsStrainRate(mp, 2));
95+
Vec3d stress(0, 0, 0);
96+
constitutive_evp(strain_rate, stress, MPsIcePressure(mp, 0), MPsRepPressure(mp, 0), MPsArea(mp, 0), elasticTimeStep, dampingTimescale);
97+
}
98+
}
99+
};
100+
p_MPs->parallel_for(setMPStress, "setMPStress");
101+
}
102+
68103
void MPMesh::calcBasis() {
69104
assert(p_mesh->getGeomType() == geom_spherical_surf);
70105

‎src/pmpo_MPMesh.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -87,6 +87,7 @@ class MPMesh{
8787

8888
void printVTP_mesh(int printVTPIndex);
8989
void calculateStrain();
90+
void calculateStress();
9091
};
9192

9293
}//namespace polyMPO end

‎src/pmpo_MPMesh_assembly.hpp‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -345,6 +345,15 @@ void MPMesh::assemblyVtx1(){
345345
if(numProcsTot>1)
346346
communicate_and_take_halo_contributions(meshField, numVertices, numEntries, 0, 0);
347347
pumipic::RecordTime("Communicate Field Values" + std::to_string(self), timer.seconds());
348+
349+
Kokkos::parallel_for("printSymmetricBlock", numVertices, KOKKOS_LAMBDA(const int vtx){
350+
if (vtx >= 10 && vtx <= 10) {
351+
printf("Field in %d: ", vtx);
352+
for (int k=0; k<numEntries; k++)
353+
printf(" %.15e ", meshField(vtx, k));
354+
printf("\n");
355+
}
356+
});
348357
}
349358

350359
//Start Communication routine

‎src/pmpo_c.cpp‎

Lines changed: 118 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -624,14 +624,10 @@ void polympo_getMPStrainRate_f(MPMesh_ptr p_mpmesh, const int nComps, const int
624624
(void)mpStrainRateHost;
625625
}
626626

627-
void polympo_setMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStressIn){
627+
void polympo_setMPStress_f(MPMesh_ptr p_mpmesh){
628628
checkMPMeshValid(p_mpmesh);
629-
std::cerr << "Error: This routine is not implemented yet\n";
630-
exit(1);
631-
(void)p_mpmesh;
632-
(void)nComps;
633-
(void)numMPs;
634-
(void)mpStressIn;
629+
auto mpMesh = ((polyMPO::MPMesh*)p_mpmesh);
630+
mpMesh->calculateStress();
635631
}
636632

637633
void polympo_getMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpStressHost){
@@ -644,6 +640,93 @@ void polympo_getMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numM
644640
(void)mpStressHost;
645641
}
646642

643+
void polympo_setAreaMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* areaMPHost){
644+
Kokkos::Timer timer;
645+
checkMPMeshValid(p_mpmesh);
646+
647+
auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
648+
//Rank information
649+
int self;
650+
MPI_Comm comm = p_MPs->getMPIComm();
651+
MPI_Comm_rank(comm, &self);
652+
//Asserts
653+
PMT_ALWAYS_ASSERT(nComps == 1);
654+
PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
655+
//MP Data
656+
auto mpArea = p_MPs->getData<polyMPO::MPF_Area>();
657+
auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
658+
//Copy to device
659+
kkViewHostU<const double**> mpAreaIn_h(areaMPHost, nComps, numMPs);
660+
Kokkos::View<double**> mpAreaIn_d("mpAreaDevice", nComps, numMPs);
661+
Kokkos::deep_copy(mpAreaIn_d, mpAreaIn_h);
662+
//Set in PS
663+
auto setMPArea = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
664+
if(mask){
665+
mpArea(mp,0) = mpAreaIn_d(0, mpAppID(mp));
666+
}
667+
};
668+
p_MPs->parallel_for(setMPArea, "setMPArea");
669+
pumipic::RecordTime("PolyMPO_setMPArea" + std::to_string(self), timer.seconds());
670+
}
671+
672+
void polympo_setIcePressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* icePressureMPHost){
673+
Kokkos::Timer timer;
674+
checkMPMeshValid(p_mpmesh);
675+
676+
auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
677+
//Rank information
678+
int self;
679+
MPI_Comm comm = p_MPs->getMPIComm();
680+
MPI_Comm_rank(comm, &self);
681+
//Asserts
682+
PMT_ALWAYS_ASSERT(nComps == 1);
683+
PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
684+
//MP Data
685+
auto mpIcePressure = p_MPs->getData<polyMPO::MPF_IcePressure>();
686+
auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
687+
//Copy to device
688+
kkViewHostU<const double**> mpIcePressure_h(icePressureMPHost, nComps, numMPs);
689+
Kokkos::View<double**> mpIcePressure_d("mpIcePressureDevice", nComps, numMPs);
690+
Kokkos::deep_copy(mpIcePressure_d, mpIcePressure_h);
691+
//Set in PS
692+
auto setMPIcePressure = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
693+
if(mask){
694+
mpIcePressure(mp,0) = mpIcePressure_d(0, mpAppID(mp));
695+
}
696+
};
697+
p_MPs->parallel_for(setMPIcePressure, "setIcePressure");
698+
pumipic::RecordTime("PolyMPO_setIcePressure" + std::to_string(self), timer.seconds());
699+
}
700+
701+
void polympo_setReplacementPressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* replacementPressureMPHost){
702+
Kokkos::Timer timer;
703+
checkMPMeshValid(p_mpmesh);
704+
705+
auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
706+
//Rank information
707+
int self;
708+
MPI_Comm comm = p_MPs->getMPIComm();
709+
MPI_Comm_rank(comm, &self);
710+
//Asserts
711+
PMT_ALWAYS_ASSERT(nComps == 1);
712+
PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
713+
//MP Data
714+
auto mpReplacementPressure = p_MPs->getData<polyMPO::MPF_ReplacementPressure>();
715+
auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
716+
//Copy to device
717+
kkViewHostU<const double**> mpReplacementPressure_h(replacementPressureMPHost, nComps, numMPs);
718+
Kokkos::View<double**> mpReplacementPressure_d("mpIcePressureDevice", nComps, numMPs);
719+
Kokkos::deep_copy(mpReplacementPressure_d, mpReplacementPressure_h);
720+
//Set in PS
721+
auto setMPReplacementPressure = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
722+
if(mask){
723+
mpReplacementPressure(mp,0) = mpReplacementPressure_d(0, mpAppID(mp));
724+
}
725+
};
726+
p_MPs->parallel_for(setMPReplacementPressure, "setReplacementPressure");
727+
pumipic::RecordTime("PolyMPO_setReplacementPressure" + std::to_string(self), timer.seconds());
728+
}
729+
647730
void polympo_startMeshFill_f(MPMesh_ptr p_mpmesh){
648731
checkMPMeshValid(p_mpmesh);
649732
((polyMPO::MPMesh*)p_mpmesh)->p_mesh->setMeshEdit(true);
@@ -1216,6 +1299,34 @@ void polyMPO_setTanLatVertexRotatedOverRadius_f(MPMesh_ptr p_mpmesh, const int n
12161299
Kokkos::deep_copy(tanLatVertexRotatedOverRadius, h_tanLatVertexRotatedOverRadius);
12171300
}
12181301

1302+
void polympo_setElasticTimeStep_f(MPMesh_ptr p_mpmesh, const double elasticTimeStep){
1303+
//chech validity
1304+
checkMPMeshValid(p_mpmesh);
1305+
auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
1306+
p_mesh->setElasticTimeStep(elasticTimeStep);
1307+
}
1308+
1309+
void polympo_setDynamicTimeStep_f(MPMesh_ptr p_mpmesh, const double dynamicTimeStep){
1310+
//chech validity
1311+
checkMPMeshValid(p_mpmesh);
1312+
auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
1313+
p_mesh->setDynamicTimeStep(dynamicTimeStep);
1314+
}
1315+
1316+
void polympo_setSolveStressMesh_f(MPMesh_ptr p_mpmesh, const int nCells, int* array){
1317+
//chech validity
1318+
checkMPMeshValid(p_mpmesh);
1319+
auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
1320+
1321+
PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells);
1322+
//copy the host array to the device
1323+
auto solveStress = p_mesh->getMeshField<polyMPO::MeshF_SolveStress>();
1324+
auto h_solveStress = Kokkos::create_mirror_view(solveStress);
1325+
for(int i=0; i<nCells; i++)
1326+
h_solveStress(i) = array[i];
1327+
Kokkos::deep_copy(solveStress, h_solveStress);
1328+
}
1329+
12191330
//Advection Calcualtions
12201331
void polympo_push_f(MPMesh_ptr p_mpmesh){
12211332
checkMPMeshValid(p_mpmesh);

‎src/pmpo_c.h‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -48,8 +48,11 @@ void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs,
4848
void polympo_getMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpVelHost);
4949
void polympo_setMPStrainRate_f(MPMesh_ptr p_mpmesh);
5050
void polympo_getMPStrainRate_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpStrainRateHost);
51-
void polympo_setMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStressIn);
51+
void polympo_setMPStress_f(MPMesh_ptr p_mpmesh);
5252
void polympo_getMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpStressHost);
53+
void polympo_setAreaMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* areaMPHost);
54+
void polympo_setIcePressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* icePressureMPHost);
55+
void polympo_setReplacementPressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* replacementPressureMPHost);
5356

5457
//Mesh info
5558
void polympo_startMeshFill_f(MPMesh_ptr p_mpmesh);
@@ -95,6 +98,9 @@ void polympo_setMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices,
9598
void polympo_getMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices, double* areaTriangle);
9699
void polympo_setGnomonicProjection_f(MPMesh_ptr p_mpmesh);
97100
void polyMPO_setTanLatVertexRotatedOverRadius_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array);
101+
void polympo_setElasticTimeStep_f(MPMesh_ptr p_mpmesh, const double elasticTimeStep);
102+
void polympo_setDynamicTimeStep_f(MPMesh_ptr p_mpmesh, const double dynamicTimeStep);
103+
void polympo_setSolveStressMesh_f(MPMesh_ptr p_mpmesh, const int nCells, int* array);
98104

99105
// Advection calculations
100106
void polympo_push_f(MPMesh_ptr p_mpmesh);

‎src/pmpo_const_relation.hpp‎

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,46 @@
1+
#ifndef POLYMPO_CONSTITUTIVE_RELATION
2+
#define POLYMPO_CONSTITUTIVE_RELATION
3+
4+
#include <stdlib.h>
5+
#include <iostream>
6+
7+
namespace polyMPO{
8+
9+
#define PUNY 1.0e-11
10+
static constexpr double eccentricity = 2.0;
11+
static constexpr double eccentricitySquared = eccentricity * eccentricity;
12+
static constexpr double dampingTimescaleParameter = 0.36;
13+
14+
KOKKOS_INLINE_FUNCTION
15+
void constitutive_evp(const Vec3d& strain, Vec3d& stress, const double& icePressure, double& replacementPressure,
16+
const double& areaMP, const double& dtElastic, const double& dampingTimescale){
17+
printf("In the Constitutive model\n");
18+
auto strainDivergence = strain[0] + strain[1];
19+
auto strainTension = strain[0] - strain[1];
20+
auto strainShearing = 2*strain[2];
21+
22+
auto stress1 = stress[0] + stress[1];
23+
auto stress2 = stress[0] - stress[1];
24+
25+
auto Delta = sqrt(strainDivergence*strainDivergence + (strainTension*strainTension + strainShearing*strainShearing)/eccentricitySquared);
26+
27+
auto pressureCoefficient = icePressure / Kokkos::max(Delta, PUNY);
28+
replacementPressure = pressureCoefficient * Delta;
29+
pressureCoefficient = (pressureCoefficient * dtElastic) / (2.0 * dampingTimescale);
30+
31+
auto denominator = 1.0 + (0.5 * dtElastic) / dampingTimescale;
32+
33+
stress1 = (stress1 + pressureCoefficient * (strainDivergence - Delta)) / denominator;
34+
stress2 = (stress2 + (pressureCoefficient / eccentricitySquared) * strainTension ) / denominator;
35+
stress[2] = (stress[2] + (pressureCoefficient / eccentricitySquared) * strainShearing * 0.5) / denominator;
36+
37+
stress[0] = 0.5 * (stress1 + stress2);
38+
stress[1] = 0.5 * (stress1 - stress2);
39+
40+
printf("Stress %.15e %.15e %.15e \n", stress[0], stress[1], stress[2]);
41+
}
42+
43+
}
44+
#endif
45+
46+

‎src/pmpo_fortran.f90‎

Lines changed: 53 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -363,15 +363,43 @@ subroutine polympo_getMPVel(mpMesh, nComps, numMPs, array) &
363363
type(c_ptr), value :: array
364364
end subroutine
365365

366+
!MP Strain
366367
subroutine polympo_setMPStrainRate(mpMesh) &
367368
bind(C, NAME='polympo_setMPStrainRate_f')
368369
use :: iso_c_binding
369370
type(c_ptr), value :: mpMesh
370371
end subroutine
372+
373+
!MP Stress
374+
subroutine polympo_setMPStress(mpMesh) &
375+
bind(C, NAME='polympo_setMPStress_f')
376+
use :: iso_c_binding
377+
type(c_ptr), value :: mpMesh
378+
end subroutine
379+
380+
subroutine polympo_setAreaMP(mpMesh, nComps, numMPs, array) &
381+
bind(C, NAME='polympo_setAreaMP_f')
382+
use :: iso_c_binding
383+
type(c_ptr), value :: mpMesh
384+
integer(c_int), value :: nComps, numMPs
385+
type(c_ptr), value :: array
386+
end subroutine
371387

372-
!getMPStrainRate
373-
!setMPStress
374-
!getMPStress
388+
subroutine polympo_setIcePressureMP(mpMesh, nComps, numMPs, array) &
389+
bind(C, NAME='polympo_setIcePressureMP_f')
390+
use :: iso_c_binding
391+
type(c_ptr), value :: mpMesh
392+
integer(c_int), value :: nComps, numMPs
393+
type(c_ptr), value :: array
394+
end subroutine
395+
396+
subroutine polympo_setReplacementPressureMP(mpMesh, nComps, numMPs, array) &
397+
bind(C, NAME='polympo_setReplacementPressureMP_f')
398+
use :: iso_c_binding
399+
type(c_ptr), value :: mpMesh
400+
integer(c_int), value :: nComps, numMPs
401+
type(c_ptr), value :: array
402+
end subroutine
375403

376404
!---------------------------------------------------------------------------
377405
!> @brief Enable the setting of mesh topology (number of entities and entity adjacencies).
@@ -900,6 +928,28 @@ subroutine polyMPO_setTanLatVertexRotatedOverRadius(mpMesh, nVertices, array) &
900928
type(c_ptr), value :: array
901929
end subroutine
902930

931+
subroutine polympo_setElasticTimeStep(mpMesh, elasticTimeStep) &
932+
bind(C, NAME='polympo_setElasticTimeStep_f')
933+
use :: iso_c_binding
934+
type(c_ptr), value :: mpMesh
935+
real(c_double), value :: elasticTimeStep
936+
end subroutine
937+
938+
subroutine polympo_setDynamicTimeStep(mpMesh, dynamicTimeStep) &
939+
bind(C, NAME='polympo_setDynamicTimeStep_f')
940+
use :: iso_c_binding
941+
type(c_ptr), value :: mpMesh
942+
real(c_double), value :: dynamicTimeStep
943+
end subroutine
944+
945+
subroutine polympo_setSolveStressMesh(mpMesh, nCells, array) &
946+
bind(C, NAME='polympo_setSolveStressMesh_f')
947+
use :: iso_c_binding
948+
type(c_ptr), value :: mpMesh
949+
integer(c_int), value :: nCells
950+
type(c_ptr), value :: array
951+
end subroutine
952+
903953
!---------------------------------------------------------------------------
904954
!> @brief calculate the MPs from given mesh vertices rotational latitude
905955
!> longitude, update the MP slices

0 commit comments

Comments
 (0)