Debris-bed friction slip law

Add debris-bed friction slip relation.
Adding two new friction parameters, could set as constants in
albany_input.yaml or read fields from MPAS mesh:
*bulkFrictionCoefficient
*basalDebrisFactor
There are corresponding changes in MPAS interface in MALI-Dev/ES3M
repository.

Author: jeremy-brooks-1

src/landIce/evaluators/LandIce_BasalFrictionCoefficient.hpp

@@ -48,14 +48,18 @@ class BasalFrictionCoefficient : public PHX::EvaluatorWithBaseImpl<Traits>,
   // Coefficients for computing beta (if not given)
   double n;                                             // [adim] exponent of Glen's law
   PHX::MDField<const ScalarT,Dim> muParam;              // [yr^q m^{-q}], friction coefficient of the power Law with exponent q
-  PHX::MDField<const ScalarT,Dim> bedRoughnessParam;    // [km],  Bed bumps avg length divided by bed bumps avg slope (for REGULARIZED_COULOMB only)
-  PHX::MDField<const ScalarT,Dim> powerParam;           // [adim], Exponent (for POWER_LAW and REGULARIZED COULOMB only)
-
+  PHX::MDField<const ScalarT,Dim> bedRoughnessParam;    // [km],  Bed bumps avg length divided by bed bumps avg slope (for REGULARIZED_COULOMB and DEBRIS_FRICTION only)
+  PHX::MDField<const ScalarT,Dim> powerParam;           // [adim], Exponent (for POWER_LAW and REGULARIZED COULOMB and DEBRIS_FRICTION only)
+  PHX::MDField<const ScalarT,Dim> bulkFrictionParam;    // [adim], bulk friction coefficient (for DEBRIS_FRICTION only)
+  PHX::MDField<const ScalarT,Dim> basalDebrisParam;	// [Pa m^-1 s] basal debris factor (for DEBRIS_FRICTION only)
+  
   ScalarT printedMu;
   ScalarT printedBedRoughness;
   ScalarT printedQ;
+  ScalarT printedBulkFriction;
+  ScalarT printedBasalDebris;
 
-  ParamScalarT mu, bedRoughness, power;
+  ParamScalarT mu, bedRoughness, power, bulkFriction, basalDebris;
 
   double beta_val;   // beta value [kPa yr/m] (for CONSTANT only)
   double flowRate_val;   // flow rate value [Pa^{-n} s^{-1}] (for CONSTANT only)
@@ -66,6 +70,8 @@ class BasalFrictionCoefficient : public PHX::EvaluatorWithBaseImpl<Traits>,
   PHX::MDField<const VelocityST>        u_norm;             // [m yr^{-1}]
   PHX::MDField<const ParamScalarT>      bedRoughnessField;  // [km], characteristic length
   PHX::MDField<const ParamScalarT>      muField;            // [yr^q m^{-q}] or [adim], Power Law with exponent q, Coulomb Friction
+  PHX::MDField<const ParamScalarT>      bulkFrictionField;  // t
+  PHX::MDField<const ParamScalarT>      basalDebrisField;   // 
   PHX::MDField<const EffPressureST>     N;                  // [kPa]
   PHX::MDField<const MeshScalarT>       coordVec;           // [km]
   PHX::MDField<const TemperatureST>     flowRate;           // [Pa^{-n} s^{-1}]
@@ -101,13 +107,13 @@ class BasalFrictionCoefficient : public PHX::EvaluatorWithBaseImpl<Traits>,
   bool nodal;
   bool is_side_equation;
   bool is_power_parameter;
-  enum class BETA_TYPE {CONSTANT, FIELD, POWER_LAW, REGULARIZED_COULOMB};
+  enum class BETA_TYPE {CONSTANT, FIELD, POWER_LAW, REGULARIZED_COULOMB, DEBRIS_FRICTION};
   enum class FIELD_TYPE {CONSTANT, FIELD, EXPONENT_OF_FIELD, EXPONENT_OF_FIELD_AT_NODES};
   enum class EFFECTIVE_PRESSURE_TYPE {CONSTANT, FIELD, HYDROSTATIC, HYDROSTATIC_AT_NODES};
   enum class FLOW_RATE_TYPE {CONSTANT, VISCOSITY_FLOW_RATE};
   BETA_TYPE beta_type;
   EFFECTIVE_PRESSURE_TYPE effectivePressure_type;
