From cb2b94db2e307ae73503557a12b7774fb6ae468d Mon Sep 17 00:00:00 2001 From: Jeremy Brooks Date: Fri, 20 Dec 2024 12:04:58 -0800 Subject: [PATCH] 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. --- .../serial/pm_cpu_gnu_modules.sh | 4 +- .../LandIce_BasalFrictionCoefficient.hpp | 18 +- .../LandIce_BasalFrictionCoefficient_Def.hpp | 214 ++++++++++++++++-- src/landIce/interfaceWithMPAS/Interface.cpp | 51 ++--- src/landIce/problems/LandIce_Hydrology.hpp | 28 +++ src/landIce/problems/LandIce_ParamEnum.hpp | 4 + src/landIce/problems/LandIce_StokesFOBase.cpp | 31 ++- src/landIce/problems/LandIce_StokesFOBase.hpp | 45 ++++ 8 files changed, 333 insertions(+), 62 deletions(-) diff --git a/doc/dashboards/perlmutter.nersc.gov/serial/pm_cpu_gnu_modules.sh b/doc/dashboards/perlmutter.nersc.gov/serial/pm_cpu_gnu_modules.sh index 57cca012c8..846707aa97 100644 --- a/doc/dashboards/perlmutter.nersc.gov/serial/pm_cpu_gnu_modules.sh +++ b/doc/dashboards/perlmutter.nersc.gov/serial/pm_cpu_gnu_modules.sh @@ -11,9 +11,9 @@ module load cray-netcdf-hdf5parallel/4.9.0.9 module load cray-parallel-netcdf/1.12.3.9 module load cmake/3.24.3 -module load cray-python/3.9.13.1 +#module load cray-python/3.9.13.1 -module load e4s/23.05 +module load e4s/23.08 spack env activate -V gcc spack load superlu diff --git a/src/landIce/evaluators/LandIce_BasalFrictionCoefficient.hpp b/src/landIce/evaluators/LandIce_BasalFrictionCoefficient.hpp index 0f1ba534ec..de39b63658 100644 --- a/src/landIce/evaluators/LandIce_BasalFrictionCoefficient.hpp +++ b/src/landIce/evaluators/LandIce_BasalFrictionCoefficient.hpp @@ -48,14 +48,18 @@ class BasalFrictionCoefficient : public PHX::EvaluatorWithBaseImpl, // Coefficients for computing beta (if not given) double n; // [adim] exponent of Glen's law PHX::MDField muParam; // [yr^q m^{-q}], friction coefficient of the power Law with exponent q - PHX::MDField bedRoughnessParam; // [km], Bed bumps avg length divided by bed bumps avg slope (for REGULARIZED_COULOMB only) - PHX::MDField powerParam; // [adim], Exponent (for POWER_LAW and REGULARIZED COULOMB only) - + PHX::MDField bedRoughnessParam; // [km], Bed bumps avg length divided by bed bumps avg slope (for REGULARIZED_COULOMB and DEBRIS_FRICTION only) + PHX::MDField powerParam; // [adim], Exponent (for POWER_LAW and REGULARIZED COULOMB and DEBRIS_FRICTION only) + PHX::MDField bulkFrictionParam; // [adim], bulk friction coefficient (for DEBRIS_FRICTION only) + PHX::MDField 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, PHX::MDField u_norm; // [m yr^{-1}] PHX::MDField bedRoughnessField; // [km], characteristic length PHX::MDField muField; // [yr^q m^{-q}] or [adim], Power Law with exponent q, Coulomb Friction + PHX::MDField bulkFrictionField; // t + PHX::MDField basalDebrisField; // PHX::MDField N; // [kPa] PHX::MDField coordVec; // [km] PHX::MDField flowRate; // [Pa^{-n} s^{-1}] @@ -101,13 +107,13 @@ class BasalFrictionCoefficient : public PHX::EvaluatorWithBaseImpl, 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 memoizer; diff --git a/src/landIce/evaluators/LandIce_BasalFrictionCoefficient_Def.hpp b/src/landIce/evaluators/LandIce_BasalFrictionCoefficient_Def.hpp index 817406a8ea..f402391556 100644 --- a/src/landIce/evaluators/LandIce_BasalFrictionCoefficient_Def.hpp +++ b/src/landIce/evaluators/LandIce_BasalFrictionCoefficient_Def.hpp @@ -51,7 +51,7 @@ BasalFrictionCoefficient (const Teuchos::ParameterList& p, //Validate Parameters Teuchos::ParameterList validPL; - validPL.set("Type", "Power Law", "Type Of Beta: Constant, Power Law, Regularized Coulomb"); + validPL.set("Type", "Power Law", "Type Of Beta: Constant, Power Law, Regularized Coulomb, Debris Friction"); validPL.set("Power Exponent", 1.0, "Exponent of Power Law or Regularized Coulomb Law"); validPL.set("Mu Type", "Field", "Mu Type: Constant,Field, Exponent Of Field, Exponent Of Field At Nodes"); validPL.set("Mu Field Name", "", "Name of the Field Mu"); @@ -69,7 +69,11 @@ BasalFrictionCoefficient (const Teuchos::ParameterList& p, validPL.set("Zero Effective Pressure On Floating Ice At Nodes", false, "Whether to zero the effective pressure on floating ice at nodes"); validPL.set("Zero Beta On Floating Ice", false, "Whether to zero beta on floating ice"); validPL.set("Exponentiate Scalar Parameters", false, "Whether the scalar parameters needs to be exponentiate"); - validPL.set("Use Pressurized Bed Above Sea Level", false, "Whether to use a Downs & Johnson (2022) type parameterization for basal water pressure"); + validPL.set("Use Pressurized Bed Above Sea Level", false, "Whether to use a Downs & Johnson (2022) type parameterization for basal water pressure"); + validPL.set("Bulk Friction Coefficient Type", "Constant", "Bulk Friction Coefficient Type: Field"); + validPL.set("Bulk Friction Coefficient", 0.0, "Constant value for Bulk Friction Coefficient (for debris friction slip)"); + validPL.set("Basal Debris Factor Type", "Constant", "Basal Debris Factor Type: Field"); + validPL.set("Basal Debris Factor", 0.0, "Constant value for Basal Debris Factor (for debris friction slip)"); beta_list.validateParameters(validPL,0); zero_on_floating = beta_list.get ("Zero Beta On Floating Ice", false); @@ -108,6 +112,8 @@ BasalFrictionCoefficient (const Teuchos::ParameterList& p, printedMu = -9999.999; printedBedRoughness = -9999.999; printedQ = -9999.999; + printedBulkFriction = -9999.999; + printedBasalDebris = -9999.999; if (betaType == "CONSTANT") { beta_type = BETA_TYPE::CONSTANT; @@ -164,6 +170,31 @@ BasalFrictionCoefficient (const Teuchos::ParameterList& p, std::endl << "Error in LandIce::BasalFrictionCoefficient: \"" << flowRateType << "\" is not a valid parameter for Is not a valid parameter for Flow Rate Type\n"); } + } else if (betaType == "DEBRIS FRICTION") { + beta_type = BETA_TYPE::DEBRIS_FRICTION; +#ifdef OUTPUT_TO_SCREEN + *output << "Velocity-dependent beta (debris friction law):\n\n" + << " beta = mu * N * |u|^{q-1} / [|u| + bedRoughness*A*N^n]^q + BulkFrictionCoefficient*N/|u| + BasalDebrisFactor*|u|\n\n" + << " with N being the effective pressure, |u| the sliding velocity\n"; +#endif + + std::string flowRateType = util::upper_case(beta_list.get("Flow Rate Type", "Viscosity Flow Rate")); + n = p.get("Viscosity Parameter List")->get("Glen's Law n"); + if(flowRateType == "VISCOSITY FLOW RATE") { + flowRate_type = FLOW_RATE_TYPE::VISCOSITY_FLOW_RATE; + flowRate = PHX::MDField(p.get("Flow Rate Variable Name"), dl->cell_scalar2); + this->addDependentField (flowRate); + } else if (flowRateType == "CONSTANT") { //"Constant" + flowRate_type = FLOW_RATE_TYPE::CONSTANT; + flowRate_val = beta_list.get("Flow Rate"); +#ifdef OUTPUT_TO_SCREEN + *output << "LandIce::BasalFrictionCoefficient: Constant Flow Rate value = " << flowRate_val << " (loaded from yaml input file).\n"; +#endif + } else { + TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, + std::endl << "Error in LandIce::BasalFrictionCoefficient: \"" << flowRateType << "\" is not a valid parameter for Is not a valid parameter for Flow Rate Type\n"); + } + } else { TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, std::endl << "Error in LandIce::BasalFrictionCoefficient: \"" << betaType << "\" is not a valid parameter for Beta Type\n"); @@ -186,7 +217,6 @@ BasalFrictionCoefficient (const Teuchos::ParameterList& p, } else if (effectivePressureType == "HYDROSTATIC") { effectivePressure_type = EFFECTIVE_PRESSURE_TYPE::HYDROSTATIC; use_pressurized_bed = beta_list.get("Use Pressurized Bed Above Sea Level", false); - save_pressure_field = p.isParameter("Effective Pressure Output Variable Name"); if(save_pressure_field && nodal) { outN = PHX::MDField(p.get ("Effective Pressure Output Variable Name"), nodal_layout); this->addEvaluatedField (outN); @@ -303,6 +333,80 @@ BasalFrictionCoefficient (const Teuchos::ParameterList& p, u_norm = PHX::MDField(p.get ("Sliding Velocity Variable Name"), layout); this->addDependentField (u_norm); + } else if (beta_type == BETA_TYPE::DEBRIS_FRICTION) { + + std::string bedRoughnessType = util::upper_case((beta_list.isParameter("Bed Roughness Type") ? beta_list.get("Bed Roughness Type") : "Field")); + + auto layout_bedRoughness_field = layout; + if(bedRoughnessType == "CONSTANT") { + bedRoughness_type = FIELD_TYPE::CONSTANT; + } else if (bedRoughnessType == "FIELD") { + bedRoughness_type = FIELD_TYPE::FIELD; + } else if (bedRoughnessType == "EXPONENT OF FIELD AT NODES") { + bedRoughness_type = FIELD_TYPE::EXPONENT_OF_FIELD_AT_NODES; + if(!nodal) { + BF = PHX::MDField(p.get ("BF Variable Name"), dl->node_qp_scalar); + layout_bedRoughness_field = nodal_layout; + this->addDependentField (BF); + } + } else if (bedRoughnessType == "EXPONENT OF FIELD") { + bedRoughness_type = FIELD_TYPE::EXPONENT_OF_FIELD; + } else { + TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, + std::endl << "Error in LandIce::BasalFrictionCoefficient: \"" << bedRoughnessType << "\" is not a valid parameter for Bed Roughness Type\n"); + } + + if(bedRoughness_type == FIELD_TYPE::CONSTANT) { + bedRoughnessParam = PHX::MDField("Bed Roughness", dl->shared_param); + this->addDependentField (bedRoughnessParam); + } else { + bedRoughnessField = PHX::MDField(p.get ("Bed Roughness Variable Name"), layout_bedRoughness_field); + this->addDependentField (bedRoughnessField); + } + + std::string bulkFrictionType = util::upper_case((beta_list.isParameter("Bulk Friction Coefficient Type") ? beta_list.get("Bulk Friction Coefficient Type") : "Field")); + + auto layout_bulkFriction_field = layout; + if (bulkFrictionType == "CONSTANT") { + bulkFriction_type = FIELD_TYPE::CONSTANT; + } else if (bulkFrictionType == "FIELD") { + bulkFriction_type = FIELD_TYPE::FIELD; + } else { + TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, + std::endl << "Error in LandIce::BasalFrictionCoefficient: \"" << bulkFrictionType << "\" is not a valid parameter for Bulk Friction Coefficient Type\n"); + } + + if (bulkFriction_type == FIELD_TYPE::CONSTANT) { + bulkFrictionParam = PHX::MDField("Bulk Friction Coefficient", dl->shared_param); + this->addDependentField (bulkFrictionParam); + } else { + bulkFrictionField = PHX::MDField(p.get ("Bulk Friction Coefficient Variable Name"), layout_bulkFriction_field); + this->addDependentField (bulkFrictionField); + } + + std::string basalDebrisType = util::upper_case((beta_list.isParameter("Basal Debris Factor Type") ? beta_list.get("Basal Debris Factor Type") : "Field")); + + auto layout_basalDebris_field = layout; + if (basalDebrisType == "CONSTANT") { + basalDebris_type = FIELD_TYPE::CONSTANT; + } else if (basalDebrisType == "FIELD") { + basalDebris_type = FIELD_TYPE::FIELD; + } else { + TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, + std::endl << "Error in LandIce:BasalFrictionCoefficient: \"" << basalDebrisType << "\" is not a valid parameter for Basal Debris Factor Type\n"); + } + + if (basalDebris_type == FIELD_TYPE::CONSTANT) { + basalDebrisParam = PHX::MDField("Basal Debris Factor", dl->shared_param); + this->addDependentField (basalDebrisParam); + } else { + basalDebrisField = PHX::MDField(p.get ("Basal Debris Factor Variable Name"), layout_basalDebris_field); + this->addDependentField (basalDebrisField); + } + + u_norm = PHX::MDField(p.get ("Sliding Velocity Variable Name"), layout); + this->addDependentField (u_norm); + } else if (beta_type == BETA_TYPE::POWER_LAW) { auto paramLib = p.get >("Parameter Library"); is_power_parameter = paramLib->isParameter("Power Exponent"); @@ -369,11 +473,11 @@ postRegistrationSetup (typename Traits::SetupData d, template KOKKOS_INLINE_FUNCTION void BasalFrictionCoefficient:: -operator() (const BasalFrictionCoefficient_Tag&, const int& cell) const { +operator() (const BasalFrictionCoefficient_Tag& tag, const int& cell) const { - ParamScalarT mu, bedRoughness, power; + ParamScalarT mu, bedRoughness, power, bulkFriction, basalDebris; - if ((beta_type == BETA_TYPE::POWER_LAW)||(beta_type==BETA_TYPE::REGULARIZED_COULOMB) ) { + if ((beta_type == BETA_TYPE::POWER_LAW) || (beta_type==BETA_TYPE::REGULARIZED_COULOMB) || (beta_type==BETA_TYPE::DEBRIS_FRICTION) ) { if (logParameters) { power = std::exp(Albany::convertScalar(powerParam(0))); if (mu_type == FIELD_TYPE::CONSTANT) @@ -385,7 +489,7 @@ operator() (const BasalFrictionCoefficient_Tag&, const int& cell) const { } } - if (beta_type== BETA_TYPE::REGULARIZED_COULOMB) { + if ( (beta_type== BETA_TYPE::REGULARIZED_COULOMB) || (beta_type == BETA_TYPE::DEBRIS_FRICTION) ) { if (logParameters) { if (bedRoughness_type == FIELD_TYPE::CONSTANT) bedRoughness = std::exp(Albany::convertScalar(bedRoughnessParam(0))); @@ -395,6 +499,13 @@ operator() (const BasalFrictionCoefficient_Tag&, const int& cell) const { } } + if (beta_type == BETA_TYPE::DEBRIS_FRICTION) { + if (bulkFriction_type == FIELD_TYPE::CONSTANT) + bulkFriction = Albany::convertScalar(bulkFrictionParam(0)); + if (basalDebris_type == FIELD_TYPE::CONSTANT) + basalDebris = Albany::convertScalar(basalDebrisParam(0)); + } + ParamScalarT muValue = 1.0; typename Albany::StrongestScalarType::type NVal = N_val; @@ -436,7 +547,7 @@ operator() (const BasalFrictionCoefficient_Tag&, const int& cell) const { thickness_field(cell,ipt)*f_p) + (1.0 - f_p)* KU::max(-1.0 * rho_w*bed_topo_field(cell,ipt),0.0) ),0.0); if(save_pressure_field) { - outN(cell,ipt) = Albany::convertScalar(NVal); + outN(cell,ipt) = static_cast(NVal); } } else { MeshScalarT thickness(0), bed_topo(0); @@ -457,7 +568,7 @@ operator() (const BasalFrictionCoefficient_Tag&, const int& cell) const { thickness_field(cell,ipt)*f_p) + (1.0 - f_p)* KU::max(-1.0 * rho_w*bed_topo_field(cell,ipt),0.0) ),0.0); if(save_pressure_field) { - outN(cell,ipt) = Albany::convertScalar(NVal); + outN(cell,ipt) = static_cast(NVal); } } else { NVal = 0; @@ -533,8 +644,65 @@ operator() (const BasalFrictionCoefficient_Tag&, const int& cell) const { } } - } - } + if (beta_type == BETA_TYPE::DEBRIS_FRICTION) { + ParamScalarT bedRoughnessValue; + ParamScalarT bulkFrictionValue; + ParamScalarT basalDebrisValue; + switch (bedRoughness_type) { + case FIELD_TYPE::FIELD: + bedRoughnessValue = bedRoughnessField(cell,ipt); + break; + case FIELD_TYPE::EXPONENT_OF_FIELD_AT_NODES: + if(nodal) + bedRoughnessValue = std::exp(bedRoughnessField(cell,ipt)); + else { + bedRoughnessValue = 0; + for (int node=0; node void BasalFrictionCoefficient:: @@ -573,7 +740,7 @@ evaluateFields (typename Traits::EvalData workset) if (memoizer.have_saved_data(workset,this->evaluatedFields())) return; - if ( (beta_type == BETA_TYPE::POWER_LAW) || (beta_type == BETA_TYPE::REGULARIZED_COULOMB) ) { + if ( (beta_type == BETA_TYPE::POWER_LAW) || (beta_type == BETA_TYPE::REGULARIZED_COULOMB) || (beta_type == BETA_TYPE::DEBRIS_FRICTION) ) { if (logParameters) { auto hostPower_view = Kokkos::create_mirror_view(powerParam.get_view()); Kokkos::deep_copy(hostPower_view, powerParam.get_view()); @@ -624,7 +791,7 @@ evaluateFields (typename Traits::EvalData workset) " Input value: " + std::to_string(Albany::ADValue(mu)) + "\n"); } - if (beta_type== BETA_TYPE::REGULARIZED_COULOMB) { + if ( (beta_type== BETA_TYPE::REGULARIZED_COULOMB) || (beta_type == BETA_TYPE::DEBRIS_FRICTION) ) { if (logParameters) { if (bedRoughness_type == FIELD_TYPE::CONSTANT) { auto hostBedRoughness_view = Kokkos::create_mirror_view(bedRoughnessParam.get_view()); @@ -656,6 +823,21 @@ evaluateFields (typename Traits::EvalData workset) "\nError in LandIce::BasalFrictionCoefficient: \"Bed Roughness\" must be >= 0.\n"); } + if (beta_type == BETA_TYPE::DEBRIS_FRICTION) { + if (bulkFriction_type == FIELD_TYPE::CONSTANT) { + auto hostBulkFriction_view = Kokkos::create_mirror_view(bulkFrictionParam.get_view()); + Kokkos::deep_copy(hostBulkFriction_view, bulkFrictionParam.get_view()); + ScalarT hostBulkFriction = hostBulkFriction_view(0); + bulkFriction = Albany::convertScalar(hostBulkFriction); + } + if (basalDebris_type == FIELD_TYPE::CONSTANT) { + auto hostBasalDebris_view = Kokkos::create_mirror_view(basalDebrisParam.get_view()); + Kokkos::deep_copy(hostBasalDebris_view, basalDebrisParam.get_view()); + ScalarT hostBasalDebris = hostBasalDebris_view(0); + basalDebris = Albany::convertScalar(hostBasalDebris); + } + } + dim = nodal ? numNodes : numQPs; if (is_side_equation) { @@ -666,8 +848,6 @@ evaluateFields (typename Traits::EvalData workset) worksetSize = workset.numCells; } - TEUCHOS_TEST_FOR_EXCEPTION ((save_pressure_field && !std::is_constructible::value), std::logic_error, - "Error! BasalFrictionCoefficient: Trying to convert a FAD type (NVal) into a double (outN).\n"); Kokkos::parallel_for(BasalFrictionCoefficient_Policy(0, worksetSize), *this); } diff --git a/src/landIce/interfaceWithMPAS/Interface.cpp b/src/landIce/interfaceWithMPAS/Interface.cpp index 9863b7cb6e..85735b40f1 100644 --- a/src/landIce/interfaceWithMPAS/Interface.cpp +++ b/src/landIce/interfaceWithMPAS/Interface.cpp @@ -106,6 +106,8 @@ void velocity_solver_solve_fo(int nLayers, int globalVerticesStride, std::vector& effectivePressureData, const std::vector& muData, const std::vector& bedRoughnessData, + const std::vector& bulkFrictionData, + const std::vector& basalDebrisData, const std::vector& temperatureDataOnPrisms, std::vector& bodyForceMagnitudeOnBasalCell, std::vector& dissipationHeatOnPrisms, @@ -154,6 +156,8 @@ void velocity_solver_solve_fo(int nLayers, int globalVerticesStride, STKFieldType* dirichletField = meshStruct->metaData->get_field (stk::topology::NODE_RANK, "dirichlet_field"); STKFieldType* muField = meshStruct->metaData->get_field (stk::topology::NODE_RANK, "mu"); STKFieldType* bedRoughnessField = meshStruct->metaData->get_field (stk::topology::NODE_RANK, "bed_roughness"); + STKFieldType* bulkFrictionField = meshStruct->metaData->get_field (stk::topology::NODE_RANK, "bulk_friction_coefficient"); + STKFieldType* basalDebrisField = meshStruct->metaData->get_field (stk::topology::NODE_RANK, "basal_debris_factor"); STKFieldType* stiffeningFactorLogField = meshStruct->metaData->get_field (stk::topology::NODE_RANK, "stiffening_factor_log"); STKFieldType* effectivePressureField = meshStruct->metaData->get_field (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& indexToVertexID, - const std::vector& indexToTriangleID, double minBeta, - const std::vector& regulThk, - const std::vector& levelsNormalizedThickness, - const std::vector& elevationData, - const std::vector& thicknessData, - std::vector& betaData, - const std::vector& bedTopographyData, - const std::vector& smbData, - const std::vector& stiffeningFactorData, - const std::vector& effectivePressureData, - const std::vector& muData, - const std::vector& temperatureDataOnPrisms, - std::vector& bodyForceMagnitudeOnBasalCell, - std::vector& dissipationHeatOnPrisms, - std::vector& velocityOnVertices, - int& error, - const double& deltat) { - std::vector effectivePressureDataNonConst(effectivePressureData); - const std::vector 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("Zero Beta On Floating Ice", zeroBetaOnShelf); basalFrictionParams.set("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("Number Of Fields",0); diff --git a/src/landIce/problems/LandIce_Hydrology.hpp b/src/landIce/problems/LandIce_Hydrology.hpp index 43a28a5761..12718db7a5 100644 --- a/src/landIce/problems/LandIce_Hydrology.hpp +++ b/src/landIce/problems/LandIce_Hydrology.hpp @@ -883,6 +883,34 @@ Hydrology::constructEvaluators (PHX::FieldManager& fm0, ptr_lambda = Teuchos::rcp(new PHAL::SharedParameter(*p,dl)); fm0.template registerEvaluator(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("Parameter Name", param_name); + p->set>>("Accessors", this->getAccessors()->template at()); + p->set< Teuchos::RCP >("Parameter Library", paramLib); + p->set("Parameters List", ¶ms->sublist("Parameters")); + p->set("Default Nominal Value", params->sublist("LandIce Bulk Friction Coefficient").get(param_name,-1.0)); + + Teuchos::RCP> ptr_bulk_friction; + ptr_bulk_friction = Teuchos::rcp(new PHAL::SharedParameter(*p,dl)); + fm0.template registerEvaluator(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("Parameter Name", param_name); + p->set>>("Accessors", this->getAccessors()->template at()); + p->set< Teuchos::RCP >("Parameter Library", paramLib); + p->set("Parameters List", ¶ms->sublist("Parameters")); + p->set("Default Nominal Value", params->sublist("LandIce Basal Debris Factor").get(param_name,-1.0)); + + Teuchos::RCP> ptr_basal_debris; + ptr_basal_debris = Teuchos::rcp(new PHAL::SharedParameter(*p,dl)); + fm0.template registerEvaluator(ptr_basal_debris); + //--- Shared Parameter for basal friction coefficient: mu ---// p = Teuchos::rcp(new Teuchos::ParameterList("Basal Friction Coefficient: mu")); diff --git a/src/landIce/problems/LandIce_ParamEnum.hpp b/src/landIce/problems/LandIce_ParamEnum.hpp index 2cc09dc605..6279779aae 100644 --- a/src/landIce/problems/LandIce_ParamEnum.hpp +++ b/src/landIce/problems/LandIce_ParamEnum.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"; diff --git a/src/landIce/problems/LandIce_StokesFOBase.cpp b/src/landIce/problems/LandIce_StokesFOBase.cpp index 2b12674bbc..c10aead26b 100644 --- a/src/landIce/problems/LandIce_StokesFOBase.cpp +++ b/src/landIce/problems/LandIce_StokesFOBase.cpp @@ -100,6 +100,8 @@ StokesFOBase (const Teuchos::RCP& 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("Bed Roughness Name" ,"bed_roughness"); + bulk_friction_name = params->sublist("Variables Names").get("Bulk Friction Coefficient Name", "bulk_friction"); + basal_debris_name = params->sublist("Variables Names").get("Basal Debris Factor Name", "basal_debris"); flow_factor_name = params->sublist("Variables Names").get("Flow Factor Name" ,"flow_factor"); stiffening_factor_log_name = params->sublist("Variables Names").get("Stiffening Factor Log Name" ,"stiffening_factor_log"); damage_factor_name = params->sublist("Variables Names").get("Damage Factor Name" ,"damage_factor"); @@ -137,7 +139,7 @@ void StokesFOBase::buildProblem (Teuchos::ArrayRCPisParameter(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 tensorCubDegrees; tensorCubDegrees = params->get>(tensorCubDegName); if(cellType->getKey() == shards::Wedge<6>::key) @@ -149,9 +151,7 @@ void StokesFOBase::buildProblem (Teuchos::ArrayRCP(*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(*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("Mu Type")); if (mu_type!="CONSTANT") { auto mu_field_name = bfc.get("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); + } } } diff --git a/src/landIce/problems/LandIce_StokesFOBase.hpp b/src/landIce/problems/LandIce_StokesFOBase.hpp index c198e3e2dc..28b5417624 100644 --- a/src/landIce/problems/LandIce_StokesFOBase.hpp +++ b/src/landIce/problems/LandIce_StokesFOBase.hpp @@ -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(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("Parameter Name", param_name); + //p->set>>("Accessors", this->getAccessors()->template at()); + //p->set< Teuchos::RCP >("Parameter Library", paramLib); + //p->set("Parameters List", ¶ms->sublist("Parameters")); + //p->set("Default Nominal Value", params->sublist("LandIce Bulk Friction Coefficient").get(param_name,-1.0)); + + //Teuchos::RCP> ptr_bulk_friction; + + //--- Shared Parameter for bulk friction coefficient: bulk_friction ---// + param_name = "Bulk Friction Coefficient"; + + if(!PHAL::is_field_evaluated(fm0, param_name, dl->shared_param)) { + p = Teuchos::rcp(new Teuchos::ParameterList("Bulk Friction Coefficient: bulk_friction")); + p->set("Parameter Name", param_name); + p->set>>("Accessors", this->getAccessors()->template at()); + p->set< Teuchos::RCP >("Parameter Library", paramLib); + p->set("Parameters List", ¶ms->sublist("Parameters")); + p->set("Default Nominal Value", pl->sublist("Bulk Friction Coefficient").get(param_name,-1.0)); + + Teuchos::RCP> ptr_bulk_friction; + ptr_bulk_friction = Teuchos::rcp(new PHAL::SharedParameter(*p,dl)); + fm0.template registerEvaluator(ptr_bulk_friction); + } + + //--- Shared Parameter for basal debris factor: basal_debris ---// + param_name = "Basal Debris Factor"; + + if(!PHAL::is_field_evaluated(fm0, param_name, dl->shared_param)) { + p = Teuchos::rcp(new Teuchos::ParameterList("Basal Debris Factor: basal_debris")); + p->set("Parameter Name", param_name); + p->set>>("Accessors", this->getAccessors()->template at()); + p->set< Teuchos::RCP >("Parameter Library", paramLib); + p->set("Parameters List", ¶ms->sublist("Parameters")); + p->set("Default Nominal Value", pl->sublist("Basal Debris Factor").get(param_name,-1.0)); + + Teuchos::RCP> ptr_basal_debris; + ptr_basal_debris = Teuchos::rcp(new PHAL::SharedParameter(*p,dl)); + fm0.template registerEvaluator(ptr_basal_debris); + } //--- Shared Parameter for basal friction coefficient: muPowerLaw ---// param_name = ParamEnumName::Mu;