Passing oceanDragCoeff and using same vel corrections as rebase 29 June

Author: dhyan1272

src/pmpo_MPMesh_assembly.hpp

@@ -342,16 +342,24 @@ void MPMesh::assemblyVtx1(){
   p_mesh->fillMeshField<meshFieldIndex>(numVtx, numEntries, 0.0);
   auto meshField = p_mesh->getMeshField<meshFieldIndex>();
 
+  auto vtxRotLon = p_mesh->getMeshField<MeshF_VtxRotLon>();
+
   //Material Points
   auto mpData = p_MPs->getData<mpfIndex>();
   auto weight = p_MPs->getData<MPF_Basis_Vals>();
   auto mpPositions = p_MPs->getData<MPF_Cur_Pos_XYZ>();
-
+  auto curPosRotLatLon = p_MPs->getData<MPF_Cur_Pos_Rot_Lat_Lon>();
+  auto MPsAppID = p_MPs->getData<MPF_MP_APP_ID>();
   //Earth Radius
   double radius = 1.0;
   if(p_mesh->getGeomType() == geom_spherical_surf)
     radius=p_mesh->getSphereRadius();
 
+  bool use_correction_term = false;
+  if constexpr (meshFieldIndex == MeshF_Vel) {
+    use_correction_term = true;
+  }
+
   //Reconstruct
   auto reconstruct = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
     if(mask) { //if material point is 'active'/'enabled'
@@ -367,9 +375,20 @@ void MPMesh::assemblyVtx1(){
                                                        VtxCoeffs_new(vID,0, 2)*CoordDiffs[2] +
                                                        VtxCoeffs_new(vID,0, 3)*CoordDiffs[3]);
 
-        for (int k=0; k<numEntries; k++){
-          auto val = factor*mpData(mp,k);
-          Kokkos::atomic_add(&meshField(vID,k), val);
+        if (use_correction_term){
+          double transport[2] = {0.0};
+          seaice_mpm_coord_parallel_transport(curPosRotLatLon(mp, 1), vtxRotLon(vID), curPosRotLatLon(mp, 0), transport);
+          const double u =  transport[0] * mpData(mp,0) + transport[1] * mpData(mp,1);
+          const double v = -transport[1] * mpData(mp,0) + transport[0] * mpData(mp,1);
+
+          Kokkos::atomic_add(&meshField(vID,0), factor * u);
+          Kokkos::atomic_add(&meshField(vID,1), factor * v);
+        }
+        else{
+          for (int k=0; k<numEntries; k++){
+            auto val = factor*mpData(mp,k);
+            Kokkos::atomic_add(&meshField(vID,k), val);
+          }
         }
       }
     }

src/pmpo_c.cpp

@@ -1207,6 +1207,25 @@ void polympo_getMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, double
   }
 }
 
+void polympo_setMeshVtxRotLon_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* longitude){
+  Kokkos::Timer timer;
+  //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 host array to the device
+  auto coordsArray = p_mesh->getMeshField<polyMPO::MeshF_VtxRotLon>();
+  auto h_coordsArray = Kokkos::create_mirror_view(coordsArray);
+  for(int i=0; i<nVertices; i++){
+    h_coordsArray(i) = longitude[i];
+  }
+  Kokkos::deep_copy(coordsArray, h_coordsArray);
+  pumipic::RecordTime("PolyMPO_setMeshVtxRotLon", timer.seconds());
+}
+
 void polympo_setMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* uVelIn, const double* vVelIn){
   //check mpMesh is valid
   checkMPMeshValid(p_mpmesh);
@@ -1674,9 +1693,9 @@ 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){
+void polympo_calculate_oceanStressCoefficient_f(MPMesh_ptr p_mpmesh, const double configIceOceanDragCoeff){
   auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
-  p_mesh->calcOceanStressCoeff();
+  p_mesh->calcOceanStressCoeff(configIceOceanDragCoeff);
 }
 
 void polympo_velocity_grid_solve_f(MPMesh_ptr p_mpmesh){

src/pmpo_c.h

@@ -89,6 +89,7 @@ int polympo_getMeshFElmType_f();
 void polympo_setMeshVtxCoords_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* xArray, const double* yArray, const double* zArray);
 void polympo_getMeshVtxCoords_f(MPMesh_ptr p_mpmesh, const int nVertices, double* xArray, double* yArray, double* zArray);
 void polympo_setMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* latitude);
