Skip to content

Commit

Permalink
LandIce: ViscosityFO: Allow evaluating temperature at quadrature points
Browse files Browse the repository at this point in the history 
Until now viscosity used P0 temperature fields.
Modified a couple of tests to exercise new capability.

Also, fix Enthalpy problem so that it uses the
pressure corrected temperature for computing the viscosity.
  • Loading branch information
@mperego
mperego committed May 12, 2023
1 parent 4aa0940 commit c139d04
Show file tree
Hide file tree
Showing 14 changed files with 202 additions and 149 deletions.
12 changes: 7 additions & 5 deletions src/landIce/evaluators/LandIce_PressureCorrectedTemperature.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,19 +43,21 @@ class PressureCorrectedTemperature: public PHX::EvaluatorWithBaseImpl<Traits>,
private:
// Input:
typedef typename Albany::StrongestScalarType<typename Albany::StrongestScalarType<TempST,MeshScalarT>::type, SurfHeightST>::type OutputScalarT;
PHX::MDField<const SurfHeightST,Cell> sHeight;
PHX::MDField<const TempST,Cell> temp;
PHX::MDField<const MeshScalarT,Cell,Dim> coord;

PHX::MDField<const SurfHeightST> sHeight;
PHX::MDField<const TempST> temp;
PHX::MDField<const MeshScalarT> coord;
// Output:
PHX::MDField<OutputScalarT,Cell> correctedTemp;
PHX::MDField<OutputScalarT> correctedTemp;

double beta; //[K Pa^{-1}]
double rho_i; //[kg m^{-3}]
double g; //[m s^{-2}]
double coeff; //[K km^{-1}]
double meltingT; //[K], 273.15

bool useP0Temp;
int numQPs;

PHAL::MDFieldMemoizer<Traits> memoizer;
};