-  FIELD_TYPE mu_type, bedRoughness_type;
+  FIELD_TYPE mu_type, bedRoughness_type, bulkFriction_type, basalDebris_type;
   FLOW_RATE_TYPE flowRate_type;
 
   PHAL::MDFieldMemoizer<Traits> memoizer;

src/landIce/evaluators/LandIce_BasalFrictionCoefficient_Def.hpp

@@ -106,6 +106,8 @@ void velocity_solver_solve_fo(int nLayers, int globalVerticesStride,
           std::vector<double>& effectivePressureData,
     const std::vector<double>& muData,
     const std::vector<double>& bedRoughnessData,
+    const std::vector<double>& bulkFrictionData,
+    const std::vector<double>& basalDebrisData,
     const std::vector<double>& temperatureDataOnPrisms,
     std::vector<double>& bodyForceMagnitudeOnBasalCell,
     std::vector<double>& dissipationHeatOnPrisms,
@@ -154,6 +156,8 @@ void velocity_solver_solve_fo(int nLayers, int globalVerticesStride,
   STKFieldType* dirichletField = meshStruct->metaData->get_field <double> (stk::topology::NODE_RANK, "dirichlet_field");
   STKFieldType* muField = meshStruct->metaData->get_field <double> (stk::topology::NODE_RANK, "mu");
   STKFieldType* bedRoughnessField = meshStruct->metaData->get_field <double> (stk::topology::NODE_RANK, "bed_roughness");
+  STKFieldType* bulkFrictionField = meshStruct->metaData->get_field <double> (stk::topology::NODE_RANK, "bulk_friction_coefficient");
+  STKFieldType* basalDebrisField = meshStruct->metaData->get_field <double> (stk::topology::NODE_RANK, "basal_debris_factor");
   STKFieldType* stiffeningFactorLogField = meshStruct->metaData->get_field <double> (stk::topology::NODE_RANK, "stiffening_factor_log");
   STKFieldType* effectivePressureField = meshStruct->metaData->get_field <double> (stk::topology::NODE_RANK, "effective_pressure");
   STKFieldType* betaField;
@@ -229,6 +233,16 @@ void velocity_solver_solve_fo(int nLayers, int globalVerticesStride,
       double* bedRoughnessVal = stk::mesh::field_data(*bedRoughnessField, node);
       bedRoughnessVal[0] = bedRoughnessData[ib];
     }
+
+    if ((il == 0) && !bulkFrictionData.empty() && (bulkFrictionField != nullptr)) {
+      double* bulkFrictionVal = stk::mesh::field_data(*bulkFrictionField, node);
+      bulkFrictionVal[0] = bulkFrictionData[ib];
+    }
+
+    if ((il == 0) && !basalDebrisData.empty() && (basalDebrisField != nullptr)) {
+      double* basalDebrisVal = stk::mesh::field_data(*basalDebrisField, node);
+      basalDebrisVal[0] = basalDebrisData[ib];
+    }
   }
 
   //In the following we import the temperature field temperatureDataOnPrisms from MPAS, 
@@ -325,6 +339,7 @@ void velocity_solver_solve_fo(int nLayers, int globalVerticesStride,
   try {
     auto model = slvrfctry->createModel(albanyApp);
     solver = slvrfctry->createSolver(model, Teuchos::null, true);
+    //solver = slvrfctry->createSolver(mpiComm, model, Teuchos::null, true);
 
     Teuchos::ParameterList solveParams;
     solveParams.set("Compute Sensitivities", false);
@@ -492,38 +507,6 @@ void velocity_solver_finalize() {
     Kokkos::finalize();
 }
 
-
-//DEPRECATED
-void velocity_solver_solve_fo(int nLayers, int globalVerticesStride,
-    int globalTrianglesStride, bool ordering, bool first_time_step,
-    const std::vector<int>& indexToVertexID,
-    const std::vector<int>& indexToTriangleID, double minBeta,
-    const std::vector<double>& regulThk,
-    const std::vector<double>& levelsNormalizedThickness,
-    const std::vector<double>& elevationData,
-    const std::vector<double>& thicknessData,
-          std::vector<double>& betaData,
-    const std::vector<double>& bedTopographyData,
-    const std::vector<double>& smbData,
-    const std::vector<double>& stiffeningFactorData,
-    const std::vector<double>& effectivePressureData,
-    const std::vector<double>& muData,
-    const std::vector<double>& temperatureDataOnPrisms,
-    std::vector<double>& bodyForceMagnitudeOnBasalCell,
-    std::vector<double>& dissipationHeatOnPrisms,
-    std::vector<double>& velocityOnVertices,
-    int& error,
-    const double& deltat) {
-      std::vector<double> effectivePressureDataNonConst(effectivePressureData);
-      const std::vector<double> bedRoughnessData;
-      velocity_solver_solve_fo(nLayers, globalVerticesStride, globalTrianglesStride, ordering, first_time_step,
-                               indexToVertexID,indexToTriangleID, minBeta, regulThk, levelsNormalizedThickness, 
-                               elevationData, thicknessData, betaData, bedTopographyData, smbData, stiffeningFactorData,
-                               effectivePressureDataNonConst, muData, bedRoughnessData, temperatureDataOnPrisms, 
-                               bodyForceMagnitudeOnBasalCell,dissipationHeatOnPrisms, velocityOnVertices, error, deltat );
-}
-
-
 /*duality:
  *
  *   mpas(F) |  albany
@@ -642,6 +625,8 @@ void velocity_solver_extrude_3d_grid(int nLayers, int globalTrianglesStride,
   basalFrictionParams.set("Mu Field Name",basalFrictionParams.get("Mu Field Name","mu"));
   basalFrictionParams.set("Bed Roughness Type",basalFrictionParams.get("Bed Roughness Type","Field"));
   basalFrictionParams.set("Effective Pressure Type",basalFrictionParams.get("Effective Pressure Type","Field"));
+  basalFrictionParams.set("Bulk Friction Coefficient Type",basalFrictionParams.get("Bulk Friction Coefficient Type","Field"));
+  basalFrictionParams.set("Basal Debris Factor Type", basalFrictionParams.get("Basal Debris Factor Type","Field"));
   basalFrictionParams.set<bool>("Zero Beta On Floating Ice", zeroBetaOnShelf);
   basalFrictionParams.set<bool>("Zero Effective Pressure On Floating Ice At Nodes", zeroEffectPressOnShelf);
 
@@ -714,7 +699,7 @@ void velocity_solver_extrude_3d_grid(int nLayers, int globalTrianglesStride,
 
   discretizationList->set("Workset Size", discretizationList->get("Workset Size", -1));
 
-  discretizationList->set("Method", discretizationList->get("Method", "STKExtruded")); //set to STKExtruded is not defined
+  discretizationList->set("Method", discretizationList->get("Method", "Extruded")); //set to Extruded is not defined
 
   auto& rfi = discretizationList->sublist("Required Fields Info");
   int fp = rfi.get<int>("Number Of Fields",0);

src/landIce/interfaceWithMPAS/Interface.cpp

@@ -883,6 +883,34 @@ Hydrology::constructEvaluators (PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
   ptr_lambda = Teuchos::rcp(new PHAL::SharedParameter<EvalT,PHAL::AlbanyTraits>(*p,dl));
   fm0.template registerEvaluator<EvalT>(ptr_lambda);
 
+  //--- Shared Parameter for basal friction coefficient: bulk friction coefficient ---//
+  p = Teuchos::rcp(new Teuchos::ParameterList("Bulk Friction Coefficient: bulk_friction"));
+
+  param_name = ParamEnumName::BulkFriction;
+  p->set<std::string>("Parameter Name", param_name);
+  p->set<Teuchos::RCP<Albany::ScalarParameterAccessors<EvalT>>>("Accessors", this->getAccessors()->template at<EvalT>());
+  p->set< Teuchos::RCP<ParamLib> >("Parameter Library", paramLib);
+  p->set<const Teuchos::ParameterList*>("Parameters List", &params->sublist("Parameters"));
+  p->set<double>("Default Nominal Value", params->sublist("LandIce Bulk Friction Coefficient").get<double>(param_name,-1.0));
+
+  Teuchos::RCP<PHAL::SharedParameter<EvalT,PHAL::AlbanyTraits>> ptr_bulk_friction;
+  ptr_bulk_friction = Teuchos::rcp(new PHAL::SharedParameter<EvalT,PHAL::AlbanyTraits>(*p,dl));
+  fm0.template registerEvaluator<EvalT>(ptr_bulk_friction);
+
+  //--- Shared Parameter for basal friction coefficient: basal debris factor ---//
+  p = Teuchos::rcp(new Teuchos::ParameterList("Basal Debris Factor: basal_debris"));
+
+  param_name = ParamEnumName::BasalDebris;
+  p->set<std::string>("Parameter Name", param_name);
+  p->set<Teuchos::RCP<Albany::ScalarParameterAccessors<EvalT>>>("Accessors", this->getAccessors()->template at<EvalT>());
+  p->set< Teuchos::RCP<ParamLib> >("Parameter Library", paramLib);
+  p->set<const Teuchos::ParameterList*>("Parameters List", &params->sublist("Parameters"));
+  p->set<double>("Default Nominal Value", params->sublist("LandIce Basal Debris Factor").get<double>(param_name,-1.0));
+
+  Teuchos::RCP<PHAL::SharedParameter<EvalT,PHAL::AlbanyTraits>> ptr_basal_debris;
+  ptr_basal_debris = Teuchos::rcp(new PHAL::SharedParameter<EvalT,PHAL::AlbanyTraits>(*p,dl));
+  fm0.template registerEvaluator<EvalT>(ptr_basal_debris);
+
   //--- Shared Parameter for basal friction coefficient: mu ---//
   p = Teuchos::rcp(new Teuchos::ParameterList("Basal Friction Coefficient: mu"));
 

src/landIce/problems/LandIce_Hydrology.hpp

@@ -18,6 +18,8 @@ enum class ParamEnum
   BedRoughness,
   Mu,
   Power,
+  BulkFriction,
+  BasalDebris,
   Homotopy,
   GLHomotopy,
   Theta_0,
@@ -31,6 +33,8 @@ namespace ParamEnumName
   static const std::string BedRoughness    = "Bed Roughness";
   static const std::string Mu              = "Mu";
   static const std::string Power           = "Power Exponent";
+  static const std::string BulkFriction    = "Bulk Friction Coefficient";
+  static const std::string BasalDebris     = "Basal Debris Factor";
   static const std::string HomotopyParam   = "Homotopy Parameter";
   static const std::string GLHomotopyParam = "Glen's Law Homotopy Parameter";
   static const std::string theta_0         = "Theta 0"; 

src/landIce/problems/LandIce_ParamEnum.hpp

@@ -100,6 +100,8 @@ StokesFOBase (const Teuchos::RCP<Teuchos::ParameterList>& params_,
   bed_topography_param_name   = bed_topography_name + "_param";
   bed_topography_observed_name= "observed_" + bed_topography_name;
   bed_roughness_name          = params->sublist("Variables Names").get<std::string>("Bed Roughness Name"    ,"bed_roughness");
+  bulk_friction_name          = params->sublist("Variables Names").get<std::string>("Bulk Friction Coefficient Name", "bulk_friction");
+  basal_debris_name           = params->sublist("Variables Names").get<std::string>("Basal Debris Factor Name", "basal_debris");
   flow_factor_name            = params->sublist("Variables Names").get<std::string>("Flow Factor Name"       ,"flow_factor");
   stiffening_factor_log_name  = params->sublist("Variables Names").get<std::string>("Stiffening Factor Log Name" ,"stiffening_factor_log");
   damage_factor_name          = params->sublist("Variables Names").get<std::string>("Damage Factor Name"     ,"damage_factor");
@@ -137,7 +139,7 @@ void StokesFOBase::buildProblem (Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpec
   std::string tensorCubDegName = "Cubature Degrees (Horiz Vert)";
   if(tensorProductCell && params->isParameter(tensorCubDegName)) {
     TEUCHOS_TEST_FOR_EXCEPTION (params->isParameter("Cubature Degree") && params->isParameter(tensorCubDegName), std::logic_error,
-                                  "Error! Either provide parameter 'Cubatur Degree' or 'Cubature Degrees (Horiz Vert)', not both\n");
+                                  "Error! Either provide parameter 'Cubature Degree' or 'Cubature Degrees (Horiz Vert)', not both\n");
     Teuchos::Array<int> tensorCubDegrees;
     tensorCubDegrees = params->get<Teuchos::Array<int>>(tensorCubDegName);
     if(cellType->getKey() == shards::Wedge<6>::key)
@@ -149,9 +151,7 @@ void StokesFOBase::buildProblem (Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpec
       cellCubature = cubFactory.create<PHX::Device, RealType, RealType>(*cellType, degree);
     }    
   } else {
-    TEUCHOS_TEST_FOR_EXCEPTION (!tensorProductCell && params->isParameter(tensorCubDegName), std::logic_error,
-                                  "Error! Parameter 'Cubature Degrees (Horiz Vert)' should be provided only for extruded meshes with Wedge or Hexaedron cells\n"
-                                  "        Use 'Cubature Degree instead'");
+    //TEUCHOS_TEST_FOR_EXCEPTION (!tensorProductCell && params->isParameter(tensorCubDegName), std::logic_error, "Error! Parameter 'Cubature Degrees (Horiz Vert)' should be provided only for extruded meshes with Wedge or Hexaedron cells\n""        Use 'Cubature Degree instead'");
     int cubDegree = params->get("Cubature Degree", defaultCubDegree);
     cellCubature = cubFactory.create<PHX::Device, RealType, RealType>(*cellType, cubDegree);
   }  
@@ -639,6 +639,7 @@ void StokesFOBase::setupEvaluatorRequests ()
       } else {
         // We have a sliding law, which requires a mu (possibly a field),
         // and possibly a bed roughness
+	// debris friction also requires bulk friction coefficient + basal debris factor
         const auto mu_type = util::upper_case(bfc.get<std::string>("Mu Type"));
         if (mu_type!="CONSTANT") {
           auto mu_field_name = bfc.get<std::string>("Mu Field Name");
@@ -658,6 +659,28 @@ void StokesFOBase::setupEvaluatorRequests ()
           }
           setSingleFieldProperties(fname, FRT::Scalar, FST::ParamScalar);
         }
+	if (type=="DEBRIS FRICTION") {
+	  // For debris friction slip law, we have bed roughness as in reg Coulomb
+	  // As well as two new parameters: bulk friction coefficient and basal debris factor
+	  auto fname = "bed_roughness";
+	  if (is_input_field[fname] || is_ss_input_field[ssName][fname]) {
+	    ss_build_interp_ev[ssName][fname][IReq::CELL_TO_SIDE] = true;
+	    ss_build_interp_ev[ssName][fname][IReq::QP_VAL] = true;
+	  }
+	  setSingleFieldProperties(fname, FRT::Scalar, FST::ParamScalar);
+	  auto bname = "bulk_friction";
+	  if (is_input_field[bname] || is_ss_input_field[ssName][bname]) {
+	    ss_build_interp_ev[ssName][bname][IReq::CELL_TO_SIDE] = true;
+	    ss_build_interp_ev[ssName][bname][IReq::QP_VAL] = true;
+	  }
+	  setSingleFieldProperties(bname, FRT::Scalar, FST::ParamScalar);
+	  auto dname = "basal_debris";
+	  if (is_input_field[dname] || is_ss_input_field[ssName][dname]) {
+	    ss_build_interp_ev[ssName][dname][IReq::CELL_TO_SIDE] = true;
+	    ss_build_interp_ev[ssName][dname][IReq::QP_VAL] = true;
+	  }
+	  setSingleFieldProperties(dname, FRT::Scalar, FST::ParamScalar);
+	}
       }
     }
 

src/landIce/problems/LandIce_StokesFOBase.cpp

@@ -321,6 +321,8 @@ class StokesFOBase : public Albany::AbstractProblem {
   std::string effective_pressure_name;
   std::string basal_friction_name;
   std::string bed_roughness_name;
+  std::string bulk_friction_name;
+  std::string basal_debris_name;
   std::string sliding_velocity_name;
   std::string vertically_averaged_velocity_name;
 
@@ -1591,6 +1593,49 @@ void StokesFOBase::constructBasalBCEvaluators (PHX::FieldManager<PHAL::AlbanyTra
       fm0.template registerEvaluator<EvalT>(ptr_lambda);
     }
 
+  //--- Shared Parameter for basal friction coefficient: bulk friction coefficient ---//
+  //p = Teuchos::rcp(new Teuchos::ParameterList("Bulk Friction Coefficient: bulk_friction"));
+
+  //param_name = ParamEnumName::BulkFriction;
+  //p->set<std::string>("Parameter Name", param_name);
+  //p->set<Teuchos::RCP<Albany::ScalarParameterAccessors<EvalT>>>("Accessors", this->getAccessors()->template at<EvalT>());
+  //p->set< Teuchos::RCP<ParamLib> >("Parameter Library", paramLib);
+  //p->set<const Teuchos::ParameterList*>("Parameters List", &params->sublist("Parameters"));
+  //p->set<double>("Default Nominal Value", params->sublist("LandIce Bulk Friction Coefficient").get<double>(param_name,-1.0));
+
+  //Teuchos::RCP<PHAL::SharedParameter<EvalT,PHAL::AlbanyTraits>> ptr_bulk_friction;
+    
+    //--- Shared Parameter for bulk friction coefficient: bulk_friction ---//
+    param_name = "Bulk Friction Coefficient";
+
+    if(!PHAL::is_field_evaluated<EvalT>(fm0, param_name, dl->shared_param)) {
+      p = Teuchos::rcp(new Teuchos::ParameterList("Bulk Friction Coefficient: bulk_friction"));
+      p->set<std::string>("Parameter Name", param_name);
+      p->set<Teuchos::RCP<Albany::ScalarParameterAccessors<EvalT>>>("Accessors", this->getAccessors()->template at<EvalT>());
+      p->set< Teuchos::RCP<ParamLib> >("Parameter Library", paramLib);
+      p->set<const Teuchos::ParameterList*>("Parameters List", &params->sublist("Parameters"));
+      p->set<double>("Default Nominal Value", pl->sublist("Bulk Friction Coefficient").get<double>(param_name,-1.0));
+
+      Teuchos::RCP<PHAL::SharedParameter<EvalT,PHAL::AlbanyTraits>> ptr_bulk_friction;
+      ptr_bulk_friction = Teuchos::rcp(new PHAL::SharedParameter<EvalT,PHAL::AlbanyTraits>(*p,dl));
+      fm0.template registerEvaluator<EvalT>(ptr_bulk_friction);
+    }
+
+    //--- Shared Parameter for basal debris factor: basal_debris ---//
+    param_name = "Basal Debris Factor";
+
+    if(!PHAL::is_field_evaluated<EvalT>(fm0, param_name, dl->shared_param)) {
+      p = Teuchos::rcp(new Teuchos::ParameterList("Basal Debris Factor: basal_debris"));
+      p->set<std::string>("Parameter Name", param_name);
+      p->set<Teuchos::RCP<Albany::ScalarParameterAccessors<EvalT>>>("Accessors", this->getAccessors()->template at<EvalT>());
+      p->set< Teuchos::RCP<ParamLib> >("Parameter Library", paramLib);
+      p->set<const Teuchos::ParameterList*>("Parameters List", &params->sublist("Parameters"));
+      p->set<double>("Default Nominal Value", pl->sublist("Basal Debris Factor").get<double>(param_name,-1.0));
+
+      Teuchos::RCP<PHAL::SharedParameter<EvalT,PHAL::AlbanyTraits>> ptr_basal_debris;
+      ptr_basal_debris = Teuchos::rcp(new PHAL::SharedParameter<EvalT,PHAL::AlbanyTraits>(*p,dl));
+      fm0.template registerEvaluator<EvalT>(ptr_basal_debris);
+    }
     //--- Shared Parameter for basal friction coefficient: muPowerLaw ---//
     param_name = ParamEnumName::Mu;