+void polympo_setMeshVtxRotLon_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* longitude);
 void polympo_getMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, double* latitude);
 void polympo_setMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* uVelocity, const double* vVelocity);
 void polympo_getMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uVelocity, double* vVelocity);
@@ -119,7 +120,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_calculate_oceanStressCoefficient_f(MPMesh_ptr p_mpmesh, const double configIceOceanDragCoeff);
 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

@@ -786,6 +786,14 @@ subroutine polympo_getMeshVtxRotLat(mpMesh, nVertices, latitude) &
     type(c_ptr), value :: latitude
   end subroutine
 
+  subroutine polympo_setMeshVtxRotLon(mpMesh, nVertices, longitude) &
+             bind(C, NAME='polympo_setMeshVtxRotLon_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nVertices
+    type(c_ptr), intent(in), value :: longitude
+  end subroutine
+
   !---------------------------------------------------------------------------
   !> @brief set the vertices velocity from a host array
   !> @param mpmesh(in/out) MPMesh object
@@ -1111,10 +1119,11 @@ subroutine polympo_set_oceanStressCoefficient(mpMesh, nVertices, array) &
     type(c_ptr), value :: array
   end subroutine
 
-  subroutine polympo_calculate_oceanStressCoefficient(mpMesh) &
+  subroutine polympo_calculate_oceanStressCoefficient(mpMesh, configIceOceanDragCoeff) &
              bind(C, NAME='polympo_calculate_oceanStressCoefficient_f')
     use :: iso_c_binding
     type(c_ptr), value :: mpMesh
+    real(c_double), value::configIceOceanDragCoeff
   end subroutine
 
   subroutine  polympo_velocity_grid_solve(mpMesh) &

src/pmpo_mesh.cpp

@@ -20,6 +20,10 @@ namespace polyMPO{
     PMT_ALWAYS_ASSERT(vtxRotLatMapEntry.first == MeshFType_VtxBased);
     vtxRotLat_ = MeshFView<MeshF_VtxRotLat>(vtxRotLatMapEntry.second,numVtxs_);
 
+    auto vtxRotLonMapEntry = meshFields2TypeAndString.at(MeshF_VtxRotLon);
+    PMT_ALWAYS_ASSERT(vtxRotLonMapEntry.first == MeshFType_VtxBased);
+    vtxRotLon_ = MeshFView<MeshF_VtxRotLon>(vtxRotLonMapEntry.second,numVtxs_);
+
     auto vtxVelMapEntry = meshFields2TypeAndString.at(MeshF_Vel);
     PMT_ALWAYS_ASSERT(vtxVelMapEntry.first == MeshFType_VtxBased);
     vtxVel_ = MeshFView<MeshF_Vel>(vtxVelMapEntry.second,numVtxs_);
@@ -160,18 +164,19 @@ namespace polyMPO{
     });
   }
 