Expand Down
36 changes: 29 additions & 7 deletions src/landIce/evaluators/LandIce_PressureCorrectedTemperature_Def.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,12 +21,28 @@ PressureCorrectedTemperature(const Teuchos::ParameterList& p, const Teuchos::RCP
coord (p.get<std::string> ("Coordinate Vector Variable Name"), dl->cell_gradient),
correctedTemp (p.get<std::string> ("Corrected Temperature Variable Name"), dl->cell_scalar2)
{
this->addDependentField(sHeight);
this->addDependentField(coord);
this->addDependentField(temp);
this->addEvaluatedField(correctedTemp);
useP0Temp = p.get<bool>("Use P0 Temperature");
if(useP0Temp) {
sHeight = decltype(sHeight)(p.get<std::string> ("Surface Height Variable Name"), dl->cell_scalar2);
temp = decltype(temp)(p.get<std::string> ("Temperature Variable Name"), dl->cell_scalar2);
coord = decltype(coord)(p.get<std::string> ("Coordinate Vector Variable Name"), dl->cell_gradient);
correctedTemp = decltype(correctedTemp)(p.get<std::string> ("Corrected Temperature Variable Name"), dl->cell_scalar2);
numQPs = 0;
} else {
sHeight = decltype(sHeight)(p.get<std::string> ("Surface Height Variable Name"), dl->qp_scalar);
temp = decltype(temp)(p.get<std::string> ("Temperature Variable Name"), dl->qp_scalar);
coord = decltype(coord)(p.get<std::string> ("Coordinate Vector Variable Name"), dl->qp_gradient);
correctedTemp = decltype(correctedTemp)(p.get<std::string> ("Corrected Temperature Variable Name"), dl->qp_scalar);
std::vector<PHX::DataLayout::size_type> dims;
dl->qp_scalar->dimensions(dims);
numQPs = dims[1];
}

this->setName("Pressure Corrected Temperature"+PHX::print<EvalT>());
this->addDependentField(sHeight);
this->addDependentField(coord);
this->addDependentField(temp);
this->addEvaluatedField(correctedTemp);
this->setName("Pressure Corrected Temperature"+PHX::print<EvalT>());

auto physicsList = p.get<Teuchos::ParameterList*>("LandIce Physical Parameters");
rho_i = physicsList->get<double>("Ice Density");
Expand Down Expand Up @@ -54,8 +70,14 @@ evaluateFields(typename Traits::EvalData workset)
{
if (memoizer.have_saved_data(workset,this->evaluatedFields())) return;

for (std::size_t cell = 0; cell < workset.numCells; ++cell)
correctedTemp(cell) = std::min(temp(cell) +coeff * (sHeight(cell) - coord(cell,2)), meltingT);
if(useP0Temp) {
for (std::size_t cell = 0; cell < workset.numCells; ++cell)
correctedTemp(cell) = std::min(temp(cell) + coeff * (sHeight(cell) - coord(cell,2)), meltingT);
} else {
for (std::size_t cell = 0; cell < workset.numCells; ++cell)
for (int qp = 0; qp < numQPs; ++qp)
correctedTemp(cell,qp) = std::min(temp(cell,qp) + coeff * (sHeight(cell,qp) - coord(cell,qp,2)), meltingT);
}
}

} // namespace LandIce
13 changes: 9 additions & 4 deletions src/landIce/evaluators/LandIce_ViscosityFO.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,9 +44,8 @@ class ViscosityFO : public PHX::EvaluatorWithBaseImpl<Traits>,
void evaluateFields(typename Traits::EvalData d);

private:
template<typename TemperatureT>
KOKKOS_INLINE_FUNCTION
TemperatureT flowRate(const TemperatureT& T) const;
TemprT flowRate(const TemprT& T) const;

const double pi, actenh, actenl, gascon, switchingT;
const double arrmlh, arrmll, scyr, k4scyr;
Expand All @@ -55,6 +54,7 @@ class ViscosityFO : public PHX::EvaluatorWithBaseImpl<Traits>,
bool extractStrainRateSq;
bool useStereographicMap;
bool useStiffeningFactor;
bool useP0Temp;
Teuchos::ParameterList* stereographicMapList;

// Coefficients for Glen's law
Expand All @@ -65,7 +65,7 @@ class ViscosityFO : public PHX::EvaluatorWithBaseImpl<Traits>,
PHX::MDField<const VelT,Cell,QuadPoint,VecDim,Dim> Ugrad; //[(k yr)^{-1}], k=1000
PHX::MDField<const VelT,Cell,QuadPoint,VecDim> U; //[m/yr]
PHX::MDField<const MeshScalarT,Cell,QuadPoint, Dim> coordVec; // [Km]
PHX::MDField<const TemprT,Cell> temperature; // [K]
PHX::MDField<const TemprT> temperature; // [K]
PHX::MDField<const RealType,Cell> flowFactorA; // [k^{-(n+1)} Pa^{-n} yr^{-1} ], k=1000. This is the coefficient A.
PHX::MDField<const ParamScalarT,Cell,QuadPoint> stiffeningFactor;
PHX::MDField<const ScalarT> homotopyParam;
Expand Down Expand Up @@ -95,12 +95,14 @@ class ViscosityFO : public PHX::EvaluatorWithBaseImpl<Traits>,
struct ViscosityFO_CONSTANT_Tag{};
struct ViscosityFO_GLENSLAW_UNIFORM_Tag{};
struct ViscosityFO_GLENSLAW_TEMPERATUREBASED_Tag{};
struct ViscosityFO_GLENSLAW_TEMPERATUREBASEDQP_Tag{};
struct ViscosityFO_GLENSLAW_FROMFILE_Tag{};

typedef Kokkos::RangePolicy<ExecutionSpace, ViscosityFO_EXPTRIG_Tag> ViscosityFO_EXPTRIG_Policy;
typedef Kokkos::RangePolicy<ExecutionSpace, ViscosityFO_CONSTANT_Tag> ViscosityFO_CONSTANT_Policy;
typedef Kokkos::RangePolicy<ExecutionSpace, ViscosityFO_GLENSLAW_UNIFORM_Tag> ViscosityFO_GLENSLAW_UNIFORM_Policy;
typedef Kokkos::RangePolicy<ExecutionSpace, ViscosityFO_GLENSLAW_TEMPERATUREBASED_Tag> ViscosityFO_GLENSLAW_TEMPERATUREBASED_Policy;
typedef Kokkos::RangePolicy<ExecutionSpace, ViscosityFO_GLENSLAW_TEMPERATUREBASEDQP_Tag> ViscosityFO_GLENSLAW_TEMPERATUREBASED_QP_Policy;
typedef Kokkos::RangePolicy<ExecutionSpace, ViscosityFO_GLENSLAW_FROMFILE_Tag> ViscosityFO_GLENSLAW_FROMFILE_Policy;

KOKKOS_INLINE_FUNCTION
Expand All @@ -115,11 +117,14 @@ class ViscosityFO : public PHX::EvaluatorWithBaseImpl<Traits>,
KOKKOS_INLINE_FUNCTION
void operator() (const ViscosityFO_GLENSLAW_TEMPERATUREBASED_Tag& tag, const int& i) const;

KOKKOS_INLINE_FUNCTION
void operator() (const ViscosityFO_GLENSLAW_TEMPERATUREBASEDQP_Tag& tag, const int& i) const;

KOKKOS_INLINE_FUNCTION
void operator() (const ViscosityFO_GLENSLAW_FROMFILE_Tag& tag, const int& i) const;

KOKKOS_INLINE_FUNCTION
void glenslaw (const ScalarT &flowFactorVec, const int& cell) const;
void glenslaw (const int& cell) const;

double R, x_0, y_0, R2;
};
Expand Down
Loading

0 comments on commit c139d04

Please sign in to comment.