Debris-friction slip law interface

Adding two fields (bulkFriction, basalDebris) to use with
debris-bed friction slip law.

Author: jeremy-brooks-1

components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.cpp

@@ -69,7 +69,7 @@ std::vector<int> indexToTriangleID,
 std::vector<int> indexToVertexID, vertexToFCell, vertexProcIDs, triangleToFVertex, indexToEdgeID, edgeToFEdge,
     fVertexToTriangleID, fCellToVertex, iceMarginEdgesLIds, dirichletNodesIDs;
 std::vector<double> dissipationHeatOnPrisms, velocityOnVertices, velocityOnCells,
-    elevationData, thicknessData, betaData, bedTopographyData, stiffnessFactorData, effecPressData, muFrictionData, bedRoughnessData, temperatureDataOnPrisms, smbData, thicknessOnCells, bodyForceOnBasalCell;
+    elevationData, thicknessData, betaData, bedTopographyData, stiffnessFactorData, effecPressData, muFrictionData, bedRoughnessData, bulkFrictionData, basalDebrisData, temperatureDataOnPrisms, smbData, thicknessOnCells, bodyForceOnBasalCell;
 std::vector<bool> isVertexBoundary, isBoundaryEdge;
 
 // only needed for creating ASCII mesh
@@ -241,7 +241,7 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower
     double const* thickness_F, double * beta_F,
     double const* smb_F, double const* temperature_F, double const* stiffnessFactor_F,
     double* effecPress_F, double const* muFriction_F,
-    double const* bedRoughnessRC_F,
+    double const* bedRoughnessRC_F, double const* bulkFriction_F, double const* basalDebris_F,
     double* const dirichletVelocityXValue, double* const dirichletVelocitYValue,
     double* u_normal_F, double* bodyForce_F, double* dissipation_heat_F,
     double* xVelocityOnCell, double* yVelocityOnCell, double const* deltat,
@@ -270,7 +270,7 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower
   if (!isDomainEmpty) {
 
     std::vector<std::pair<int, int> > marineBdyExtensionMap;
-    importFields(marineBdyExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, stiffnessFactor_F, effecPress_F, muFriction_F, bedRoughnessRC_F, temperature_F, smb_F,  minThickness);
+    importFields(marineBdyExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, stiffnessFactor_F, effecPress_F, muFriction_F, bedRoughnessRC_F, bulkFriction_F, basalDebris_F, temperature_F, smb_F,  minThickness);
 
     std::vector<double> regulThk(thicknessData);
     for (int index = 0; index < nVertices; index++)
@@ -287,7 +287,7 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower
         Ordering, first_time_step, indexToVertexID, indexToTriangleID, minBeta,
         regulThk, levelsNormalizedThickness, elevationData, thicknessData,
         betaData, bedTopographyData, smbData,
-        stiffnessFactorData, effecPressData, muFrictionData, bedRoughnessData,
+        stiffnessFactorData, effecPressData, muFrictionData, bedRoughnessData, bulkFrictionData, basalDebrisData,
         temperatureDataOnPrisms, bodyForceOnBasalCell, dissipationHeatOnPrisms, velocityOnVertices,
         albany_error, dt);
     *error=albany_error;
