Skip to content

Commit

Permalink
Debris-bed friction slip law
Browse files Browse the repository at this point in the history 
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.
  • Loading branch information
@jeremy-brooks-1
jeremy-brooks-1 committed Jan 9, 2025
1 parent 2e14b42 commit 9badc68
Show file tree
Hide file tree
Showing 7 changed files with 331 additions and 60 deletions.
18 changes: 12 additions & 6 deletions src/landIce/evaluators/LandIce_BasalFrictionCoefficient.hpp
View file
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 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)
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; // 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}]
Expand Down Expand Up @@ -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;
Expand Down
214 changes: 197 additions & 17 deletions src/landIce/evaluators/LandIce_BasalFrictionCoefficient_Def.hpp
View file

Large diffs are not rendered by default.

51 changes: 18 additions & 33 deletions src/landIce/interfaceWithMPAS/Interface.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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);

Expand Down Expand Up @@ -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);
Expand Down
28 changes: 28 additions & 0 deletions src/landIce/problems/LandIce_Hydrology.hpp
View file
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
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
31 changes: 27 additions & 4 deletions src/landIce/problems/LandIce_StokesFOBase.cpp
View file
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");
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 @@ -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)
Expand All @@ -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);
}
Expand Down Expand Up @@ -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");
Expand All @@ -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);
}
}
}

Expand Down
45 changes: 45 additions & 0 deletions src/landIce/problems/LandIce_StokesFOBase.hpp
View file
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 @@ -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;

Expand Down

0 comments on commit 9badc68

Please sign in to comment.