From ea2d3293a87dafdcb347937d328acc845fa014da Mon Sep 17 00:00:00 2001 From: Jeremy Brooks Date: Fri, 20 Dec 2024 12:04:58 -0800 Subject: [PATCH 1/4] 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. Note: N0 and u0 are hard-coded, but need to be definable via YAML file (as scalars, not fields) --- .../LandIce_BasalFrictionCoefficient.hpp | 18 +- .../LandIce_BasalFrictionCoefficient_Def.hpp | 217 ++++++++++++++++-- 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 ++++ 7 files changed, 334 insertions(+), 60 deletions(-) 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..8ae934b222 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,12 @@ 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; + ParamScalarT N0, u0; // declare regularization values for effective pressure and sliding velocity - 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 +490,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 +500,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 +548,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 +569,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 +645,67 @@ 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 +743,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 +794,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 +826,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 +851,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 17021bb32c..1605fc6486 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, @@ -327,6 +341,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); @@ -494,38 +509,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 @@ -641,6 +624,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); @@ -713,7 +698,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 2da02828db..33dd49557a 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"); @@ -138,7 +140,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) @@ -150,9 +152,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); } @@ -643,6 +643,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"); @@ -662,6 +663,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 6b0ca11ee2..b142251806 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; From 544532df3ff98077b91767fb5b3e18cb8e538f5d Mon Sep 17 00:00:00 2001 From: Mauro Perego Date: Wed, 17 Dec 2025 16:51:11 -0700 Subject: [PATCH 2/4] LandIce: fixes and improvements for the debris friction sliding law Add possiblity to set N0 and u0 via input file Fix parameter logic Added a test to exercise new sliding law --- .../LandIce_BasalFrictionCoefficient.hpp | 12 +- .../LandIce_BasalFrictionCoefficient_Def.hpp | 87 ++-- src/landIce/interfaceWithMPAS/Interface.cpp | 210 +++++--- src/landIce/problems/LandIce_StokesFOBase.cpp | 6 +- src/landIce/problems/LandIce_StokesFOBase.hpp | 26 +- tests/landIce/FO_GIS/CMakeLists.txt | 15 + ...t_fo_humboldt_depthInt_debrisFriction.yaml | 484 ++++++++++++++++++ 7 files changed, 710 insertions(+), 130 deletions(-) create mode 100644 tests/landIce/FO_GIS/input_fo_humboldt_depthInt_debrisFriction.yaml diff --git a/src/landIce/evaluators/LandIce_BasalFrictionCoefficient.hpp b/src/landIce/evaluators/LandIce_BasalFrictionCoefficient.hpp index de39b63658..a73ecf8cc9 100644 --- a/src/landIce/evaluators/LandIce_BasalFrictionCoefficient.hpp +++ b/src/landIce/evaluators/LandIce_BasalFrictionCoefficient.hpp @@ -49,10 +49,10 @@ class BasalFrictionCoefficient : public PHX::EvaluatorWithBaseImpl, 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 and DEBRIS_FRICTION only) - PHX::MDField powerParam; // [adim], Exponent (for POWER_LAW and REGULARIZED COULOMB and DEBRIS_FRICTION only) + PHX::MDField powerParam; // [adim], Exponent (for POWER_LAW and REGULARIZED COULOMB and DEBRIS_FRICTION) 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) - + PHX::MDField basalDebrisParam; // [kPa m^-1 s] basal debris factor (for DEBRIS_FRICTION only) + ScalarT printedMu; ScalarT printedBedRoughness; ScalarT printedQ; @@ -70,8 +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 bulkFrictionField; // [adim], bulk friction coefficient (for DEBRIS_FRICTION only) + PHX::MDField basalDebrisField; // [kPa m^-1 s] basal debris factor (for DEBRIS_FRICTION only) PHX::MDField N; // [kPa] PHX::MDField coordVec; // [km] PHX::MDField flowRate; // [Pa^{-n} s^{-1}] @@ -97,6 +97,8 @@ class BasalFrictionCoefficient : public PHX::EvaluatorWithBaseImpl, 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; diff --git a/src/landIce/evaluators/LandIce_BasalFrictionCoefficient_Def.hpp b/src/landIce/evaluators/LandIce_BasalFrictionCoefficient_Def.hpp index 8ae934b222..4ddc56d746 100644 --- a/src/landIce/evaluators/LandIce_BasalFrictionCoefficient_Def.hpp +++ b/src/landIce/evaluators/LandIce_BasalFrictionCoefficient_Def.hpp @@ -62,6 +62,8 @@ BasalFrictionCoefficient (const Teuchos::ParameterList& p, validPL.set("Flow Rate", 3.17e-24, "Constant Value for Flow Rate"); validPL.set("Effective Pressure Type", "Field", "Type of N: One, Field, Hydrostatic, Hydrostatic Computed At Nodes"); validPL.set("Effective Pressure", 1.0, "Effective Pressure [kPa]"); + validPL.set("Effective Pressure Regularization", 1000.0, "Effective Pressure Regularization [kPa]"); + validPL.set("Sliding Velocity Regularization", 500.0, "Sliding Velocity Regularization [m yr^{-1}]"); validPL.set("Minimum Fraction Overburden Pressure", 1.0, "Minimum Fraction Overburden Pressure"); validPL.set("Length Scale Factor", 1.0, "Length Scale Factor [km]"); validPL.set("Beta Field Name", "", "Name of the Field Mu"); @@ -174,7 +176,7 @@ BasalFrictionCoefficient (const Teuchos::ParameterList& p, 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" + << " beta = mu * N * |u|^{q-1} / [|u| + bedRoughness*A*N^n]^q + BulkFrictionCoefficient*N/(1+N/N0)/|u| + BasalDebrisFactor/(1+|u|/u0)\n\n" << " with N being the effective pressure, |u| the sliding velocity\n"; #endif @@ -217,6 +219,7 @@ 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); @@ -337,6 +340,11 @@ BasalFrictionCoefficient (const Teuchos::ParameterList& p, std::string bedRoughnessType = util::upper_case((beta_list.isParameter("Bed Roughness Type") ? beta_list.get("Bed Roughness Type") : "Field")); + N0 = beta_list.get("Effective Pressure Regularization"); + u0 = beta_list.get("Sliding Velocity Regularization"); + // N0 = 1.0e3; // Regularization value for effective pressure [kPa] + // u0 = 500.0; // Regularization value for sliding velocity [m yr^-1] + auto layout_bedRoughness_field = layout; if(bedRoughnessType == "CONSTANT") { bedRoughness_type = FIELD_TYPE::CONSTANT; @@ -473,12 +481,11 @@ postRegistrationSetup (typename Traits::SetupData d, template KOKKOS_INLINE_FUNCTION void BasalFrictionCoefficient:: -operator() (const BasalFrictionCoefficient_Tag& tag, const int& cell) const { +operator() (const BasalFrictionCoefficient_Tag&, const int& cell) const { ParamScalarT mu, bedRoughness, power, bulkFriction, basalDebris; - ParamScalarT N0, u0; // declare regularization values for effective pressure and sliding velocity - if ((beta_type == BETA_TYPE::POWER_LAW) || (beta_type==BETA_TYPE::REGULARIZED_COULOMB) || (beta_type==BETA_TYPE::DEBRIS_FRICTION) ) { + 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) @@ -548,9 +555,9 @@ operator() (const BasalFrictionCoefficient_Tag& 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) = static_cast(NVal); + outN(cell,ipt) = Albany::convertScalar(NVal); } - } else { + } else { MeshScalarT thickness(0), bed_topo(0); for (int node=0; node(NVal); + outN(cell,ipt) = Albany::convertScalar(NVal); } - } else { + } else { NVal = 0; for (int node=0; node @@ -759,7 +762,7 @@ evaluateFields (typename Traits::EvalData workset) auto hostPower_view = Kokkos::create_mirror_view(powerParam.get_view()); Kokkos::deep_copy(hostPower_view, powerParam.get_view()); ScalarT hostPower = hostPower_view(0); - power = std::exp(Albany::convertScalar(hostPower)); + power = Albany::convertScalar(hostPower); if (mu_type == FIELD_TYPE::CONSTANT) { auto hostMu_view = Kokkos::create_mirror_view(muParam.get_view()); Kokkos::deep_copy(hostMu_view, muParam.get_view()); @@ -851,6 +854,8 @@ 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 1605fc6486..b663ffea8c 100644 --- a/src/landIce/interfaceWithMPAS/Interface.cpp +++ b/src/landIce/interfaceWithMPAS/Interface.cpp @@ -94,14 +94,14 @@ typedef struct TET_ { 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& 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&, // smbData, const std::vector& stiffeningFactorData, std::vector& effectivePressureData, const std::vector& muData, @@ -152,7 +152,7 @@ void velocity_solver_solve_fo(int nLayers, int globalVerticesStride, STKFieldType* surfaceHeightField = meshStruct->metaData->get_field (stk::topology::NODE_RANK, "surface_height"); STKFieldType* thicknessField = meshStruct->metaData->get_field (stk::topology::NODE_RANK, "ice_thickness"); STKFieldType* bedTopographyField = meshStruct->metaData->get_field (stk::topology::NODE_RANK, "bed_topography"); - STKFieldType* smbField = meshStruct->metaData->get_field (stk::topology::NODE_RANK, "surface_mass_balance"); + //STKFieldType* smbField = meshStruct->metaData->get_field (stk::topology::NODE_RANK, "surface_mass_balance"); 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"); @@ -209,11 +209,12 @@ void velocity_solver_solve_fo(int nLayers, int globalVerticesStride, double* effectivePressure = stk::mesh::field_data(*effectivePressureField, node); effectivePressure[0] = effectivePressureData[ib]; } - - if(smbField != NULL) { - double* smb = stk::mesh::field_data(*smbField, node); - smb[0] = smbData[ib]; - } + + // Not needed + // if(smbField != NULL) { + // double* smb = stk::mesh::field_data(*smbField, node); + // smb[0] = smbData[ib]; + // } double* sol = stk::mesh::field_data(*solutionField, node); sol[0] = velocityOnVertices[j]; sol[1] = velocityOnVertices[j + numVertices3D]; @@ -341,7 +342,6 @@ 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); @@ -509,6 +509,40 @@ 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; + const std::vector bulkFrictionData; + const std::vector basalDebrisData; + 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, bulkFrictionData, basalDebrisData, + temperatureDataOnPrisms, bodyForceMagnitudeOnBasalCell,dissipationHeatOnPrisms, velocityOnVertices, error, deltat ); +} + + /*duality: * * mpas(F) | albany @@ -622,10 +656,20 @@ void velocity_solver_extrude_3d_grid(int nLayers, int globalTrianglesStride, basalFrictionParams.set("Power Exponent", basalFrictionParams.get("Power Exponent",1.0)); basalFrictionParams.set("Mu Type",basalFrictionParams.get("Mu Type","Field")); 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")); + + bool requireBedRoughness(false), requireBulkFriction(false), requireBasalDebris(false); + if((betaType == "REGULARIZED COULOMB") || (betaType == "DEBRIS FRICTION")) { + basalFrictionParams.set("Bed Roughness Type",basalFrictionParams.get("Bed Roughness Type","Field")); + requireBedRoughness = true; + if(betaType == "DEBRIS FRICTION") { + 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")); + requireBulkFriction = true; + requireBasalDebris = true; + } + } + basalFrictionParams.set("Zero Beta On Floating Ice", zeroBetaOnShelf); basalFrictionParams.set("Zero Effective Pressure On Floating Ice At Nodes", zeroEffectPressOnShelf); @@ -698,81 +742,113 @@ 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", "Extruded")); //set to Extruded is not defined + discretizationList->set("Method", discretizationList->get("Method", "STKExtruded")); //set to STKExtruded is not defined auto& rfi = discretizationList->sublist("Required Fields Info"); int fp = rfi.get("Number Of Fields",0); - discretizationList->sublist("Required Fields Info").set("Number Of Fields",fp+11); - Teuchos::ParameterList& field0 = discretizationList->sublist("Required Fields Info").sublist(util::strint("Field",0+fp)); - Teuchos::ParameterList& field1 = discretizationList->sublist("Required Fields Info").sublist(util::strint("Field",1+fp)); - Teuchos::ParameterList& field2 = discretizationList->sublist("Required Fields Info").sublist(util::strint("Field",2+fp)); - Teuchos::ParameterList& field3 = discretizationList->sublist("Required Fields Info").sublist(util::strint("Field",3+fp)); - Teuchos::ParameterList& field4 = discretizationList->sublist("Required Fields Info").sublist(util::strint("Field",4+fp)); - Teuchos::ParameterList& field5 = discretizationList->sublist("Required Fields Info").sublist(util::strint("Field",5+fp)); - Teuchos::ParameterList& field6 = discretizationList->sublist("Required Fields Info").sublist(util::strint("Field",6+fp)); - Teuchos::ParameterList& field7 = discretizationList->sublist("Required Fields Info").sublist(util::strint("Field",7+fp)); - Teuchos::ParameterList& field8 = discretizationList->sublist("Required Fields Info").sublist(util::strint("Field",8+fp)); - Teuchos::ParameterList& field9 = discretizationList->sublist("Required Fields Info").sublist(util::strint("Field",9+fp)); - Teuchos::ParameterList& field10 = discretizationList->sublist("Required Fields Info").sublist(util::strint("Field",10+fp)); //set temperature - field0.set("Field Name", "temperature"); - if(depthIntegratedModel) - field0.set("Field Type", "QuadPoint Scalar"); - else - field0.set("Field Type", "Elem Scalar"); - field0.set("Field Origin", "Mesh"); + { + Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); + field.set("Field Name", "temperature"); + if(depthIntegratedModel) + field.set("Field Type", "QuadPoint Scalar"); + else + field.set("Field Type", "Elem Scalar"); + field.set("Field Origin", "Mesh"); + } //set ice thickness - field1.set("Field Name", "ice_thickness"); - field1.set("Field Type", "Node Scalar"); - field1.set("Field Origin", "Mesh"); + { + Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); + field.set("Field Name", "ice_thickness"); + field.set("Field Type", "Node Scalar"); + field.set("Field Origin", "Mesh"); + } //set surface_height - field2.set("Field Name", "surface_height"); - field2.set("Field Type", "Node Scalar"); - field2.set("Field Origin", "Mesh"); + { + Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); + field.set("Field Name", "surface_height"); + field.set("Field Type", "Node Scalar"); + field.set("Field Origin", "Mesh"); + } //set bed topography - field3.set("Field Name", "bed_topography"); - field3.set("Field Type", "Node Scalar"); - field3.set("Field Origin", "Mesh"); + { + Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); + field.set("Field Name", "bed_topography"); + field.set("Field Type", "Node Scalar"); + field.set("Field Origin", "Mesh"); + } //set mu field - field4.set("Field Name", "mu"); - field4.set("Field Type", "Node Scalar"); - field4.set("Field Origin", "Mesh"); - - //set mu field - field5.set("Field Name", "bed_roughness"); - field5.set("Field Type", "Node Scalar"); - field5.set("Field Origin", "Mesh"); - - //set surface mass balance - field6.set("Field Name", "surface_mass_balance"); - field6.set("Field Type", "Node Scalar"); - field6.set("Field Origin", "Mesh"); + { + Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); + field.set("Field Name", "mu"); + field.set("Field Type", "Node Scalar"); + field.set("Field Origin", "Mesh"); + } //set dirichlet field - field7.set("Field Name", "dirichlet_field"); - field7.set("Field Type", "Node Vector"); - field7.set("Field Origin", "Mesh"); + { + Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); + field.set("Field Name", "dirichlet_field"); + field.set("Field Type", "Node Vector"); + field.set("Field Origin", "Mesh"); + } //set stiffening factor - field8.set("Field Name", "stiffening_factor_log"); - field8.set("Field Type", "Node Scalar"); - field8.set("Field Origin", "Mesh"); + { + Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); + field.set("Field Name", "stiffening_factor_log"); + field.set("Field Type", "Node Scalar"); + field.set("Field Origin", "Mesh"); + } //set effective pressure - field9.set("Field Name", "effective_pressure"); - field9.set("Field Type", "Node Scalar"); - field9.set("Field Origin", "Mesh"); - field9.set("Field Usage", "Input-Output"); + { + Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); + field.set("Field Name", "effective_pressure"); + field.set("Field Type", "Node Scalar"); + field.set("Field Origin", "Mesh"); + field.set("Field Usage", "Input-Output"); + } + + //set bed roughness + if(requireBedRoughness) { + Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); + field.set("Field Name", "bed_roughness"); + field.set("Field Type", "Node Scalar"); + field.set("Field Origin", "Mesh"); + } + + //set basal debris factor + if(requireBasalDebris){ + Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); + field.set("Field Name", "basal_debris"); + field.set("Field Type", "Node Scalar"); + field.set("Field Origin", "Mesh"); + } + + //set bulk friction + if(requireBulkFriction) { + Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); + field.set("Field Name", "bulk_friction"); + field.set("Field Type", "Node Scalar"); + field.set("Field Origin", "Mesh"); + } // Outputs - field10.set("Field Name", "body_force"); - field10.set("Field Type", "Elem Vector"); - field10.set("Field Usage", "Output"); + { + Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); + field.set("Field Name", "body_force"); + field.set("Field Type", "Elem Vector"); + field.set("Field Usage", "Output"); + } + + // set number of requested fields + rfi.set("Number Of Fields",fp); // Side set outputs auto& ss_pl =discretizationList->sublist("Side Set Discretizations"); diff --git a/src/landIce/problems/LandIce_StokesFOBase.cpp b/src/landIce/problems/LandIce_StokesFOBase.cpp index 33dd49557a..5e416c2b1a 100644 --- a/src/landIce/problems/LandIce_StokesFOBase.cpp +++ b/src/landIce/problems/LandIce_StokesFOBase.cpp @@ -152,7 +152,9 @@ 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); } @@ -643,7 +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 + basal debris factor + // debris friction also requires bulk friction coefficient and 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"); diff --git a/src/landIce/problems/LandIce_StokesFOBase.hpp b/src/landIce/problems/LandIce_StokesFOBase.hpp index b142251806..e112999323 100644 --- a/src/landIce/problems/LandIce_StokesFOBase.hpp +++ b/src/landIce/problems/LandIce_StokesFOBase.hpp @@ -1458,7 +1458,9 @@ 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"; @@ -1614,7 +1604,8 @@ void StokesFOBase::constructBasalBCEvaluators (PHX::FieldManagerset>>("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::ParameterList beta_list = pl->sublist("Basal Friction Coefficient"); + p->set("Default Nominal Value", beta_list.get(param_name,-1.0)); Teuchos::RCP> ptr_bulk_friction; ptr_bulk_friction = Teuchos::rcp(new PHAL::SharedParameter(*p,dl)); @@ -1630,7 +1621,8 @@ void StokesFOBase::constructBasalBCEvaluators (PHX::FieldManagerset>>("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::ParameterList beta_list = pl->sublist("Basal Friction Coefficient"); + p->set("Default Nominal Value", beta_list.get(param_name,-1.0)); Teuchos::RCP> ptr_basal_debris; ptr_basal_debris = Teuchos::rcp(new PHAL::SharedParameter(*p,dl)); @@ -1680,6 +1672,8 @@ void StokesFOBase::constructBasalBCEvaluators (PHX::FieldManagerset("Effective Pressure QP Variable Name", effective_pressure_side_name); p->set("Flow Rate Variable Name", side_fname(flow_factor_name, ssName),flow_factor_side_name); p->set("Bed Roughness Variable Name", bed_roughness_side_name); + p->set("Basal Debris Factor Variable Name", basal_debris_side_name); + p->set("Bulk Friction Coefficient Variable Name", bulk_friction_side_name); p->set("Side Set Name", ssName); p->set("Coordinate Vector Variable Name", side_fname(Albany::coord_vec_name, ssName)); p->set("Parameter List", &pl->sublist("Basal Friction Coefficient")); @@ -1796,6 +1790,8 @@ void StokesFOBase::constructSMBEvaluators (PHX::FieldManager 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 ------------------------- // diff --git a/tests/landIce/FO_GIS/CMakeLists.txt b/tests/landIce/FO_GIS/CMakeLists.txt index 39be54ad39..7a773d2776 100644 --- a/tests/landIce/FO_GIS/CMakeLists.txt +++ b/tests/landIce/FO_GIS/CMakeLists.txt @@ -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 @@ -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 @@ -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() diff --git a/tests/landIce/FO_GIS/input_fo_humboldt_depthInt_debrisFriction.yaml b/tests/landIce/FO_GIS/input_fo_humboldt_depthInt_debrisFriction.yaml new file mode 100644 index 0000000000..ecc7feec92 --- /dev/null +++ b/tests/landIce/FO_GIS/input_fo_humboldt_depthInt_debrisFriction.yaml @@ -0,0 +1,484 @@ +%YAML 1.1 +--- +ANONYMOUS: + Debug Output: + Write Solution to MatrixMarket: 0 + Problem: + Phalanx Graph Visualization Detail: 0 + Depth Integrated Model: true + Solution Method: Steady + Name: LandIce Stokes First Order 3D + LandIce Rigid Body Modes For Preconditioner: + Compute Constant Modes: true + Compute Rotation Modes: true + Compute Sensitivities: true + Basal Side Name: basalside + Surface Side Name: upperside + Flat Bed Approximation: false + Cubature Degrees (Horiz Vert): [4, 6] + Basal Cubature Degree: 4 + Surface Cubature Degree: 4 + LandIce Field Norm: + sliding_velocity_basalside: + Regularization Type: Given Value + Regularization Value: 1.0e-05 + Dirichlet BCs: + SDBC on NS extruded_boundary_node_set_3 for DOF U0 prescribe Field: velocity + SDBC on NS extruded_boundary_node_set_3 for DOF U1 prescribe Field: velocity + LandIce BCs: + Number: 2 + BC 0: + Type: Basal Friction + Side Set Name: basalside + Basal Friction Coefficient: + Type: Debris Friction + Mu Type: Exponent of Field + Mu Field Name: basal_friction_log + Power Exponent: 0.3333333333 + Effective Pressure Type: Hydrostatic Computed at Nodes + Bed Roughness Type: Field + Bulk Friction Coefficient Type: Field + Basal Debris Factor Type: Field + Effective Pressure Regularization: 1.0e3 + Sliding Velocity Regularization: 5.0e2 + #Zero Beta On Floating Ice: true + Flow Rate Type: Viscosity Flow Rate + BC 1: + Type: Lateral + Cubature Degree: 3 + Side Set Name: extruded_boundary_side_set_1 + Melange Force: 6.0e+07 #[N/m] + Melange Submerged Thickness Threshold: 0.1 #[km] + #Immersed Ratio: 0.893 + Parameters: + Number Of Parameters: 0 + LandIce Viscosity: + Use P0 Temperature: false + Extract Strain Rate Sq: true + Type: Glen's Law + Glen's Law Homotopy Parameter: 0.3 + Continuous Homotopy With Constant Initial Viscosity: true + Coefficient For Continuous Homotopy: 8.0 + Glen's Law A: 2.400304414e-24 # [Pa^-n s^-1 + Glen's Law n: 3.0 + Flow Rate Type: Temperature Based + LandIce Physical Parameters: + Conductivity of ice: 2.1 + Diffusivity temperate ice: 1.1e-08 + Heat capacity of ice: 2.009e+03 + Water Density: 1.028e+03 + Ice Density: 9.10e+02 + Gravity Acceleration: 9.81e+00 + Reference Temperature: 2.65e+02 + Clausius-Clapeyron Coefficient: 7.9e-08 + Ice Latent Heat Of Fusion: 3.34e+05 + Permeability factor: 1.0e-12 #1e-12 + Viscosity of water: 1.8e-03 + Omega exponent alpha: 2.0e+00 + Diffusivity homotopy exponent: -1.1e+00 + Body Force: + Type: FO INTERP SURF GRAD + Response Functions: + Number Of Responses: 1 + Response 0: + Type: Sum Of Responses + Number Of Responses: 2 + Response 1: + Name: Surface Mass Balance Mismatch + SMB Coefficient: 1.0 + Regularization Coefficient: 1.0e+00 + Scaling Coefficient: 5.8824e-07 + H Coefficient: 100.0e+00 + Response 0: + Scaling Coefficient: 5.8824e-07 + Asinh Scaling: 1.0e+01 + Name: Surface Velocity Mismatch + Regularization Coefficient: 1.0e+00 + Discretization: + Workset Size: -1 + Method: STKExtruded + Surface Height Field Name: surface_height + Number Of Time Derivatives: 0 + Exodus Output File Name: output/humboldt.exo + Use Serial Mesh: false + Columnwise Ordering: true + NumLayers: 1 + Thickness Field Name: ice_thickness + Extrude Basal Node Fields: [ice_thickness, surface_height, basal_friction_log, bed_topography, apparent_mass_balance, apparent_mass_balance_RMS] + Basal Node Fields Ranks: [1, 1, 1, 1, 1, 1] + Interpolate Basal Node Layered Fields: [temperature, velocity] + Basal Node Layered Fields Ranks: [1, 2] + Use Glimmer Spacing: true + Required Fields Info: + Number Of Fields: 10 + Field 0: + Field Name: ice_thickness + Field Type: Node Scalar + Field Origin: Mesh + Field 1: + Field Name: surface_height + Field Type: Node Scalar + Field Origin: Mesh + Field 2: + Field Name: basal_friction_log + Field Type: Node Scalar + Field Origin: Mesh + Field 3: + Field Name: temperature + Field Type: Node Scalar + Field Origin: Mesh + Field 4: + Field Name: bed_topography + Field Type: Node Scalar + Field Origin: Mesh + Field 5: + Field Name: apparent_mass_balance + Field Type: Node Scalar + Field Origin: Mesh + Field 6: + Field Name: apparent_mass_balance_RMS + Field Type: Node Scalar + Field Origin: Mesh + Field 7: + Field Name: bed_topography + Field Type: Node Scalar + Field Origin: Mesh + Field Usage: Output + Field 8: + Field Name: surface_height + Field Type: Node Scalar + Field Origin: Mesh + Field Usage: Output + Field 9: + Field Name: velocity + Field Type: Node Vector + Field Origin: Mesh + Side Set Discretizations: + Side Sets: [basalside, upperside] + basalside: + Workset Size: -1 + Method: Ioss + Number Of Time Derivatives: 0 + Restart Index: 1 + Use Serial Mesh: ${USE_SERIAL_MESH} + Exodus Input File Name: ../AsciiMeshes/Humboldt/humboldt_2d.exo + Exodus Output File Name: output/humboldt_basal.exo + Required Fields Info: + Number Of Fields: 15 + Field 0: + Field Name: ice_thickness + Field Type: Node Scalar + Field Origin: File + File Name: ../AsciiMeshes/Humboldt/thickness.ascii + Field 1: + Field Name: observed_ice_thickness + Field Type: Node Scalar + Field Origin: File + File Name: ../AsciiMeshes/Humboldt/thickness.ascii + Field 2: + Field Name: observed_ice_thickness_RMS + Field Type: Node Scalar + Field Origin: File + File Name: ../AsciiMeshes/Humboldt/thickness_uncertainty.ascii + Field 3: + Field Name: surface_height + Field Type: Node Scalar + Field Origin: File + File Name: ../AsciiMeshes/Humboldt/surface_height.ascii + Field 4: + Field Name: bed_topography + Field Type: Node Scalar + Field Origin: File + File Name: ../AsciiMeshes/Humboldt/bed_topography.ascii + Field 5: + Field Name: flux_divergence + Field Type: Elem Scalar + Field Usage: Output + Field 6: + Field Name: basal_friction_log + Field Type: Node Scalar + Field Origin: File + File Name: ../AsciiMeshes/Humboldt/basal_friction_log.ascii + Field 7: + Field Name: temperature + Field Type: Node Layered Scalar + Number Of Layers: 8 + Field Origin: File + File Name: ../AsciiMeshes/Humboldt/nodal_temperature.ascii + Field 8: + Field Name: apparent_mass_balance + Field Type: Node Scalar + Field Origin: File + File Name: ../AsciiMeshes/Humboldt/apparent_mass_balance.ascii + Field 9: + Field Name: apparent_mass_balance_RMS + Field Type: Node Scalar + Field Origin: File + File Name: ../AsciiMeshes/Humboldt/apparent_mass_balance_uncertainty.ascii + Field 10: + Field Name: velocity + Field Origin: File + Field Type: Node Layered Vector + Number Of Layers: 2 + Vector Dim: 2 + File Name: ../AsciiMeshes/Humboldt/extruded_surface_velocity.ascii + Field 11: + Field Name: bed_roughness + Field Type: Node Scalar + Field Origin: File + Field Value: 0.01 + Field 12: + Field Name: effective_pressure + Field Type: Node Scalar + Field Origin: Mesh + Field Usage: Output + Field 13: + Field Name: bulk_friction + Field Type: Node Scalar + Field Origin: File + Field Value: 0.0 + Field 14: + Field Name: basal_debris + Field Type: Node Scalar + Field Origin: File + Field Value: 0.1 + upperside: + Method: SideSetSTK + Number Of Time Derivatives: 0 + Exodus Output File Name: output/humboldt_upper.exo + Required Fields Info: + Number Of Fields: 2 + Field 0: + Field Name: observed_surface_velocity + Field Type: Node Vector + Vector Dim: 2 + Field Origin: File + File Name: ../AsciiMeshes/Humboldt/surface_velocity.ascii + Field 1: + Field Name: observed_surface_velocity_RMS + Field Type: Node Scalar + Field Origin: File + File Name: ../AsciiMeshes/Humboldt/surface_velocity_uncertainty.ascii + Piro: + Sensitivity Method: Adjoint + Write Only Converged Solution: false + NOX: + Thyra Group Options: + Function Scaling: None + Update Row Sum Scaling: Before Each Nonlinear Solve + Status Tests: + Test Type: Combo + Combo Type: OR + Number of Tests: 2 + Test 0: + Test Type: Combo + Combo Type: AND + Number of Tests: 1 + Test 0: + Test Type: NormF + Norm Type: Two Norm + Scale Type: Scaled + Tolerance: 1.0e-08 + Test 1: + Test Type: MaxIters + Maximum Iterations: 50 + Direction: + Method: Newton + Newton: + Forcing Term Method: Constant + Linear Solver: + Write Linear System: false + Tolerance: 1.0e-7 + Stratimikos Linear Solver: + NOX Stratimikos Options: {} + Stratimikos: + Linear Solver Type: Belos + Linear Solver Types: + Belos: + Solver Type: Block GMRES + Solver Types: + Block GMRES: + Output Frequency: 20 + Output Style: 1 + Verbosity: 33 + Maximum Iterations: 200 + Block Size: 1 + Num Blocks: 200 + Flexible Gmres: false + VerboseObject: + Verbosity Level: low + Preconditioner Type: MueLu + Preconditioner Types: + Ifpack2: + Overlap: 0 + Prec Type: Amesos2 + MueLu: + Matrix: + PDE equations: 2 + Factories: + myLineDetectionFact: + factory: LineDetectionFactory + 'linedetection: orientation': coordinates + mySemiCoarsenPFact1: + factory: SemiCoarsenPFactory + 'semicoarsen: coarsen rate': 14 + UncoupledAggregationFact2: + factory: UncoupledAggregationFactory + 'aggregation: ordering': graph + 'aggregation: max selected neighbors': 0 + 'aggregation: min agg size': 3 + 'aggregation: phase3 avoid singletons': true + MyCoarseMap2: + factory: CoarseMapFactory + Aggregates: UncoupledAggregationFact2 + myTentativePFact2: + 'tentative: calculate qr': true + factory: TentativePFactory + Aggregates: UncoupledAggregationFact2 + CoarseMap: MyCoarseMap2 + mySaPFact2: + 'sa: eigenvalue estimate num iterations': 10 + 'sa: damping factor': 1.3333333e+00 + factory: SaPFactory + P: myTentativePFact2 + myTransferCoordinatesFact: + factory: CoordinatesTransferFactory + CoarseMap: MyCoarseMap2 + Aggregates: UncoupledAggregationFact2 + myTogglePFact: + factory: TogglePFactory + 'semicoarsen: number of levels': 2 + TransferFactories: + P1: mySemiCoarsenPFact1 + P2: mySaPFact2 + Ptent1: mySemiCoarsenPFact1 + Ptent2: myTentativePFact2 + Nullspace1: mySemiCoarsenPFact1 + Nullspace2: myTentativePFact2 + myRestrictorFact: + factory: TransPFactory + P: myTogglePFact + myToggleTransferCoordinatesFact: + factory: ToggleCoordinatesTransferFactory + Chosen P: myTogglePFact + TransferFactories: + Coordinates1: mySemiCoarsenPFact1 + Coordinates2: myTransferCoordinatesFact + myRAPFact: + 'rap: fix zero diagonals': true + factory: RAPFactory + P: myTogglePFact + R: myRestrictorFact + TransferFactories: + For Coordinates: myToggleTransferCoordinatesFact + myRepartitionHeuristicFact: + factory: RepartitionHeuristicFactory + A: myRAPFact + 'repartition: min rows per proc': 1000 + 'repartition: max imbalance': 1.327e+00 + 'repartition: start level': 4 + myZoltanInterface: + factory: ZoltanInterface + A: myRAPFact + Coordinates: myToggleTransferCoordinatesFact + number of partitions: myRepartitionHeuristicFact + myRepartitionFact: + factory: RepartitionFactory + A: myRAPFact + Partition: myZoltanInterface + 'repartition: remap parts': true + number of partitions: myRepartitionHeuristicFact + myRebalanceProlongatorFact: + factory: RebalanceTransferFactory + type: Interpolation + P: myTogglePFact + Coordinates: myToggleTransferCoordinatesFact + Nullspace: myTogglePFact + myRebalanceRestrictionFact: + factory: RebalanceTransferFactory + type: Restriction + R: myRestrictorFact + myRebalanceAFact: + factory: RebalanceAcFactory + A: myRAPFact + TransferFactories: { } + mySmoother1: + factory: TrilinosSmoother + #type: LINESMOOTHING_BANDEDRELAXATION + type: LINESMOOTHING_TRIDIRELAXATION + 'smoother: pre or post': both + ParameterList: + 'relaxation: type': Gauss-Seidel + 'relaxation: sweeps': 1 + 'relaxation: damping factor': 1.0e+00 + 'block relaxation: decouple dofs': true + 'partitioner: PDE equations': 2 + mySmoother3: + factory: TrilinosSmoother + type: RELAXATION + 'smoother: pre or post': both + ParameterList: + 'relaxation: type': Gauss-Seidel + 'relaxation: sweeps': 1 + 'relaxation: damping factor': 1.0e+00 + mySmoother4: + factory: TrilinosSmoother + type: RELAXATION + 'smoother: pre or post': pre + ParameterList: + 'relaxation: type': Gauss-Seidel + 'relaxation: sweeps': 4 + 'relaxation: damping factor': 1.0e+00 + Hierarchy: + max levels: 7 + 'coarse: max size': 2000 + verbosity: Low + use kokkos refactor: false + Finest: + Smoother: mySmoother1 + CoarseSolver: mySmoother4 + P: myRebalanceProlongatorFact + Nullspace: myRebalanceProlongatorFact + CoarseNumZLayers: myLineDetectionFact + LineDetection_Layers: myLineDetectionFact + LineDetection_VertLineIds: myLineDetectionFact + A: myRebalanceAFact + Coordinates: myRebalanceProlongatorFact + Importer: myRepartitionFact + All: + startLevel: 1 + Smoother: mySmoother4 + CoarseSolver: mySmoother4 + P: myRebalanceProlongatorFact + Nullspace: myRebalanceProlongatorFact + CoarseNumZLayers: myLineDetectionFact + LineDetection_Layers: myLineDetectionFact + LineDetection_VertLineIds: myLineDetectionFact + A: myRebalanceAFact + Coordinates: myRebalanceProlongatorFact + Importer: myRepartitionFact + Rescue Bad Newton Solve: true + Line Search: + Full Step: + Full Step: 1.0 + Method: Backtrack + Nonlinear Solver: Line Search Based + Printing: + Output Precision: 3 + Output Processor: 0 + Output Information: + Error: true + Warning: true + Outer Iteration: true + Parameters: false + Details: false + Linear Solver Details: false + Stepper Iteration: true + Stepper Details: true + Stepper Parameters: true + Solver Options: + Status Test Check Type: Minimal + Regression For Response 0: + Test Value: 7.312004304193e+01 + Relative Tolerance: 1.0e-06 +... From 326b29157e137f794641e8e2bdcd3ee8b7bc47b2 Mon Sep 17 00:00:00 2001 From: Mauro Perego Date: Wed, 25 Feb 2026 18:08:41 -0700 Subject: [PATCH 3/4] LandIce: Fix naming of Basal Debris Factor field. --- src/landIce/interfaceWithMPAS/Interface.cpp | 4 ++-- src/landIce/problems/LandIce_StokesFOBase.cpp | 10 +++++----- .../input_fo_humboldt_depthInt_debrisFriction.yaml | 2 +- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/landIce/interfaceWithMPAS/Interface.cpp b/src/landIce/interfaceWithMPAS/Interface.cpp index b663ffea8c..c92cab165d 100644 --- a/src/landIce/interfaceWithMPAS/Interface.cpp +++ b/src/landIce/interfaceWithMPAS/Interface.cpp @@ -156,7 +156,7 @@ 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* bulkFrictionField = meshStruct->metaData->get_field (stk::topology::NODE_RANK, "bulk_friction"); 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"); @@ -826,7 +826,7 @@ void velocity_solver_extrude_3d_grid(int nLayers, int globalTrianglesStride, //set basal debris factor if(requireBasalDebris){ Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); - field.set("Field Name", "basal_debris"); + field.set("Field Name", "basal_debris_factor"); field.set("Field Type", "Node Scalar"); field.set("Field Origin", "Mesh"); } diff --git a/src/landIce/problems/LandIce_StokesFOBase.cpp b/src/landIce/problems/LandIce_StokesFOBase.cpp index 5e416c2b1a..cacf7ffd5e 100644 --- a/src/landIce/problems/LandIce_StokesFOBase.cpp +++ b/src/landIce/problems/LandIce_StokesFOBase.cpp @@ -101,7 +101,7 @@ StokesFOBase (const Teuchos::RCP& params_, 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"); + basal_debris_name = params->sublist("Variables Names").get("Basal Debris Factor Name", "basal_debris_factor"); 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"); @@ -658,7 +658,7 @@ 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; @@ -668,19 +668,19 @@ void StokesFOBase::setupEvaluatorRequests () 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"; + 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"; + 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"; + 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; diff --git a/tests/landIce/FO_GIS/input_fo_humboldt_depthInt_debrisFriction.yaml b/tests/landIce/FO_GIS/input_fo_humboldt_depthInt_debrisFriction.yaml index ecc7feec92..86449923cd 100644 --- a/tests/landIce/FO_GIS/input_fo_humboldt_depthInt_debrisFriction.yaml +++ b/tests/landIce/FO_GIS/input_fo_humboldt_depthInt_debrisFriction.yaml @@ -238,7 +238,7 @@ ANONYMOUS: Field Origin: File Field Value: 0.0 Field 14: - Field Name: basal_debris + Field Name: basal_debris_factor Field Type: Node Scalar Field Origin: File Field Value: 0.1 From 0def4aed27c8b8ac24e7ac9d23d1bb25201dae11 Mon Sep 17 00:00:00 2001 From: Mauro Perego Date: Tue, 3 Mar 2026 16:19:54 -0700 Subject: [PATCH 4/4] LandIce: Improve logic in the MPAS Interface Avoid importing fields as 3D fields when not needed --- src/landIce/interfaceWithMPAS/Interface.cpp | 255 ++++++++++++-------- 1 file changed, 154 insertions(+), 101 deletions(-) diff --git a/src/landIce/interfaceWithMPAS/Interface.cpp b/src/landIce/interfaceWithMPAS/Interface.cpp index c92cab165d..980c58531d 100644 --- a/src/landIce/interfaceWithMPAS/Interface.cpp +++ b/src/landIce/interfaceWithMPAS/Interface.cpp @@ -76,6 +76,8 @@ bool depthIntegratedModel(false); std::vector dirichletNodesIdsDepthInt; +std::set noupdate_fields; + Teuchos::RCP > solver; bool keptMesh =false; @@ -133,6 +135,8 @@ void velocity_solver_solve_fo(int nLayers, int globalVerticesStride, *MPAS_dt = deltat; + auto& probParamList = paramList->sublist("Problem"); + auto& layerThicknessRatio = meshStruct->layers_data.dz_ref; if(depthIntegratedModel) @@ -151,28 +155,15 @@ void velocity_solver_solve_fo(int nLayers, int globalVerticesStride, STKFieldType* surfaceHeightField = meshStruct->metaData->get_field (stk::topology::NODE_RANK, "surface_height"); STKFieldType* thicknessField = meshStruct->metaData->get_field (stk::topology::NODE_RANK, "ice_thickness"); - STKFieldType* bedTopographyField = meshStruct->metaData->get_field (stk::topology::NODE_RANK, "bed_topography"); - //STKFieldType* smbField = meshStruct->metaData->get_field (stk::topology::NODE_RANK, "surface_mass_balance"); 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"); - 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; - STKFieldType* effectivePressureBasalField; - - auto& probParamList = paramList->sublist("Problem"); - const auto& landiceBcList = probParamList.sublist("LandIce BCs"); - const auto& basalParams = landiceBcList.sublist("BC 0"); - const auto& basalFrictionParams = basalParams.sublist("Basal Friction Coefficient"); - const auto betaType = util::upper_case(basalFrictionParams.get("Type")); - auto ss_mesh = meshStruct->sideSetMeshStructs.at("basalside"); - auto ss_mesh_stk = Teuchos::rcp_dynamic_cast(ss_mesh,true); - betaField = ss_mesh_stk->metaData->get_field (stk::topology::NODE_RANK, "beta"); - effectivePressureBasalField = ss_mesh_stk->metaData->get_field (stk::topology::NODE_RANK, "effective_pressure"); + TEUCHOS_TEST_FOR_EXCEPTION(noupdate_fields.find("surface_height") != noupdate_fields.end(), std::runtime_error, + "Error! surface_height field needs to be set by MPAS.\n"); + TEUCHOS_TEST_FOR_EXCEPTION(noupdate_fields.find("ice_thickness") != noupdate_fields.end(), std::runtime_error, + "Error! ice_thickness field needs to be set by MPAS.\n"); + bool update_dirichlet = !velocityOnVertices.empty() && (dirichletField != nullptr) && (noupdate_fields.find("dirichlet_field") == noupdate_fields.end()); + bool update_stiffeningFactorLog = !stiffeningFactorData.empty() && (stiffeningFactorLogField != nullptr) && (noupdate_fields.find("stiffening_factor_log") == noupdate_fields.end()); for (int j = 0; j < numVertices3D; ++j) { int ib = (ordering == 0) * (j % lVertexColumnShift) @@ -195,60 +186,92 @@ void velocity_solver_solve_fo(int nLayers, int globalVerticesStride, double* coord = stk::mesh::field_data(*meshStruct->getCoordinatesField(), node); coord[2] = elevationData[ib] + (levelsNormalizedThickness[il]-1.0) * thicknessData[ib]; - double* thickness = stk::mesh::field_data(*thicknessField, node); thickness[0] = thicknessData[ib]; double* sHeight = stk::mesh::field_data(*surfaceHeightField, node); sHeight[0] = elevationData[ib]; - double* bedTopography = stk::mesh::field_data(*bedTopographyField, node); - bedTopography[0] = bedTopographyData[ib]; - double* stiffeningFactorLog = stk::mesh::field_data(*stiffeningFactorLogField, node); - stiffeningFactorLog[0] = std::log(stiffeningFactorData[ib]); - - if(!effectivePressureData.empty() && (effectivePressureField != nullptr)) { - double* effectivePressure = stk::mesh::field_data(*effectivePressureField, node); - effectivePressure[0] = effectivePressureData[ib]; + + if(update_stiffeningFactorLog) { + double* stiffeningFactorLog = stk::mesh::field_data(*stiffeningFactorLogField, node); + stiffeningFactorLog[0] = std::log(stiffeningFactorData[ib]); } - - // Not needed - // if(smbField != NULL) { - // double* smb = stk::mesh::field_data(*smbField, node); - // smb[0] = smbData[ib]; - // } + double* sol = stk::mesh::field_data(*solutionField, node); sol[0] = velocityOnVertices[j]; sol[1] = velocityOnVertices[j + numVertices3D]; if(neq==3) { sol[2] = thicknessData[ib]; } - double* dirichletVel = stk::mesh::field_data(*dirichletField, node); - dirichletVel[0]=velocityOnVertices[j]; //velocityOnVertices stores initial guess and dirichlet velocities. - dirichletVel[1]=velocityOnVertices[j + numVertices3D]; + if(update_dirichlet) { + double* dirichletVel = stk::mesh::field_data(*dirichletField, node); + dirichletVel[0]=velocityOnVertices[j]; //velocityOnVertices stores initial guess and dirichlet velocities. + dirichletVel[1]=velocityOnVertices[j + numVertices3D]; + } + } + + + // basal fields + const auto& landiceBcList = probParamList.sublist("LandIce BCs"); + const auto& basalParams = landiceBcList.sublist("BC 0"); + const auto& basalFrictionParams = basalParams.sublist("Basal Friction Coefficient"); + const auto betaType = util::upper_case(basalFrictionParams.get("Type")); + auto ss_mesh = meshStruct->sideSetMeshStructs.at("basalside"); + auto ss_mesh_stk = Teuchos::rcp_dynamic_cast(ss_mesh,true); + + STKFieldType* bedTopographyField = ss_mesh_stk->metaData->get_field (stk::topology::NODE_RANK, "bed_topography"); + STKFieldType* muField = ss_mesh_stk->metaData->get_field (stk::topology::NODE_RANK, "mu"); + STKFieldType* bedRoughnessField = ss_mesh_stk->metaData->get_field (stk::topology::NODE_RANK, "bed_roughness"); + STKFieldType* bulkFrictionField = ss_mesh_stk->metaData->get_field (stk::topology::NODE_RANK, "bulk_friction"); + STKFieldType* basalDebrisField = ss_mesh_stk->metaData->get_field (stk::topology::NODE_RANK, "basal_debris_factor"); + STKFieldType* betaField = ss_mesh_stk->metaData->get_field (stk::topology::NODE_RANK, "beta"); + STKFieldType* effectivePressureField = ss_mesh_stk->metaData->get_field (stk::topology::NODE_RANK, "effective_pressure"); + + bool update_bedTopography = !bedTopographyData.empty() && (bedTopographyField != nullptr) && (noupdate_fields.find("bed_topography") == noupdate_fields.end()); + bool update_mu = !muData.empty() && (muField != nullptr) && (noupdate_fields.find("mu") == noupdate_fields.end()); + bool update_bedRoughness = !bedRoughnessData.empty() && (bedRoughnessField != nullptr) && (noupdate_fields.find("bed_roughness") == noupdate_fields.end()); + bool update_bulkFriction = !bulkFrictionData.empty() && (bulkFrictionField != nullptr) && (noupdate_fields.find("bulk_friction") == noupdate_fields.end()); + bool update_basalDebris = !basalDebrisData.empty() && (basalDebrisField != nullptr) && (noupdate_fields.find("basal_debris_factor") == noupdate_fields.end()); + bool update_effectivePressure = !effectivePressureData.empty() && (effectivePressureField!=nullptr) && (noupdate_fields.find("effective_pressure") == noupdate_fields.end()); + for(int i = 0; i < (int) indexToVertexID.size(); ++i) { + stk::mesh::Entity node = ss_mesh_stk->bulkData->get_entity(stk::topology::NODE_RANK, indexToVertexID[i]); + int ib = vertexLayerShift * i; + + if (update_bedTopography) { + double* bedTopography = stk::mesh::field_data(*bedTopographyField, node); + bedTopography[0] = bedTopographyData[ib]; + } - if ((il == 0) && !muData.empty() && (muField != nullptr)) { + if (update_mu) { double* muVal = stk::mesh::field_data(*muField, node); muVal[0] = muData[ib]; } - if ((il == 0) && !bedRoughnessData.empty() && (bedRoughnessField != nullptr)) { + if (update_bedRoughness) { double* bedRoughnessVal = stk::mesh::field_data(*bedRoughnessField, node); bedRoughnessVal[0] = bedRoughnessData[ib]; } - if ((il == 0) && !bulkFrictionData.empty() && (bulkFrictionField != nullptr)) { + if (update_bulkFriction) { double* bulkFrictionVal = stk::mesh::field_data(*bulkFrictionField, node); bulkFrictionVal[0] = bulkFrictionData[ib]; } - if ((il == 0) && !basalDebrisData.empty() && (basalDebrisField != nullptr)) { + if (update_basalDebris) { double* basalDebrisVal = stk::mesh::field_data(*basalDebrisField, node); basalDebrisVal[0] = basalDebrisData[ib]; } + + if (update_effectivePressure) { + double* effectivePressureVal = stk::mesh::field_data(*effectivePressureField,node); + effectivePressureVal[0] = effectivePressureData[ib]; + } } //In the following we import the temperature field temperatureDataOnPrisms from MPAS, //which is stored as a constant in each Prism (wedge), into a temperature filed in Albany mesh. - if(depthIntegratedModel) { + STKFieldType* temperatureField = meshStruct->metaData->get_field(stk::topology::ELEMENT_RANK, "temperature"); + bool update_temperature = !temperatureDataOnPrisms.empty() && (temperatureField!=nullptr) && (noupdate_fields.find("temperature") == noupdate_fields.end()); + if(depthIntegratedModel && update_temperature) { //In this case the Albany mesh only has one layer, but the MPAS mesh can still have multiple layers bool usePOTemp = probParamList.sublist("LandIce Viscosity").get("Use P0 Temperature"); if(!usePOTemp) { @@ -279,11 +302,9 @@ void velocity_solver_solve_fo(int nLayers, int globalVerticesStride, } // Populate the temperature mesh field, defined at quad points - STKFieldType* temperature_field = meshStruct->metaData->get_field (stk::topology::ELEMENT_RANK, "temperature"); - for(int ib=0; ib < (int) indexToTriangleID.size(); ++ib ) { stk::mesh::Entity elem = meshStruct->bulkData->get_entity(stk::topology::ELEMENT_RANK, indexToTriangleID[ib]); - double* temperature = stk::mesh::field_data(*temperature_field, elem); + double* temperature = stk::mesh::field_data(*temperatureField, elem); for(int qp=0; qpmetaData->get_field(stk::topology::ELEMENT_RANK, "temperature"); for(int ib=0; ib < (int) indexToTriangleID.size(); ++ib ) { stk::mesh::Entity elem = meshStruct->bulkData->get_entity(stk::topology::ELEMENT_RANK, indexToTriangleID[ib]); - double* temperature = stk::mesh::field_data(*temperature_field, elem); + double* temperature = stk::mesh::field_data(*temperatureField, elem); for(int il=0; ilmetaData->get_field(stk::topology::ELEMENT_RANK, "temperature"); + } else if(update_temperature){ // Here we copy the temperature on Prisms from MPAS into a P0 field in the Albany mesh. for(int ib=0; ib < (int) indexToTriangleID.size(); ++ib ) { for(int il=0; ilbulkData->get_entity(stk::topology::ELEMENT_RANK, gId); - double* temperature = stk::mesh::field_data(*temperature_field, elem); + double* temperature = stk::mesh::field_data(*temperatureField, elem); temperature[0] = temperatureDataOnPrisms[lId]; } } @@ -413,10 +432,10 @@ void velocity_solver_solve_fo(int nLayers, int globalVerticesStride, } } - if (Teuchos::nonnull(ss_mesh_stk) && !effectivePressureData.empty() && (effectivePressureBasalField!=nullptr)) { + if (Teuchos::nonnull(ss_mesh_stk) && !effectivePressureData.empty() && (effectivePressureField!=nullptr)) { for(int ib = 0; ib < (int) indexToVertexID.size(); ++ib) { stk::mesh::Entity node = ss_mesh_stk->bulkData->get_entity(stk::topology::NODE_RANK, indexToVertexID[ib]); - const double* effectivePressureVal = stk::mesh::field_data(*effectivePressureBasalField,node); + const double* effectivePressureVal = stk::mesh::field_data(*effectivePressureField,node); effectivePressureData[vertexLayerShift * ib] = effectivePressureVal[0]; } } @@ -657,6 +676,8 @@ void velocity_solver_extrude_3d_grid(int nLayers, int globalTrianglesStride, basalFrictionParams.set("Mu Type",basalFrictionParams.get("Mu Type","Field")); basalFrictionParams.set("Mu Field Name",basalFrictionParams.get("Mu Field Name","mu")); basalFrictionParams.set("Effective Pressure Type",basalFrictionParams.get("Effective Pressure Type","Field")); + std::string effectivePressureType = util::upper_case(basalFrictionParams.get("Effective Pressure Type")); + bool isEPressureOutput = effectivePressureType.find("HYDROSTATIC")!=std::string::npos; bool requireBedRoughness(false), requireBulkFriction(false), requireBasalDebris(false); if((betaType == "REGULARIZED COULOMB") || (betaType == "DEBRIS FRICTION")) { @@ -746,8 +767,19 @@ void velocity_solver_extrude_3d_grid(int nLayers, int globalTrianglesStride, auto& rfi = discretizationList->sublist("Required Fields Info"); int fp = rfi.get("Number Of Fields",0); + noupdate_fields.clear(); + std::set existing_fields; + //first process fields added in the input yaml file + for (int i=0; i("Field Name"); + existing_fields.insert(field_name); + //these fields are either read from file or outputs, they should not be overwritten with MPAS fields. + if(rfi.sublist(util::strint("Field",i)).get("Field Origin", "Mesh") != "Mesh") + noupdate_fields.insert(field_name); + } //set temperature + if (existing_fields.find("temperature") == existing_fields.end()) { Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); field.set("Field Name", "temperature"); @@ -759,6 +791,7 @@ void velocity_solver_extrude_3d_grid(int nLayers, int globalTrianglesStride, } //set ice thickness + if (existing_fields.find("ice_thickness") == existing_fields.end()) { Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); field.set("Field Name", "ice_thickness"); @@ -767,6 +800,7 @@ void velocity_solver_extrude_3d_grid(int nLayers, int globalTrianglesStride, } //set surface_height + if (existing_fields.find("surface_height") == existing_fields.end()) { Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); field.set("Field Name", "surface_height"); @@ -774,102 +808,121 @@ void velocity_solver_extrude_3d_grid(int nLayers, int globalTrianglesStride, field.set("Field Origin", "Mesh"); } - //set bed topography + //set dirichlet field + if (existing_fields.find("dirichlet_field") == existing_fields.end()) { Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); - field.set("Field Name", "bed_topography"); - field.set("Field Type", "Node Scalar"); + field.set("Field Name", "dirichlet_field"); + field.set("Field Type", "Node Vector"); field.set("Field Origin", "Mesh"); } - //set mu field + //set stiffening factor + if (existing_fields.find("stiffening_factor_log") == existing_fields.end()) { Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); - field.set("Field Name", "mu"); + field.set("Field Name", "stiffening_factor_log"); field.set("Field Type", "Node Scalar"); field.set("Field Origin", "Mesh"); } - //set dirichlet field + // Outputs + if (existing_fields.find("body_force") == existing_fields.end()) { Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); - field.set("Field Name", "dirichlet_field"); - field.set("Field Type", "Node Vector"); - field.set("Field Origin", "Mesh"); + field.set("Field Name", "body_force"); + field.set("Field Type", "Elem Vector"); + field.set("Field Usage", "Output"); } - //set stiffening factor + // set number of requested fields + rfi.set("Number Of Fields",fp); + + // Side set outputs + auto& ss_pl =discretizationList->sublist("Side Set Discretizations"); + Teuchos::Array bsn (1,"basalside"); + ss_pl.set>("Side Sets", bsn); + auto& basal_pl = ss_pl.sublist("basalside"); + basal_pl.set("Method","SideSetSTK"); + auto& basal_req = basal_pl.sublist("Required Fields Info"); + + int bfp = basal_req.get("Number Of Fields",0); + existing_fields.clear(); + for (int i=0; i("Field Name"); + existing_fields.insert(field_name); + //these fields are either read from file or outputs, they should not be overwritten with MPAS fields. + if(basal_req.sublist(util::strint("Field",i)).get("Field Origin", "Mesh") != "Mesh") + noupdate_fields.insert(field_name); + } + + + //set bed topography + if (existing_fields.find("bed_topography") == existing_fields.end()) { - Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); - field.set("Field Name", "stiffening_factor_log"); + Teuchos::ParameterList& field = basal_req.sublist(util::strint("Field",bfp++)); + field.set("Field Name", "bed_topography"); field.set("Field Type", "Node Scalar"); field.set("Field Origin", "Mesh"); } - //set effective pressure + //set mu field + if (existing_fields.find("mu") == existing_fields.end()) { - Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); - field.set("Field Name", "effective_pressure"); + Teuchos::ParameterList& field = basal_req.sublist(util::strint("Field",bfp++)); + field.set("Field Name", "mu"); field.set("Field Type", "Node Scalar"); field.set("Field Origin", "Mesh"); - field.set("Field Usage", "Input-Output"); } //set bed roughness - if(requireBedRoughness) { - Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); + if(requireBedRoughness && (existing_fields.find("bed_roughness") == existing_fields.end())) { + Teuchos::ParameterList& field = basal_req.sublist(util::strint("Field",bfp++)); field.set("Field Name", "bed_roughness"); field.set("Field Type", "Node Scalar"); field.set("Field Origin", "Mesh"); } //set basal debris factor - if(requireBasalDebris){ - Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); + if(requireBasalDebris && (existing_fields.find("basal_debris_factor") == existing_fields.end())){ + Teuchos::ParameterList& field = basal_req.sublist(util::strint("Field",bfp++)); field.set("Field Name", "basal_debris_factor"); field.set("Field Type", "Node Scalar"); field.set("Field Origin", "Mesh"); } //set bulk friction - if(requireBulkFriction) { - Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); + if(requireBulkFriction && (existing_fields.find("bulk_friction") == existing_fields.end())) { + Teuchos::ParameterList& field = basal_req.sublist(util::strint("Field",bfp++)); field.set("Field Name", "bulk_friction"); field.set("Field Type", "Node Scalar"); field.set("Field Origin", "Mesh"); } - // Outputs + //set effective pressure + if (existing_fields.find("effective_pressure") == existing_fields.end()) { - Teuchos::ParameterList& field = rfi.sublist(util::strint("Field",fp++)); - field.set("Field Name", "body_force"); - field.set("Field Type", "Elem Vector"); - field.set("Field Usage", "Output"); + Teuchos::ParameterList& field = basal_req.sublist(util::strint("Field",bfp++)); + field.set("Field Name", "effective_pressure"); + field.set("Field Type", "Node Scalar"); + if(isEPressureOutput) { + field.set("Field Usage", "Output"); + noupdate_fields.insert("effective_pressure"); + } + else + field.set("Field Origin", "Mesh"); } - // set number of requested fields - rfi.set("Number Of Fields",fp); - - // Side set outputs - auto& ss_pl =discretizationList->sublist("Side Set Discretizations"); - Teuchos::Array bsn (1,"basalside"); - ss_pl.set>("Side Sets", bsn); - auto& basal_pl = ss_pl.sublist("basalside"); - basal_pl.set("Method","SideSetSTK"); - auto& basal_req = basal_pl.sublist("Required Fields Info"); - - int bfp = basal_req.get("Number Of Fields",0); - basal_req.set("Number Of Fields",bfp+2); - - auto& ss_field0 = basal_req.sublist(util::strint("Field",0+bfp)); - auto& ss_field1 = basal_req.sublist(util::strint("Field",1+bfp)); - ss_field0.set("Field Name", "beta"); - ss_field0.set("Field Type", "Node Scalar"); - ss_field0.set("Field Usage", "Output"); + //set beta + if (existing_fields.find("beta") == existing_fields.end()) + { + Teuchos::ParameterList& field = basal_req.sublist(util::strint("Field",bfp++)); + field.set("Field Name", "beta"); + field.set("Field Type", "Node Scalar"); + field.set("Field Usage", "Output"); + } - ss_field1.set("Field Name", "effective_pressure"); - ss_field1.set("Field Type", "Node Scalar"); - ss_field1.set("Field Usage", "Output"); + basal_req.set("Number Of Fields",bfp); // Register LandIce problems auto& pb_factories = Albany::FactoriesContainer::instance();