First commit to port

Author: dhyan1272

src/pmpo_c.cpp

@@ -1448,6 +1448,20 @@ void polympo_setSolveVelocityMesh_f(MPMesh_ptr p_mpmesh, const int nVertices, in
   Kokkos::deep_copy(solveVelocity, h_solveVelocity);
 }
 
+void polympo_setIceAreaVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array){
+  //chech validity
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+
+  PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices);
+  //copy the host array to the device
+  auto iceArea = p_mesh->getMeshField<polyMPO::MeshF_VtxMass>();
+  auto h_iceArea = Kokkos::create_mirror_view(iceArea);
+  for(int i=0; i<nVertices; i++)
+    h_iceArea(i, 0) = array[i];
+  Kokkos::deep_copy(iceArea, h_iceArea);
+}
+
 void polympo_calculateStressDivergence_f(MPMesh_ptr p_mpmesh){
   //chech validity
   checkMPMeshValid(p_mpmesh);
@@ -1570,6 +1584,11 @@ void polympo_set_oceanStressCoefficient_f(MPMesh_ptr p_mpmesh, const int nVertic
   Kokkos::deep_copy(oceanStressCoeff, h_oceanStressCoeff);
 }
 
+void polympo_calculate_oceanStressCoefficient_f(MPMesh_ptr p_mpmesh){
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+  p_mesh->calcOceanStressCoeff();
+}
+
 void polympo_velocity_grid_solve_f(MPMesh_ptr p_mpmesh){
   //Temporary grid Solve after calculateDivergence
   auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;

src/pmpo_c.h

@@ -106,6 +106,7 @@ 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_setSolveVelocityMesh_f(MPMesh_ptr p_mpmesh, const int nVertices, int* array);
+void polympo_setIceAreaVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array);
 void polympo_calculateStressDivergence_f(MPMesh_ptr p_mpmesh);
 void polympo_getStressDivergence_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uArray, double* vArray);
 void polympo_setTotalMassVtx_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array);
@@ -114,6 +115,7 @@ void polympo_set_surfaceTiltForce_f(MPMesh_ptr p_mpmesh, const int nVertices, do
 void polympo_set_totalMassVertexfVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array);
 void polympo_set_oceanStress_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uArray, double* vArray);
 void polympo_set_oceanStressCoefficient_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array);
+void polympo_calculate_oceanStressCoefficient_f(MPMesh_ptr p_mpmesh);
 void polympo_velocity_grid_solve_f(MPMesh_ptr p_mpmesh);
 void polympo_set_boundary_normal_vertex_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* uArray, double* vArray);
 void polympo_set_free_slip_bc_f(MPMesh_ptr p_mpmesh);

src/pmpo_fortran.f90

@@ -1011,6 +1011,14 @@ subroutine polympo_setSolveVelocityMesh(mpMesh, nVertices, array) &
     type(c_ptr), value :: array
   end subroutine
 
+ subroutine polympo_setIceAreaVertex(mpMesh, nVertices, array) &
+             bind(C, NAME='polympo_setIceAreaVertex_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nVertices
+    type(c_ptr), value :: array
+  end subroutine
+
   subroutine polympo_calculateStressDivergence(mpMesh) &
              bind(C, NAME='polympo_calculateStressDivergence_f')
     use :: iso_c_binding
@@ -1073,6 +1081,12 @@ subroutine polympo_set_oceanStressCoefficient(mpMesh, nVertices, array) &
     type(c_ptr), value :: array
   end subroutine
 
+  subroutine polympo_calculate_oceanStressCoefficient(mpMesh) &
+             bind(C, NAME='polympo_calculate_oceanStressCoefficient_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+  end subroutine
+
   subroutine  polympo_velocity_grid_solve(mpMesh) &
               bind(C, NAME='polympo_velocity_grid_solve_f')
     use :: iso_c_binding

src/pmpo_mesh.cpp

@@ -157,6 +157,18 @@ namespace polyMPO{
     });
   }
 
+  void Mesh::calcOceanStressCoeff(){
+    int numVerticesOwned = getNumVerticesOwned();
+    auto iceAreaVtx = getMeshField<polyMPO::MeshF_VtxMass>();
+    auto oceanStressCoeff = getMeshField<MeshF_OceanStressCoeff>();
+    auto velocity = getMeshField<MeshF_Vel>();
+    auto solve_velocity = getMeshField<MeshF_SolveVelocity>();
+
+    Kokkos::parallel_for("calcOceanStressCoeff", numVerticesOwned, KOKKOS_LAMBDA(const int vtx){
+      if(solve_velocity(vtx) == 0) return;
+      oceanStressCoeff(vtx, 0) = 0.00536 * 1026.0 * iceAreaVtx(vtx, 0) * sqrt(velocity(vtx, 0)*velocity(vtx, 0) + velocity(vtx, 1)*velocity(vtx, 1));
+    });
+  }
 
   void Mesh::gridSolveGPU(){
     //Mesh Fields

src/pmpo_mesh.hpp

@@ -271,6 +271,8 @@ class Mesh {
     double getDynamicTimeStep(){
       return dynamicTimeStep_;
     }
+
+    void calcOceanStressCoeff(); 
     void gridSolveGPU();
     void aggregateDeluDyn();
     void applyFreeSlipBC();