-  void Mesh::calcOceanStressCoeff(){
+  void Mesh::calcOceanStressCoeff(const double configIceOceanDragCoeff){
     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>();
     auto oceanVelocity = getMeshField<MeshF_OceanVelocity>();
+    auto seaiceDensitySeaWater_ = polyMPO::seaiceDensitySeaWater;
 
     Kokkos::parallel_for("calcOceanStressCoeff", numVerticesOwned, KOKKOS_LAMBDA(const int vtx){
       if(solve_velocity(vtx) == 0) return;
       auto relVelSq = pow(oceanVelocity(vtx, 0) - velocity(vtx, 0), 2) +  pow(oceanVelocity(vtx, 1) - velocity(vtx, 1), 2);
-      oceanStressCoeff(vtx, 0) = 0.00536 * 1026.0 * iceAreaVtx(vtx, 0) * sqrt(relVelSq);
+      oceanStressCoeff(vtx, 0) = configIceOceanDragCoeff * seaiceDensitySeaWater_ * iceAreaVtx(vtx, 0) * sqrt(relVelSq);
     });
   }
 

src/pmpo_mesh.hpp

@@ -19,6 +19,7 @@ enum MeshFieldIndex{
     MeshF_Unsupported,
     MeshF_VtxCoords,
     MeshF_VtxRotLat,
+    MeshF_VtxRotLon,
     MeshF_ElmCenterXYZ,
     MeshF_DualTriangleArea,
     MeshF_Vel,
@@ -55,6 +56,7 @@ enum MeshFieldType{
 template <MeshFieldIndex> struct meshFieldToType;
 template <> struct meshFieldToType < MeshF_VtxCoords         > { using type = Kokkos::View<vec3d_t*>; };
 template <> struct meshFieldToType < MeshF_VtxRotLat         > { using type = DoubleView; };
+template <> struct meshFieldToType < MeshF_VtxRotLon         > { using type = DoubleView; };
 template <> struct meshFieldToType < MeshF_ElmCenterXYZ      > { using type = Kokkos::View<vec3d_t*>; };
 template <> struct meshFieldToType < MeshF_DualTriangleArea  > { using type = Kokkos::View<doubleSclr_t*>; };
 template <> struct meshFieldToType < MeshF_Vel               > { using type = Kokkos::View<vec2d_t*>; };
@@ -88,6 +90,7 @@ const std::map<MeshFieldIndex, std::pair<MeshFieldType, std::string>> meshFields
         {MeshF_Unsupported,      {MeshFType_Unsupported,"MeshField_Unsupported"}},
         {MeshF_VtxCoords,        {MeshFType_VtxBased,"MeshField_VerticesCoords"}},
         {MeshF_VtxRotLat,        {MeshFType_VtxBased,"MeshField_VerticesLatitude"}},
+        {MeshF_VtxRotLon,        {MeshFType_VtxBased,"MeshField_VerticesLongitude"}},
         {MeshF_ElmCenterXYZ,     {MeshFType_ElmBased,"MeshField_ElementCenterXYZ"}},
         {MeshF_DualTriangleArea, {MeshFType_VtxBased,"MeshField_DualTriangleArea"}},
         {MeshF_Vel,              {MeshFType_VtxBased,"MeshField_Velocity"}},
@@ -144,6 +147,7 @@ class Mesh {
     //start of meshFields
     MeshFView<MeshF_VtxCoords> vtxCoords_;
     MeshFView<MeshF_VtxRotLat> vtxRotLat_;
+    MeshFView<MeshF_VtxRotLon> vtxRotLon_;
     MeshFView<MeshF_ElmCenterXYZ> elmCenterXYZ_;
     MeshFView<MeshF_DualTriangleArea> dualTriangleArea_;
 
@@ -159,7 +163,7 @@ class Mesh {
     //GnomonicProjection
     MeshFView<MeshF_VtxGnomProj> vtxGnomProj_;
     MeshFView<MeshF_ElmCenterGnomProj> elmCenterGnomProj_;
-    
+
     MeshFView<MeshF_TanLatVertexRotatedOverRadius> tanLatVertexRotatedOverRadius_;
     MeshFView<MeshF_SolveStress> solveStress_;
     MeshFView<MeshF_SolveVelocity> solveVelocity_;
@@ -276,7 +280,7 @@ class Mesh {
       return dynamicTimeStep_;
     }
 
-    void calcOceanStressCoeff(); 
+    void calcOceanStressCoeff(const double configIceOceanDragCoeff);
     void gridSolveGPU();
     void aggregateDeluDyn();
     void applyFreeSlipBC();
@@ -298,6 +302,9 @@ auto Mesh::getMeshField(){
     else if constexpr (index==MeshF_VtxRotLat){
         return vtxRotLat_;
     }
+    else if constexpr (index==MeshF_VtxRotLon){
+        return vtxRotLon_;
+    }
     else if constexpr (index==MeshF_ElmCenterXYZ){
         return elmCenterXYZ_;
     }

src/pmpo_utils.hpp

@@ -45,6 +45,9 @@ using DoubleView = Kokkos::View<double*>;
 using IntView = Kokkos::View<int*>;
 using BoolView = Kokkos::View<bool*>;
 
+//CONSTANTS
+inline constexpr double seaiceDensitySeaWater = 1026.0;
+
 class Vec2d {
   private:
     vec2d_t coords_;
@@ -554,6 +557,24 @@ void lat_lon_from_xyz(double& lat, double& lon, Vec3d& xyz, double r){
   lat = Kokkos::asin(xyz[2]/r);
 }
 
+KOKKOS_INLINE_FUNCTION
+double seaice_mpm_wrap_longitude(const double longitude){
+  const auto PI = 3.141592653589;
+  const auto TAU = 2*PI;
+  const double wrapped = longitude - Kokkos::floor((longitude + PI) / TAU) * TAU;
+  return wrapped;
+}
+
+KOKKOS_INLINE_FUNCTION
+void seaice_mpm_coord_parallel_transport(const double long_start, const double long_end,
+                                         const double latitude, double transport[2]){
+  double dLon = long_end - long_start;
+  dLon = seaice_mpm_wrap_longitude(dLon);
+  const double psi = dLon * Kokkos::sin(latitude);
+  transport[0] = Kokkos::cos(psi);
+  transport[1] = Kokkos::sin(psi);
+}
+
 }//namespace polyMPO end
 
 #endif

src/pmpo_wachspressBasis.hpp

@@ -15,25 +15,48 @@ void sphericalInterpolation(MPMesh& mpMesh){
   auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
   int numVtxs = p_mesh->getNumVertices();
   auto elm2VtxConn = p_mesh->getElm2VtxConn();
+  auto vtxRotLon = p_mesh->getMeshField<MeshF_VtxRotLon>();
 
   auto p_MPs = mpMesh.p_MPs;
   auto MPsPosition = p_MPs->getPositions();
   auto MPsBasis = p_MPs->getData<MPF_Basis_Vals>();
+  auto curPosRotLatLon = p_MPs->getData<MPF_Cur_Pos_Rot_Lat_Lon>();
+  auto MPsAppID = p_MPs->getData<MPF_MP_APP_ID>();
 
   constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice<meshFieldIndex>;
   auto mpField = p_MPs->getData<mpfIndex>();
 
   const int numEntries = mpSliceToNumEntries<mpfIndex>();
   auto meshField = p_mesh->getMeshField<meshFieldIndex>(); 
 
+  bool use_correction_term = false;
+  if constexpr (meshFieldIndex == MeshF_OnSurfVeloIncr) {
+    use_correction_term = true;
+  }
+
   auto interpolation = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
     if(mask) { //if material point is 'active'/'enabled'
       int numVtx = elm2VtxConn(elm,0);
-      for(int entry=0; entry<numEntries; entry++){
-        double mpValue = 0.0;
-        for(int i=1; i<= numVtx; i++)
-          mpValue += meshField(elm2VtxConn(elm,i)-1,entry)*MPsBasis(mp,i-1);
-        mpField(mp,entry) = mpValue;
+
+      if(use_correction_term){
+        double transport[2] = {0.0};
+        double mpField_l[2] = {0.0};
+        for(int i=1; i<=numVtx; i++){
+          int vID = elm2VtxConn(elm,i)-1;
+          seaice_mpm_coord_parallel_transport(vtxRotLon(vID), curPosRotLatLon(mp, 1), curPosRotLatLon(mp, 0), transport);
+          mpField_l[0] += ( transport[0] * meshField(vID, 0) + transport[1] * meshField(vID, 1)) * MPsBasis(mp,i-1);
+          mpField_l[1] += (-transport[1] * meshField(vID, 0) + transport[0] * meshField(vID, 1)) * MPsBasis(mp,i-1);
+        }
+        mpField(mp, 0) = mpField_l[0];
+        mpField(mp, 1) = mpField_l[1];
+      }
+      else{
+        for(int entry=0; entry<numEntries; entry++){
+          double mpValue = 0.0;
+          for(int i=1; i<= numVtx; i++)
+            mpValue += meshField(elm2VtxConn(elm,i)-1,entry)*MPsBasis(mp,i-1);
+          mpField(mp,entry) = mpValue;
+        }
       }
     }
   };