@@ -1120,7 +1120,7 @@ double signedTriangleAreaOnSphere(const double* x, const double* y,
 
 void importFields(std::vector<std::pair<int, int> >& marineBdyExtensionMap,  double const* bedTopography_F, double const * lowerSurface_F, double const * thickness_F,
     double const * beta_F, double const* stiffnessFactor_F, double const* effecPress_F, double const* muFriction_F,
-    double const* bedRoughnessRC_F,
+    double const* bedRoughnessRC_F, double const* bulkFriction_F, double const* basalDebris_F,
     double const * temperature_F, double const * smb_F, double eps) {
 
   int vertexLayerShift = (Ordering == 0) ? 1 : nLayers + 1;
@@ -1137,6 +1137,10 @@ void importFields(std::vector<std::pair<int, int> >& marineBdyExtensionMap,  dou
     muFrictionData.assign(nVertices, 1e10);
   if (bedRoughnessRC_F!= 0)
     bedRoughnessData.assign(nVertices, 1e10);
+  if (bulkFriction_F != 0)
+    bulkFrictionData.assign(nVertices, 1e10);
+  if (basalDebris_F != 0)
+    basalDebrisData.assign(nVertices, 1e10);
   if(temperature_F != 0)
     temperatureDataOnPrisms.assign(nLayers * nTriangles, 1e10);
   if (smb_F != 0)
@@ -1161,6 +1165,10 @@ void importFields(std::vector<std::pair<int, int> >& marineBdyExtensionMap,  dou
       muFrictionData[index] = muFriction_F[iCell];
     if (bedRoughnessRC_F != 0)
       bedRoughnessData[index] = bedRoughnessRC_F[iCell] / unit_length;
+    if (bulkFriction_F != 0)
+      bulkFrictionData[index] = bulkFriction_F[iCell];
+    if (basalDebris_F != 0)
+      basalDebrisData[index] = basalDebris_F[iCell];
   }
 
   int lElemColumnShift = (Ordering == 1) ? 1 : nTriangles;
@@ -1650,7 +1658,7 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep
     double const* surfaceAirTemperature_F, double const* basalHeatFlux_F,
     double const* stiffnessFactor_F,
     double const* effecPress_F, double const* muFriction_F,
-    double const* bedRoughnessRC_F,
+    double const* bedRoughnessRC_F, double const* bulkFriction_F, double const* basalDebris_F,
     double const* thickness_F, double const* thicknessUncertainty_F,
     double const* smb_F, double const* smbUncertainty_F,
     double const* bmb_F, double const* bmbUncertainty_F,
@@ -1736,7 +1744,7 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep
 
     std::vector<std::pair<int, int> > marineBdyExtensionMap;  // local map to be created by importFields
     importFields(marineBdyExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F,
-                   stiffnessFactor_F, effecPress_F, muFriction_F, bedRoughnessRC_F,
+                   stiffnessFactor_F, effecPress_F, muFriction_F, bedRoughnessRC_F, bulkFriction_F, basalDebris_F,
                    temperature_F, smb_F, minThickness);
 
     import2DFieldsObservations(marineBdyExtensionMap,
@@ -1803,6 +1811,7 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep
 
     write_ascii_mesh_field(effecPressData, "effective_pressure");
 
+
     // These two fields are more complicated than the others so cannot use the function to write out
 
     std::cout << "Writing temperature.ascii." << std::endl;

components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.hpp

@@ -90,7 +90,7 @@ void velocity_solver_solve_l1l2(double const* lowerSurface_F,
 
 void velocity_solver_solve_fo(double const* bedTopography_F, double const* lowerSurface_F,
     double const* thickness_F, double * beta_F, double const* smb_F, double const* temperature_F, double const* stiffnessFactor_F,
-    double* effecPress_F, double const* muFriction_F, double const* bedRoughnessRC_F,
+    double* effecPress_F, double const* muFriction_F, double const* bedRoughnessRC_F, double const* bulkFriction_F, double const* basalDebris_F,
     double* const dirichletVelocityXValue = 0, double* const dirichletVelocitYValue = 0,
     double* u_normal_F = 0, double* bodyForce_F = 0, double* dissipation_heat_F = 0,
     double* xVelocityOnCell = 0, double* yVelocityOnCell = 0, double const * deltat = 0,
@@ -134,7 +134,7 @@ void write_ascii_mesh(int const* indexToCellID_F,
     double const* surfaceAirTemperature_F, double const* basalHeatFlux_F,
     double const* stiffnessFactor_F,
     double const* effecPress_F, double const* muFriction_F,
-    double const* bedRoughnessRC_F,
+    double const* bedRoughnessRC_F, double const* bulkFriction_F, double const* basalDebris_F, 
     double const* thickness_F, double const* thicknessUncertainty_F,
     double const* smb_F, double const* smbUncertainty_F,
     double const* bmb_F, double const* bmbUncertainty_F,
@@ -165,6 +165,8 @@ extern void velocity_solver_solve_fo__(int nLayers, int nGlobalVertices,
           std::vector<double>& effecPressData,
     const std::vector<double>& muFrictionData,
     const std::vector<double>& bedRoughnessData,
+    const std::vector<double>& bulkFrictionData,
+    const std::vector<double>& basalDebrisData,
     const std::vector<double>& temperatureDataOnPrisms,
     std::vector<double>& bodyForceOnBasalCell,
     std::vector<double>& dissipationHeatOnPrisms,
@@ -209,7 +211,8 @@ void createReducedMPI(int nLocalEntities, MPI_Comm& reduced_comm_id);
 
 void importFields(std::vector<std::pair<int, int> >& marineBdyExtensionMap,
                 double const* bedTopography_F, double const* lowerSurface_F, double const* thickness_F,
-    double const* beta_F = 0, double const* stiffnessFactor_F = 0, double const* effecPress_F = 0, double const* muFriction_F = 0, double const* bedRoughness_F = 0, double const* temperature_F = 0, double const* smb_F = 0, double eps = 0);
+    double const* beta_F = 0, double const* stiffnessFactor_F = 0, double const* effecPress_F = 0, double const* muFriction_F = 0, double const* bedRoughness_F = 0,
+    double const* bulkFriction_F = 0, double const* basalDebris_F = 0, double const* temperature_F = 0, double const* smb_F = 0, double eps = 0);
 
 void import2DFieldsObservations(std::vector<std::pair<int, int> >& marineBdyExtensionMap,
             double const * lowerSurface_F,