Skip to content

Debris Friction sliding law #1184

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mperego merged 4 commits into from
Mar 4, 2026
Merged

Debris Friction sliding law #1184

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
ea2d329
Debris-bed friction slip law
jeremy-brooks-1 Dec 20, 2024
544532d
LandIce: fixes and improvements for the debris friction sliding law
mperego Dec 17, 2025
326b291
LandIce: Fix naming of Basal Debris Factor field.
mperego Feb 26, 2026
0def4ae
LandIce: Improve logic in the MPAS Interface
mperego Mar 3, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 13 additions & 5 deletions src/landIce/evaluators/LandIce_BasalFrictionCoefficient.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
PHX::MDField<const ScalarT,Dim> bulkFrictionParam; // [adim], bulk friction coefficient (for DEBRIS_FRICTION only)
PHX::MDField<const ScalarT,Dim> basalDebrisParam; // [kPa 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)
Expand All @@ -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; // [adim], bulk friction coefficient (for DEBRIS_FRICTION only)
PHX::MDField<const ParamScalarT> basalDebrisField; // [kPa m^-1 s] basal debris factor (for DEBRIS_FRICTION only)
PHX::MDField<const EffPressureST> N; // [kPa]
PHX::MDField<const MeshScalarT> coordVec; // [km]
PHX::MDField<const TemperatureST> flowRate; // [Pa^{-n} s^{-1}]
Expand All @@ -91,6 +97,8 @@ class BasalFrictionCoefficient : public PHX::EvaluatorWithBaseImpl<Traits>,
bool save_pressure_field;
double overburden_fraction; // [adim]
double pressure_smoothing_length_scale; //[km]
double N0; // [kPa] Effective Pressure Regularization;
double u0; // [m/yr] Sliding Velocity Regularization;

int numNodes;
int numQPs;
Expand All @@ -101,13 +109,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;
Expand Down
214 changes: 201 additions & 13 deletions src/landIce/evaluators/LandIce_BasalFrictionCoefficient_Def.hpp
View file
Open in desktop

Large diffs are not rendered by default.

354 changes: 234 additions & 120 deletions src/landIce/interfaceWithMPAS/Interface.cpp
View file
Open in desktop

Large diffs are not rendered by default.

28 changes: 28 additions & 0 deletions src/landIce/problems/LandIce_Hydrology.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"));

Expand Down
4 changes: 4 additions & 0 deletions src/landIce/problems/LandIce_ParamEnum.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,8 @@ enum class ParamEnum
BedRoughness,
Mu,
Power,
BulkFriction,
BasalDebris,
Homotopy,
GLHomotopy,
Theta_0,
Expand All @@ -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";
Expand Down
29 changes: 27 additions & 2 deletions src/landIce/problems/LandIce_StokesFOBase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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_factor");
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");
Expand Down Expand Up @@ -138,7 +140,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)
Expand Down Expand Up @@ -643,6 +645,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 and 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");
Expand All @@ -655,13 +658,35 @@ void StokesFOBase::setupEvaluatorRequests ()
if (type=="REGULARIZED COULOMB") {
// For Coulomb we *may* have another distributed parameter field.
// We interpolate (and possibly project on ss) only if it is an inputs
auto fname = "bed_roughness";
auto fname = bed_roughness_name;
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);
}
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_name;
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_name;
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_name;
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);
}
}
}

