StressDivergence within Elastic Time Loop compared

Author: dhyan1272

src/pmpo_MPMesh.cpp

@@ -52,9 +52,6 @@ void MPMesh::calculateStrain(){
         v22 = v22 + MPsBasisGrads(mp, i*2 + 1) * velField(iVertex, 1);
         uTanOverR = uTanOverR + MPsBasis(mp, i) * tanLatVertexRotatedOverRadius(iVertex, 0) * velField(iVertex, 0);
         vTanOverR = vTanOverR + MPsBasis(mp, i) * tanLatVertexRotatedOverRadius(iVertex, 0) * velField(iVertex, 1);
-        MPsStrainRate(mp, 0) =  v11 - vTanOverR;
-        MPsStrainRate(mp, 1) =  v22;
-        MPsStrainRate(mp, 2) =  0.5*(v12 + v21 + uTanOverR);
         //Debugging
         if(MPsAppID(mp)==0){
           printf("Strain Calc: iVertex %d vel field %.15e %.15e \n", iVertex, velField(iVertex, 0), velField(iVertex, 1));
@@ -86,7 +83,7 @@ void MPMesh::calculateStress(){
     if(mask){
 
       Vec3d strain_rate (MPsStrainRate(mp, 0), MPsStrainRate(mp, 1), MPsStrainRate(mp, 2));
-      Vec3d stress(MPsStress(mp, 0), MPsStress(mp, 1), MPsStress(mp, 1));
+      Vec3d stress(MPsStress(mp, 0), MPsStress(mp, 1), MPsStress(mp, 2));
       constitutive_evp(strain_rate, stress, MPsIcePressure(mp, 0), MPsRepPressure(mp, 0), MPsArea(mp, 0), elasticTimeStep, dampingTimescale);
       for (int m=0 ; m<3; m++)
         MPsStress(mp, m) = stress[m];
@@ -172,26 +169,24 @@ void MPMesh::calculateStressDivergence(){
   };
   p_MPs->parallel_for(stress_div, "assembly");
 
-  //TODO
-  //COMMUNICATE THE VERTEX FIELDS
+  //TODO COMMUNICATE THE VERTEX FIELDS
 
-  //TODO put as mesh field
-  Kokkos::View<vec2d_t*> stressDivergence("stressDivergence", p_mesh->getNumVertices());
+ 
+  auto stressDivergence = p_mesh->getMeshField<MeshF_StressDivergence>();
 
-  Kokkos::parallel_for("calculate_divergence", numVtx, KOKKOS_LAMBDA(const int vtx){
-   
+  Kokkos::parallel_for("calculate_divergence", numVtx, KOKKOS_LAMBDA(const int vtx){  
     stressDivergence(vtx, 0) = stress_divU(vtx);
     stressDivergence(vtx, 1) = stress_divV(vtx); 
     //Debugging
     if (vtx >= 10 && vtx <= 11) {
-      printf("Vtx %d Divergence %.15e %.15e %.15e %.15e %.15e %.15e \n", vtx, nearAnEdge_l(vtx), vtxMatrixMass_l(vtx),
-                                                                         stress_divU(vtx), stress_divV(vtx),
-                                                                         divU_edge(vtx), divV_edge(vtx));
+      //printf("Vtx %d Divergence %.15e %.15e %.15e %.15e %.15e %.15e \n", vtx, nearAnEdge_l(vtx), vtxMatrixMass_l(vtx),
+      //                                                                   stress_divU(vtx), stress_divV(vtx),
+      //                                                                   divU_edge(vtx), divV_edge(vtx));
+      printf("Vtx %d Divergence %.15e %.15e \n", vtx, stressDivergence(vtx, 0), stressDivergence(vtx, 1));
     }
   });
-   
-  //TODO
-  //COMMUNICATE THE VERTEX FIELDS
+  
+  //TODO COMMUNICATE THE VERTEX FIELDS
 
 }
 

src/pmpo_MPMesh_assembly.hpp

@@ -382,15 +382,6 @@ void MPMesh::assemblyVtx1(){
   if(numProcsTot>1) 
     communicate_and_take_halo_contributions(meshField, numVertices, numEntries, 0, 0);
   pumipic::RecordTime("Communicate Field Values" + std::to_string(self), timer.seconds());
-
-  Kokkos::parallel_for("printSymmetricBlock", numVertices, KOKKOS_LAMBDA(const int vtx){
-    if (vtx >= 10 && vtx <= 10) {
-      printf("Field in %d:  ", vtx);
-      for (int k=0; k<numEntries; k++)
-        printf("  %.15e ", meshField(vtx, k));
-      printf("\n");
-    }
-  });
 }
 
 //Start Communication routine

src/pmpo_c.cpp

@@ -1349,6 +1349,24 @@ void polympo_calculateStressDivergence_f(MPMesh_ptr p_mpmesh){
   ((polyMPO::MPMesh*)p_mpmesh) -> calculateStressDivergence();
 }
 
+void polympo_getStressDivergence_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uArray, double* vArray){
+  //chech validity
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+
+  //check the size
+  PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); 
+
+  //copy the device to host 
+  auto stressDivergence = p_mesh->getMeshField<polyMPO::MeshF_StressDivergence>();
+  auto h_stressDivergence = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), stressDivergence);
+  for(int i=0; i<nVertices; i++){
+    uArray[i] = h_stressDivergence(i,0);
+    vArray[i] = h_stressDivergence(i,1);
+  }
+}
+
+
 //Advection Calcualtions
 void polympo_push_f(MPMesh_ptr p_mpmesh){
   checkMPMeshValid(p_mpmesh);

src/pmpo_c.h

@@ -103,6 +103,8 @@ void polympo_setElasticTimeStep_f(MPMesh_ptr p_mpmesh, const double elasticTimeS
 void polympo_setDynamicTimeStep_f(MPMesh_ptr p_mpmesh, const double dynamicTimeStep);
 void polympo_setSolveStressMesh_f(MPMesh_ptr p_mpmesh, const int nCells, int* array);
 void polympo_calculateStressDivergence_f(MPMesh_ptr p_mpmesh);
+void polympo_getStressDivergence_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uArray, double* vArray);
+
 
 // Advection calculations
 void polympo_push_f(MPMesh_ptr p_mpmesh);

src/pmpo_fortran.f90

@@ -966,6 +966,15 @@ subroutine polympo_calculateStressDivergence(mpMesh) &
     type(c_ptr), value :: mpMesh
   end subroutine
 
+  subroutine polympo_getStressDivergence(mpMesh, nVertices, uArray, vArray) &
+             bind(C, NAME='polympo_getStressDivergence_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nVertices
+    type(c_ptr), value :: uArray, vArray
+  end subroutine
+
+
   !---------------------------------------------------------------------------
   !> @brief calculate the MPs from given mesh vertices rotational latitude
   !>        longitude, update the MP slices

src/pmpo_mesh.cpp

@@ -51,6 +51,10 @@ namespace polyMPO{
     auto interiorVertexEntry = meshFields2TypeAndString.at(MeshF_InteriorVertex);
     PMT_ALWAYS_ASSERT(interiorVertexEntry.first == MeshFType_VtxBased);
     interiorVertex_ = MeshFView<MeshF_InteriorVertex>(interiorVertexEntry.second, numVtxs_);
+
+    auto stressDivergenceEntry = meshFields2TypeAndString.at(MeshF_StressDivergence);
+    PMT_ALWAYS_ASSERT(stressDivergenceEntry.first == MeshFType_VtxBased);
+    stressDivergence_ = MeshFView<MeshF_StressDivergence>(stressDivergenceEntry.second, numVtxs_);
   }
 
   void Mesh::setMeshElmBasedFieldSize(){

src/pmpo_mesh.hpp

@@ -31,7 +31,8 @@ enum MeshFieldIndex{
     MeshF_ElmCenterGnomProj,
     MeshF_TanLatVertexRotatedOverRadius,
     MeshF_SolveStress,
-    MeshF_InteriorVertex
+    MeshF_InteriorVertex,
+    MeshF_StressDivergence
 };
 enum MeshFieldType{
     MeshFType_Invalid = -2,
@@ -56,6 +57,7 @@ template <> struct meshFieldToType < MeshF_ElmCenterGnomProj > { using type = Ko
 template <> struct meshFieldToType < MeshF_TanLatVertexRotatedOverRadius > { using type = Kokkos::View<doubleSclr_t*>; };
 template <> struct meshFieldToType < MeshF_SolveStress       > { using type = IntView; };
 template <> struct meshFieldToType < MeshF_InteriorVertex    > { using type = IntView; };
+template <> struct meshFieldToType < MeshF_StressDivergence  > { using type = Kokkos::View<vec2d_t*>; };
 
 template <MeshFieldIndex index>
 using MeshFView = typename meshFieldToType<index>::type;
@@ -77,7 +79,8 @@ const std::map<MeshFieldIndex, std::pair<MeshFieldType, std::string>> meshFields
         {MeshF_ElmCenterGnomProj,{MeshFType_ElmBased,"MeshField_ElementCenterGnomonicprojection"}},
         {MeshF_TanLatVertexRotatedOverRadius, {MeshFType_VtxBased,"MeshField_TanLatVertexRotatedOverRadius"}},
         {MeshF_SolveStress,      {MeshFType_ElmBased,"MeshField_SolveStress"}},
-        {MeshF_InteriorVertex,   {MeshFType_VtxBased,"MeshField_InteriorVertex"}}
+        {MeshF_InteriorVertex,   {MeshFType_VtxBased,"MeshField_InteriorVertex"}},
+        {MeshF_StressDivergence, {MeshFType_VtxBased,"MeshField_StressDivergence"}},
 };
 
 enum mesh_type {mesh_unrecognized_lower = -1,
@@ -111,6 +114,7 @@ class Mesh {
     MeshFView<MeshF_VtxRotLat> vtxRotLat_;
     MeshFView<MeshF_ElmCenterXYZ> elmCenterXYZ_;
     MeshFView<MeshF_DualTriangleArea> dualTriangleArea_;
+    MeshFView<MeshF_StressDivergence> stressDivergence_;
     MeshFView<MeshF_Vel> vtxVel_;
     MeshFView<MeshF_VtxMass> vtxMass_;
     MeshFView<MeshF_ElmMass> elmMass_;
@@ -281,6 +285,9 @@ auto Mesh::getMeshField(){
     else if constexpr (index==MeshF_InteriorVertex){
         return interiorVertex_;
     }
+    else if constexpr (index==MeshF_StressDivergence){
+        return stressDivergence_;
+    }
     fprintf(stderr,"Mesh Field Index error!\n");
     exit(1);
 }