Expand Down
43 changes: 42 additions & 1 deletion src/landIce/problems/LandIce_StokesFOBase.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -1456,7 +1458,9 @@ void StokesFOBase::constructBasalBCEvaluators (PHX::FieldManager<PHAL::AlbanyTra
std::string ice_thickness_side_name = side_fname(ice_thickness_name, ssName);
std::string ice_overburden_side_name = side_fname("ice_overburden", ssName);
std::string effective_pressure_side_name = side_fname(effective_pressure_name, ssName);
std::string bed_roughness_side_name = side_fname("bed_roughness", ssName);
std::string bed_roughness_side_name = side_fname(bed_roughness_name, ssName);
std::string bulk_friction_side_name = side_fname(bulk_friction_name, ssName);
std::string basal_debris_side_name = side_fname(basal_debris_name, ssName);
std::string bed_topography_side_name = side_fname(bed_topography_name, ssName);
std::string flow_factor_side_name = side_fname(flow_factor_name, ssName);

Expand Down Expand Up @@ -1591,6 +1595,39 @@ void StokesFOBase::constructBasalBCEvaluators (PHX::FieldManager<PHAL::AlbanyTra
fm0.template registerEvaluator<EvalT>(ptr_lambda);
}

//--- 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"));
Teuchos::ParameterList beta_list = pl->sublist("Basal Friction Coefficient");
p->set<double>("Default Nominal Value", beta_list.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"));
Teuchos::ParameterList beta_list = pl->sublist("Basal Friction Coefficient");
p->set<double>("Default Nominal Value", beta_list.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;

Expand Down Expand Up @@ -1635,6 +1672,8 @@ void StokesFOBase::constructBasalBCEvaluators (PHX::FieldManager<PHAL::AlbanyTra
p->set<std::string>("Effective Pressure QP Variable Name", effective_pressure_side_name);
p->set<std::string>("Flow Rate Variable Name", side_fname(flow_factor_name, ssName),flow_factor_side_name);
p->set<std::string>("Bed Roughness Variable Name", bed_roughness_side_name);
p->set<std::string>("Basal Debris Factor Variable Name", basal_debris_side_name);
p->set<std::string>("Bulk Friction Coefficient Variable Name", bulk_friction_side_name);
p->set<std::string>("Side Set Name", ssName);
p->set<std::string>("Coordinate Vector Variable Name", side_fname(Albany::coord_vec_name, ssName));
p->set<Teuchos::ParameterList*>("Parameter List", &pl->sublist("Basal Friction Coefficient"));
Expand Down Expand Up @@ -1751,6 +1790,8 @@ void StokesFOBase::constructSMBEvaluators (PHX::FieldManager<PHAL::AlbanyTraits>
std::string effective_pressure_side_name = basal_fname(effective_pressure_name);
std::string vertically_averaged_velocity_side_name = basal_fname(vertically_averaged_velocity_name);
std::string bed_roughness_side_name = basal_fname(bed_roughness_name);
std::string bulk_friction_side_name = basal_fname(bulk_friction_name);
std::string basal_debris_side_name = basal_fname(basal_debris_name);

// -------------------------------- LandIce evaluators ------------------------- //

Expand Down
15 changes: 15 additions & 0 deletions tests/landIce/FO_GIS/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -250,6 +250,10 @@ if (ALBANY_MUELU)
${CMAKE_CURRENT_BINARY_DIR}/input_fo_humboldt_depthInt_regCoulomb_constA.yaml)
add_test(${testNameRoot}_Humboldt_DepthIntegrated_RegCoulomb_ConstA ${Albany.exe} input_fo_humboldt_depthInt_regCoulomb_constA.yaml)

configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input_fo_humboldt_depthInt_debrisFriction.yaml
${CMAKE_CURRENT_BINARY_DIR}/input_fo_humboldt_depthInt_debrisFriction.yaml)
add_test(${testNameRoot}_Humboldt_DepthIntegrated_DebrisFriction ${Albany.exe} input_fo_humboldt_depthInt_debrisFriction.yaml)

if (NOT ALBANY_PARALELL_EXODUS)
set_tests_properties(${testNameRoot}_Humboldt_MueLu
PROPERTIES
Expand All @@ -260,14 +264,21 @@ if (ALBANY_MUELU)
PROPERTIES
LABELS "LandIce;Forward"
FIXTURES_REQUIRED humboldtMeshSetup)

set_tests_properties(${testNameRoot}_Humboldt_DepthIntegrated_RegCoulomb
PROPERTIES
LABELS "LandIce;Forward"
FIXTURES_REQUIRED humboldtMeshSetup)

set_tests_properties(${testNameRoot}_Humboldt_DepthIntegrated_RegCoulomb_ConstA
PROPERTIES
LABELS "LandIce;Forward"
FIXTURES_REQUIRED humboldtMeshSetup)

set_tests_properties(${testNameRoot}_Humboldt_DepthIntegrated_DebrisFriction
PROPERTIES
LABELS "LandIce;Forward"
FIXTURES_REQUIRED humboldtMeshSetup)
else ()
set_tests_properties(${testNameRoot}_Humboldt_MueLu
PROPERTIES
Expand All @@ -284,6 +295,10 @@ if (ALBANY_MUELU)
set_tests_properties(${testNameRoot}_Humboldt_DepthIntegrated_RegCoulomb_ConstA
PROPERTIES
LABELS "LandIce;Forward")

set_tests_properties(${testNameRoot}_Humboldt_DepthIntegrated_DebrisFriction
PROPERTIES
LABELS "LandIce;Forward")
endif()

endif()
Expand Down
Loading
Loading