From c060d2630ae8e4143dd77eb6cbcd45c7021cbdd6 Mon Sep 17 00:00:00 2001 From: Bo Gao <83304414+gaobhub@users.noreply.github.com> Date: Fri, 21 Jul 2023 16:39:46 -0400 Subject: [PATCH] Move calculation of soil resistance from landcover to sources; move WRM related parameters from evaluators to model parameters --- .../soil_resistance_model_partition.cc | 46 +++++ .../soil_resistance_model_partition.hh | 47 +++++ ...oil_resistance_sakagucki_zeng_evaluator.cc | 167 ++++++++++++++++++ ...oil_resistance_sakagucki_zeng_evaluator.hh | 81 +++++++++ ...resistance_sakagucki_zeng_evaluator_reg.hh | 25 +++ .../soil_resistance_sakagucki_zeng_model.hh | 129 ++++++++++++++ .../soil_resistance_sellers_evaluator.cc | 124 +++++++++++++ .../soil_resistance_sellers_evaluator.hh | 77 ++++++++ .../soil_resistance_sellers_evaluator_reg.hh | 25 +++ .../wrm/rel_perm_evaluator.cc | 4 +- .../wrm/rel_perm_frzBC_evaluator.cc | 4 +- .../wrm/rel_perm_sutraice_evaluator.cc | 4 +- .../wrm/suction_head_evaluator.cc | 4 +- .../wrm/wrm_evaluator.cc | 4 +- .../wrm/wrm_linear_relperm.cc | 4 +- .../wrm/wrm_permafrost_evaluator.cc | 4 +- .../wrm/wrm_zero_relperm.cc | 4 +- src/pks/surface_balance/CMakeLists.txt | 1 - .../land_cover/LandCover.cc | 7 - .../land_cover/LandCover.hh | 4 - .../evaporation_downregulation_evaluator.cc | 47 ++--- .../evaporation_downregulation_evaluator.hh | 4 +- .../evaporation_downregulation_model.cc | 95 ---------- .../evaporation_downregulation_model.hh | 32 ++-- .../land_cover/seb_physics_defs.hh | 14 +- .../land_cover/seb_physics_funcs.cc | 49 +---- .../land_cover/seb_physics_funcs.hh | 9 - .../seb_threecomponent_evaluator.cc | 38 +--- .../seb_threecomponent_evaluator.hh | 2 +- .../land_cover/seb_twocomponent_evaluator.cc | 35 +--- .../land_cover/seb_twocomponent_evaluator.hh | 2 +- 31 files changed, 783 insertions(+), 309 deletions(-) create mode 100644 src/pks/flow/constitutive_relations/sources/soil_resistance_model_partition.cc create mode 100644 src/pks/flow/constitutive_relations/sources/soil_resistance_model_partition.hh create mode 100644 src/pks/flow/constitutive_relations/sources/soil_resistance_sakagucki_zeng_evaluator.cc create mode 100644 src/pks/flow/constitutive_relations/sources/soil_resistance_sakagucki_zeng_evaluator.hh create mode 100644 src/pks/flow/constitutive_relations/sources/soil_resistance_sakagucki_zeng_evaluator_reg.hh create mode 100644 src/pks/flow/constitutive_relations/sources/soil_resistance_sakagucki_zeng_model.hh create mode 100644 src/pks/flow/constitutive_relations/sources/soil_resistance_sellers_evaluator.cc create mode 100644 src/pks/flow/constitutive_relations/sources/soil_resistance_sellers_evaluator.hh create mode 100644 src/pks/flow/constitutive_relations/sources/soil_resistance_sellers_evaluator_reg.hh delete mode 100644 src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_model.cc diff --git a/src/pks/flow/constitutive_relations/sources/soil_resistance_model_partition.cc b/src/pks/flow/constitutive_relations/sources/soil_resistance_model_partition.cc new file mode 100644 index 000000000..641c919aa --- /dev/null +++ b/src/pks/flow/constitutive_relations/sources/soil_resistance_model_partition.cc @@ -0,0 +1,46 @@ +/* + Copyright 2010-202x held jointly by participating institutions. + ATS is released under the three-clause BSD License. + The terms of use and "as is" disclaimer for this license are + provided in the top-level COPYRIGHT file. + + Authors: Ethan Coon (ecoon@lanl.gov) +*/ + +/* + A collection of WRMs along with a Mesh Partition. +*/ + +#include "dbc.hh" +#include "soil_resistance_model_partition.hh" + + +namespace Amanzi { +namespace Flow { + +// Non-member factory +Teuchos::RCP +createSoilResistanceModelPartition(Teuchos::ParameterList& plist) +{ + SoilResistanceSakaguckiZengModelList rs_list; + std::vector region_list; + + for (Teuchos::ParameterList::ConstIterator lcv = plist.begin(); lcv != plist.end(); ++lcv) { + std::string name = lcv->first; + if (plist.isSublist(name)) { + Teuchos::ParameterList sublist = plist.sublist(name); + region_list.push_back(sublist.get("region")); + rs_list.push_back(Teuchos::rcp(new SoilResistanceSakaguckiZengModel(sublist))); + } else { + AMANZI_ASSERT(0); + } + } + + Teuchos::RCP part = + Teuchos::rcp(new Functions::MeshPartition(AmanziMesh::CELL, region_list)); + + return Teuchos::rcp(new SoilResistanceModelPartition(part, rs_list)); +} + +} // namespace Flow +} // namespace Amanzi diff --git a/src/pks/flow/constitutive_relations/sources/soil_resistance_model_partition.hh b/src/pks/flow/constitutive_relations/sources/soil_resistance_model_partition.hh new file mode 100644 index 000000000..dab21f72a --- /dev/null +++ b/src/pks/flow/constitutive_relations/sources/soil_resistance_model_partition.hh @@ -0,0 +1,47 @@ +/* + Copyright 2010-202x held jointly by participating institutions. + ATS is released under the three-clause BSD License. + The terms of use and "as is" disclaimer for this license are + provided in the top-level COPYRIGHT file. + + Authors: Ethan Coon (ecoon@lanl.gov) +*/ + +//! A collection of WRMs along with a Mesh Partition. +/*! + +A WRM partition is a list of (region, WRM) pairs, where the regions partition +the mesh. + +.. _wrm-partition-typedinline-spec +.. admonition:: wrm-partition-typedinline-spec + + * `"region`" ``[string]`` Region on which the WRM is valid. + * `"WRM type`" ``[string]`` Name of the WRM type. + * `"_WRM_type_ parameters`" ``[_WRM_type_-spec]`` See below for the required + parameter spec for each type. + +*/ + +#ifndef AMANZI_FLOW_RELATIONS_SOIL_RESISTANCE_PARTITION_ +#define AMANZI_FLOW_RELATIONS_SOIL_RESISTANCE_PARTITION_ + +#include "soil_resistance_sakagucki_zeng_model.hh" +#include "MeshPartition.hh" + +namespace Amanzi { +namespace Flow { + +typedef std::vector> + SoilResistanceSakaguckiZengModelList; +typedef std::pair, SoilResistanceSakaguckiZengModelList> + SoilResistanceModelPartition; + +// Non-member factory +Teuchos::RCP +createSoilResistanceModelPartition(Teuchos::ParameterList& plist); + +} // namespace Flow +} // namespace Amanzi + +#endif diff --git a/src/pks/flow/constitutive_relations/sources/soil_resistance_sakagucki_zeng_evaluator.cc b/src/pks/flow/constitutive_relations/sources/soil_resistance_sakagucki_zeng_evaluator.cc new file mode 100644 index 000000000..0c5a73141 --- /dev/null +++ b/src/pks/flow/constitutive_relations/sources/soil_resistance_sakagucki_zeng_evaluator.cc @@ -0,0 +1,167 @@ +/* + Copyright 2010-202x held jointly by participating institutions. + ATS is released under the three-clause BSD License. + The terms of use and "as is" disclaimer for this license are + provided in the top-level COPYRIGHT file. + + Authors: Ethan Coon (ecoon@lanl.gov) +*/ + +/* + Evaluates the porosity, given a small compressibility of rock. + +*/ + +#include "Mesh_Algorithms.hh" +#include "soil_resistance_sakagucki_zeng_evaluator.hh" +#include "soil_resistance_sakagucki_zeng_model.hh" + +namespace Amanzi { +namespace Flow { + +SoilResistanceSakaguckiZengEvaluator::SoilResistanceSakaguckiZengEvaluator + (Teuchos::ParameterList& plist) : EvaluatorSecondaryMonotypeCV(plist) +{ + std::string domain_surf = Keys::getDomain(my_keys_.front().first); + Tag tag = my_keys_.front().second; + Key domain_ss = Keys::readDomainHint(plist, domain_surf, "surface", "subsurface"); + + sat_gas_key_ = Keys::readKey(plist_, domain_ss, "gas saturation", "saturation_gas"); + dependencies_.insert(KeyTag{ sat_gas_key_, tag }); + + poro_key_ = Keys::readKey(plist_, domain_ss, "porosity", "porosity"); + dependencies_.insert(KeyTag{ poro_key_, tag }); + + Teuchos::ParameterList& sublist = plist_.sublist("WRM parameters"); + sublist.remove("model type", false); + models_ = createSoilResistanceModelPartition(sublist); +} + +Teuchos::RCP +SoilResistanceSakaguckiZengEvaluator::Clone() const +{ + return Teuchos::rcp(new SoilResistanceSakaguckiZengEvaluator(*this)); +} + + +// Required methods from EvaluatorSecondaryMonotypeCV +void +SoilResistanceSakaguckiZengEvaluator::Evaluate_(const State& S, + const std::vector& result) +{ + Tag tag = my_keys_.front().second; + // Initialize the MeshPartition + if (!models_->first->initialized()) { + models_->first->Initialize(result[0]->Mesh()->parent(), -1); + models_->first->Verify(); + } + + Teuchos::RCP sat_gas = S.GetPtr(sat_gas_key_, tag); + Teuchos::RCP poro = S.GetPtr(poro_key_, tag); + Teuchos::RCP mesh = sat_gas->Mesh(); + Teuchos::RCP surf_mesh = result[0]->Mesh(); + + // evaluate the model + for (CompositeVector::name_iterator comp = result[0]->begin(); + comp != result[0]->end(); ++comp) { + AMANZI_ASSERT(*comp == "cell"); // partition on cell only + const Epetra_MultiVector& sat_gas_v = *(sat_gas->ViewComponent(*comp, false)); + const Epetra_MultiVector& poro_v = *(poro->ViewComponent(*comp, false)); + Epetra_MultiVector& result_v = *(result[0]->ViewComponent(*comp, false)); + + int count = result[0]->size(*comp); + for (int sc = 0; sc != count; ++sc) { + int f = surf_mesh->entity_get_parent(AmanziMesh::Entity_kind::CELL, sc); + int c = AmanziMesh::getFaceOnBoundaryInternalCell(*mesh, f); + result_v[0][sc] = models_->second[(*models_->first)[c]]-> + RsoilbySakagickiZeng(sat_gas_v[0][c], poro_v[0][c]); + } + } +} + + +void +SoilResistanceSakaguckiZengEvaluator::EvaluatePartialDerivative_( + const State& S, + const Key& wrt_key, + const Tag& wrt_tag, + const std::vector& result) +{ + // Initialize the MeshPartition + if (!models_->first->initialized()) { + models_->first->Initialize(result[0]->Mesh(), -1); + models_->first->Verify(); + } + + Tag tag = my_keys_.front().second; + Teuchos::RCP sat_gas = S.GetPtr(sat_gas_key_, tag); + Teuchos::RCP poro = S.GetPtr(poro_key_, tag); + Teuchos::RCP mesh = sat_gas->Mesh(); + Teuchos::RCP surf_mesh = result[0]->Mesh(); + + if (wrt_key == sat_gas_key_) { + // evaluate the model + for (CompositeVector::name_iterator comp = result[0]->begin(); + comp != result[0]->end(); ++comp) { + AMANZI_ASSERT(*comp == "cell"); // partition on cell only + const Epetra_MultiVector& sat_gas_v = *(sat_gas->ViewComponent(*comp, false)); + const Epetra_MultiVector& poro_v = *(poro->ViewComponent(*comp, false)); + Epetra_MultiVector& result_v = *(result[0]->ViewComponent(*comp, false)); + + int count = result[0]->size(*comp); + for (int sc = 0; sc != count; ++sc) { + int f = surf_mesh->entity_get_parent(AmanziMesh::Entity_kind::CELL, sc); + int c = AmanziMesh::getFaceOnBoundaryInternalCell(*mesh, f); + result_v[0][sc] = models_->second[(*models_->first)[c]]-> + DRsoilbySakagickiZengDSatGas(sat_gas_v[0][c], poro_v[0][c]); + } + } + + } else if (wrt_key == poro_key_) { + // evaluate the model + for (CompositeVector::name_iterator comp = result[0]->begin(); + comp != result[0]->end(); ++comp) { + AMANZI_ASSERT(*comp == "cell"); // partition on cell only + const Epetra_MultiVector& sat_gas_v = *(sat_gas->ViewComponent(*comp, false)); + const Epetra_MultiVector& poro_v = *(poro->ViewComponent(*comp, false)); + Epetra_MultiVector& result_v = *(result[0]->ViewComponent(*comp, false)); + + int count = result[0]->size(*comp); + for (int sc = 0; sc != count; ++sc) { + int f = surf_mesh->entity_get_parent(AmanziMesh::Entity_kind::CELL, sc); + int c = AmanziMesh::getFaceOnBoundaryInternalCell(*mesh, f); + result_v[0][sc] = models_->second[(*models_->first)[c]]-> + DRsoilbySakagickiZengDPorosity(sat_gas_v[0][c], poro_v[0][c]); + } + } + + } else { + AMANZI_ASSERT(0); + } +} + + +void +SoilResistanceSakaguckiZengEvaluator::EnsureCompatibility_ToDeps_(State& S) +{ + const auto& fac = S.Require(my_keys_.front().first, + my_keys_.front().second); + if (fac.Mesh() != Teuchos::null) { + CompositeVectorSpace dep_fac; + dep_fac.SetMesh(fac.Mesh()->parent()) + ->SetGhosted(true) + ->AddComponent("cell", AmanziMesh::CELL, 1); + + for (const auto& dep : dependencies_) { + if (Keys::getDomain(dep.first) == Keys::getDomain(my_keys_.front().first)) { + S.Require(dep.first, dep.second).Update(fac); + } else { + S.Require(dep.first, dep.second).Update(dep_fac); + } + } + } +} + + +} // namespace Flow +} // namespace Amanzi diff --git a/src/pks/flow/constitutive_relations/sources/soil_resistance_sakagucki_zeng_evaluator.hh b/src/pks/flow/constitutive_relations/sources/soil_resistance_sakagucki_zeng_evaluator.hh new file mode 100644 index 000000000..fc2eb3196 --- /dev/null +++ b/src/pks/flow/constitutive_relations/sources/soil_resistance_sakagucki_zeng_evaluator.hh @@ -0,0 +1,81 @@ +/* + Copyright 2010-202x held jointly by participating institutions. + ATS is released under the three-clause BSD License. + The terms of use and "as is" disclaimer for this license are + provided in the top-level COPYRIGHT file. + + Authors: Ethan Coon (ecoon@lanl.gov) +*/ + +/* + Evaluates the porosity, given a small compressibility of rock. + + Compressible grains are both physically realistic (based on bulk modulus) + and a simple way to provide a non-elliptic, diagonal term for helping + solvers to converge. + +*/ + +/*! + +Compressible grains are both physically realistic (based on bulk modulus) and a +simple way to provide a non-elliptic, diagonal term for helping solvers to +converge. + +`"evaluator type`" = `"compressible porosity`" + +.. _compressible-porosity-evaluator-spec +.. admonition:: compressible-porosity-evaluator-spec + + * `"compressible porosity model parameters`" ``[compressible-porosity-standard-model-spec-list]`` + + KEYS: + + - `"pressure`" **DOMAIN-pressure** + - `"base porosity`" **DOMAIN-base_porosity** + +*/ + + +#ifndef AMANZI_FLOWRELATIONS_SOIL_RESISTANCE_SAKAGUCKI_ZENG_EVALUATOR_HH_ +#define AMANZI_FLOWRELATIONS_SOIL_RESISTANCE_SAKAGUCKI_ZENG_EVALUATOR_HH_ + +#include "Factory.hh" +#include "EvaluatorSecondaryMonotype.hh" +#include "soil_resistance_model_partition.hh" + +namespace Amanzi { +namespace Flow { + +class SoilResistanceSakaguckiZengEvaluator : public EvaluatorSecondaryMonotypeCV { + public: + explicit SoilResistanceSakaguckiZengEvaluator(Teuchos::ParameterList& plist); + SoilResistanceSakaguckiZengEvaluator(const SoilResistanceSakaguckiZengEvaluator& other) = default; + Teuchos::RCP Clone() const override; + + Teuchos::RCP get_Models() { return models_; } + + protected: + // Required methods from EvaluatorSecondaryMonotypeCV + virtual void Evaluate_(const State& S, const std::vector& result) override; + virtual void EvaluatePartialDerivative_(const State& S, + const Key& wrt_key, + const Tag& wrt_tag, + const std::vector& result) override; + + virtual void EnsureCompatibility_ToDeps_(State& S) override; + + protected: + Key sat_gas_key_; + Key poro_key_; + + Teuchos::RCP models_; + + private: + static Utils::RegisteredFactory fac_; +}; + +} // namespace Flow +} // namespace Amanzi + +#endif diff --git a/src/pks/flow/constitutive_relations/sources/soil_resistance_sakagucki_zeng_evaluator_reg.hh b/src/pks/flow/constitutive_relations/sources/soil_resistance_sakagucki_zeng_evaluator_reg.hh new file mode 100644 index 000000000..73ff5dd6e --- /dev/null +++ b/src/pks/flow/constitutive_relations/sources/soil_resistance_sakagucki_zeng_evaluator_reg.hh @@ -0,0 +1,25 @@ +/* + Copyright 2010-202x held jointly by participating institutions. + ATS is released under the three-clause BSD License. + The terms of use and "as is" disclaimer for this license are + provided in the top-level COPYRIGHT file. + + Authors: Ethan Coon (ecoon@lanl.gov) +*/ + +/* + Evaluates the porosity, given a small compressibility of rock. + +*/ + +#include "soil_resistance_sakagucki_zeng_evaluator.hh" + +namespace Amanzi { +namespace Flow { + +// registry of method +Utils::RegisteredFactory + SoilResistanceSakaguckiZengEvaluator::fac_("sakagucki-zeng soil resistance"); + +} // namespace Flow +} // namespace Amanzi diff --git a/src/pks/flow/constitutive_relations/sources/soil_resistance_sakagucki_zeng_model.hh b/src/pks/flow/constitutive_relations/sources/soil_resistance_sakagucki_zeng_model.hh new file mode 100644 index 000000000..8346d393f --- /dev/null +++ b/src/pks/flow/constitutive_relations/sources/soil_resistance_sakagucki_zeng_model.hh @@ -0,0 +1,129 @@ +/* + Copyright 2010-202x held jointly by participating institutions. + ATS is released under the three-clause BSD License. + The terms of use and "as is" disclaimer for this license are + provided in the top-level COPYRIGHT file. + + Authors: Ethan Coon (ecoon@lanl.gov) +*/ + +//! A simple model for allowing porosity to vary with pressure. +/*! + +Based on a linear increase, i.e. + +.. math:: + + \phi = \phi_{base} + H(p - p_{atm}) * \alpha + +where :math:`H` is the heaviside function and :math:`\alpha` is the provided +compressibility. If the inflection point is set to zero, the above function is +exact. However, then the porosity function is not smooth (has discontinuous +derivatives), so the inflection point smooths this with a quadratic that +matches the value and derivative at the inflection point and is 0 with 0 slope +at atmospheric pressure. + +.. _compressible-porosity-standard-model-spec +.. admonition:: compressible-porosity-standard-model-spec + + * `"region`" ``[string]`` Region on which this is applied. + * `"pore compressibility [Pa^-1]`" ``[double]`` :math:`\alpha` as described above + * `"pore compressibility inflection point [Pa]`" ``[double]`` **1000** + + The inflection point above which the function is linear. + +Example: + +.. code-block:: xml + + + + + + + +*/ + +#ifndef AMANZI_FLOWRELATIONS_SOIL_RESISTANCE_SAKAGUCKI_ZENG_MODEL_HH_ +#define AMANZI_FLOWRELATIONS_SOIL_RESISTANCE_SAKAGUCKI_ZENG_MODEL_HH_ + +#include "Teuchos_ParameterList.hpp" +#include "dbc.hh" + +namespace Amanzi { +namespace Flow { + +class SoilResistanceSakaguckiZengModel { + public: + explicit SoilResistanceSakaguckiZengModel(Teuchos::ParameterList& plist) : plist_(plist) + { + InitializeFromPlist_(); + } + + + double RsoilbySakagickiZeng(double saturation_gas, double porosity) + { + double r_soil; + if (saturation_gas == 0.) { + r_soil = 0.; // ponded water + } else { + double vp_diffusion = 2.2e-5 * std::pow(porosity, 2) * std::pow(1 - sr_, 2 + 3 * b_); + double L_Rsoil = d_ * (std::exp(std::pow(saturation_gas, 5)) - 1) / (std::exp(1) - 1); + r_soil = L_Rsoil / vp_diffusion; + } + AMANZI_ASSERT(r_soil >= 0); + return r_soil; + } + + + double DRsoilbySakagickiZengDSatGas(double saturation_gas, double porosity) + { + double vp_diffusion = 2.2e-5 * std::pow(porosity, 2) * std::pow(1 - sr_, 2 + 3 * b_); + double coef = d_ / (std::exp(1) - 1) / vp_diffusion; + return coef * std::exp(std::pow(saturation_gas, 5)) * 5 * std::pow(saturation_gas, 4); + } + + + double DRsoilbySakagickiZengDPorosity(double saturation_gas, double porosity) + { + double L_Rsoil = d_ * (std::exp(std::pow(saturation_gas, 5)) - 1) / (std::exp(1) - 1); + double coef = L_Rsoil / 2.2e-5 * std::pow(1 - sr_, -2 - 3 * b_); + return -2 * coef * std::pow(porosity, -3); + } + + protected: + void InitializeFromPlist_() + { + d_ = plist_.get("dessicated zone thickness [m]", 0.1); + sr_ = plist_.get("residual saturation [-]"); + + if (plist_.get("WRM Type") == "van Genuchten") { + if (plist_.isParameter("van Genuchten m [-]")) { + double m = plist_.get("van Genuchten m [-]"); + double n = 1.0 / (1.0 - m); + double lambda = (n - 1) * (1 - std::pow(0.5, n / (n - 1))); + b_ = 1. / lambda; + } else { + double n = plist_.get("van Genuchten n [-]"); + double lambda = (n - 1) * (1 - std::pow(0.5, n / (n - 1))); + b_ = 1. / lambda; + } + } else if (plist_.get("WRM Type") == "Brooks-Corey") { + double lambda = plist_.get("Brooks Corey lambda [-]"); + b_ = 1. / lambda; + } else if (plist_.get("WRM Type") == "Clapp-Hornberger") { + b_ = plist_.get("Clapp-Hornberger b [-]"); + } + } + + protected: + Teuchos::ParameterList plist_; + double sr_; + double d_; + double b_; +}; + +} // namespace Flow +} // namespace Amanzi + +#endif diff --git a/src/pks/flow/constitutive_relations/sources/soil_resistance_sellers_evaluator.cc b/src/pks/flow/constitutive_relations/sources/soil_resistance_sellers_evaluator.cc new file mode 100644 index 000000000..fd005015b --- /dev/null +++ b/src/pks/flow/constitutive_relations/sources/soil_resistance_sellers_evaluator.cc @@ -0,0 +1,124 @@ +/* + Copyright 2010-202x held jointly by participating institutions. + ATS is released under the three-clause BSD License. + The terms of use and "as is" disclaimer for this license are + provided in the top-level COPYRIGHT file. + + Authors: Ethan Coon (ecoon@lanl.gov) +*/ + +/* + Evaluates the porosity, given a small compressibility of rock. + +*/ + +#include "Mesh_Algorithms.hh" +#include "soil_resistance_sellers_evaluator.hh" + +namespace Amanzi { +namespace Flow { + +SoilResistanceSellersEvaluator::SoilResistanceSellersEvaluator(Teuchos::ParameterList& plist) + : EvaluatorSecondaryMonotypeCV(plist) +{ + std::string domain_surf = Keys::getDomain(my_keys_.front().first); + Tag tag = my_keys_.front().second; + Key domain_ss = Keys::readDomainHint(plist, domain_surf, "surface", "subsurface"); + + // my dependencies + sat_key_ = Keys::readKey(plist_, domain_ss, "liquid saturation", "saturation_liquid"); + dependencies_.insert(KeyTag{ sat_key_, tag }); +} + + +Teuchos::RCP +SoilResistanceSellersEvaluator::Clone() const +{ + return Teuchos::rcp(new SoilResistanceSellersEvaluator(*this)); +} + + +// Required methods from EvaluatorSecondaryMonotypeCV +void +SoilResistanceSellersEvaluator::Evaluate_(const State& S, + const std::vector& result) +{ + Tag tag = my_keys_.front().second; + Teuchos::RCP sat = S.GetPtr(sat_key_, tag); + Teuchos::RCP mesh = sat->Mesh(); + Teuchos::RCP surf_mesh = result[0]->Mesh(); + + // evaluate the model + for (CompositeVector::name_iterator comp = result[0]->begin(); + comp != result[0]->end(); ++comp) { + AMANZI_ASSERT(*comp == "cell"); // partition on cell only + const Epetra_MultiVector& sat_v = *(sat->ViewComponent(*comp, false)); + Epetra_MultiVector& result_v = *(result[0]->ViewComponent(*comp, false)); + + int count = result[0]->size(*comp); + for (int sc = 0; sc != count; ++sc) { + int f = surf_mesh->entity_get_parent(AmanziMesh::Entity_kind::CELL, sc); + int c = AmanziMesh::getFaceOnBoundaryInternalCell(*mesh, f); + result_v[0][sc] = std::exp(8.206 - 4.255 * sat_v[0][c]); + } + } +} + + +void +SoilResistanceSellersEvaluator::EvaluatePartialDerivative_( + const State& S, + const Key& wrt_key, + const Tag& wrt_tag, + const std::vector& result) +{ + Tag tag = my_keys_.front().second; + Teuchos::RCP sat = S.GetPtr(sat_key_, tag); + Teuchos::RCP mesh = sat->Mesh(); + Teuchos::RCP surf_mesh = result[0]->Mesh(); + + if (wrt_key == sat_key_) { + // evaluate the model + for (CompositeVector::name_iterator comp = result[0]->begin(); + comp != result[0]->end(); ++comp) { + AMANZI_ASSERT(*comp == "cell"); // partition on cell only + const Epetra_MultiVector& sat_v = *(sat->ViewComponent(*comp, false)); + Epetra_MultiVector& result_v = *(result[0]->ViewComponent(*comp, false)); + + int count = result[0]->size(*comp); + for (int sc = 0; sc != count; ++sc) { + int f = surf_mesh->entity_get_parent(AmanziMesh::Entity_kind::CELL, sc); + int c = AmanziMesh::getFaceOnBoundaryInternalCell(*mesh, f); + result_v[0][sc] = -4.255 * std::exp(8.206 - 4.255 * sat_v[0][c]); + } + } + } else { + AMANZI_ASSERT(0); + } +} + + +void +SoilResistanceSellersEvaluator::EnsureCompatibility_ToDeps_(State& S) +{ + const auto& fac = S.Require(my_keys_.front().first, + my_keys_.front().second); + if (fac.Mesh() != Teuchos::null) { + CompositeVectorSpace dep_fac; + dep_fac.SetMesh(fac.Mesh()->parent()) + ->SetGhosted(true) + ->AddComponent("cell", AmanziMesh::CELL, 1); + + for (const auto& dep : dependencies_) { + if (Keys::getDomain(dep.first) == Keys::getDomain(my_keys_.front().first)) { + S.Require(dep.first, dep.second).Update(fac); + } else { + S.Require(dep.first, dep.second).Update(dep_fac); + } + } + } +} + + +} // namespace Flow +} // namespace Amanzi diff --git a/src/pks/flow/constitutive_relations/sources/soil_resistance_sellers_evaluator.hh b/src/pks/flow/constitutive_relations/sources/soil_resistance_sellers_evaluator.hh new file mode 100644 index 000000000..39dd5a247 --- /dev/null +++ b/src/pks/flow/constitutive_relations/sources/soil_resistance_sellers_evaluator.hh @@ -0,0 +1,77 @@ +/* + Copyright 2010-202x held jointly by participating institutions. + ATS is released under the three-clause BSD License. + The terms of use and "as is" disclaimer for this license are + provided in the top-level COPYRIGHT file. + + Authors: Ethan Coon (ecoon@lanl.gov) +*/ + +/* + Evaluates the porosity, given a small compressibility of rock. + + Compressible grains are both physically realistic (based on bulk modulus) + and a simple way to provide a non-elliptic, diagonal term for helping + solvers to converge. + +*/ + +/*! + +Compressible grains are both physically realistic (based on bulk modulus) and a +simple way to provide a non-elliptic, diagonal term for helping solvers to +converge. + +`"evaluator type`" = `"compressible porosity`" + +.. _compressible-porosity-evaluator-spec +.. admonition:: compressible-porosity-evaluator-spec + + * `"compressible porosity model parameters`" ``[compressible-porosity-standard-model-spec-list]`` + + KEYS: + + - `"pressure`" **DOMAIN-pressure** + - `"base porosity`" **DOMAIN-base_porosity** + +*/ + + +#ifndef AMANZI_FLOWRELATIONS_SOIL_RESISTANCE_SELLERS_EVALUATOR_HH_ +#define AMANZI_FLOWRELATIONS_SOIL_RESISTANCE_SELLERS_EVALUATOR_HH_ + +#include "Factory.hh" +#include "EvaluatorSecondaryMonotype.hh" + +namespace Amanzi { +namespace Flow { + +class SoilResistanceSellersModel; + +class SoilResistanceSellersEvaluator : public EvaluatorSecondaryMonotypeCV { + public: + explicit SoilResistanceSellersEvaluator(Teuchos::ParameterList& plist); + SoilResistanceSellersEvaluator(const SoilResistanceSellersEvaluator& other) = default; + Teuchos::RCP Clone() const override; + + protected: + // Required methods from EvaluatorSecondaryMonotypeCV + virtual void Evaluate_(const State& S, const std::vector& result) override; + virtual void EvaluatePartialDerivative_(const State& S, + const Key& wrt_key, + const Tag& wrt_tag, + const std::vector& result) override; + + virtual void EnsureCompatibility_ToDeps_(State& S) override; + + protected: + Key sat_key_; + + private: + static Utils::RegisteredFactory fac_; +}; + +} // namespace Flow +} // namespace Amanzi + +#endif diff --git a/src/pks/flow/constitutive_relations/sources/soil_resistance_sellers_evaluator_reg.hh b/src/pks/flow/constitutive_relations/sources/soil_resistance_sellers_evaluator_reg.hh new file mode 100644 index 000000000..a889f37de --- /dev/null +++ b/src/pks/flow/constitutive_relations/sources/soil_resistance_sellers_evaluator_reg.hh @@ -0,0 +1,25 @@ +/* + Copyright 2010-202x held jointly by participating institutions. + ATS is released under the three-clause BSD License. + The terms of use and "as is" disclaimer for this license are + provided in the top-level COPYRIGHT file. + + Authors: Ethan Coon (ecoon@lanl.gov) +*/ + +/* + Evaluates the porosity, given a small compressibility of rock. + +*/ + +#include "soil_resistance_sellers_evaluator.hh" + +namespace Amanzi { +namespace Flow { + +// registry of method +Utils::RegisteredFactory + SoilResistanceSellersEvaluator::fac_("sellers soil resistance"); + +} // namespace Flow +} // namespace Amanzi diff --git a/src/pks/flow/constitutive_relations/wrm/rel_perm_evaluator.cc b/src/pks/flow/constitutive_relations/wrm/rel_perm_evaluator.cc index db7d0e7a4..0e4af71a3 100644 --- a/src/pks/flow/constitutive_relations/wrm/rel_perm_evaluator.cc +++ b/src/pks/flow/constitutive_relations/wrm/rel_perm_evaluator.cc @@ -16,8 +16,8 @@ namespace Flow { RelPermEvaluator::RelPermEvaluator(Teuchos::ParameterList& plist) : EvaluatorSecondaryMonotypeCV(plist), min_val_(0.) { - AMANZI_ASSERT(plist_.isSublist("WRM parameters")); - Teuchos::ParameterList sublist = plist_.sublist("WRM parameters"); + Teuchos::ParameterList& sublist = plist_.sublist("WRM parameters"); + sublist.remove("model type", false); wrms_ = createWRMPartition(sublist); InitializeFromPlist_(); } diff --git a/src/pks/flow/constitutive_relations/wrm/rel_perm_frzBC_evaluator.cc b/src/pks/flow/constitutive_relations/wrm/rel_perm_frzBC_evaluator.cc index f4100b889..87ea98782 100644 --- a/src/pks/flow/constitutive_relations/wrm/rel_perm_frzBC_evaluator.cc +++ b/src/pks/flow/constitutive_relations/wrm/rel_perm_frzBC_evaluator.cc @@ -19,8 +19,8 @@ namespace Flow { RelPermFrzBCEvaluator::RelPermFrzBCEvaluator(Teuchos::ParameterList& plist) : EvaluatorSecondaryMonotypeCV(plist), min_val_(0.) { - AMANZI_ASSERT(plist_.isSublist("WRM parameters")); - Teuchos::ParameterList sublist = plist_.sublist("WRM parameters"); + Teuchos::ParameterList& sublist = plist_.sublist("freezing rel perm parameters"); + sublist.remove("model type", false); wrms_ = createWRMPartition(sublist); InitializeFromPlist_(); } diff --git a/src/pks/flow/constitutive_relations/wrm/rel_perm_sutraice_evaluator.cc b/src/pks/flow/constitutive_relations/wrm/rel_perm_sutraice_evaluator.cc index 2724bb439..90dae4ff2 100644 --- a/src/pks/flow/constitutive_relations/wrm/rel_perm_sutraice_evaluator.cc +++ b/src/pks/flow/constitutive_relations/wrm/rel_perm_sutraice_evaluator.cc @@ -18,8 +18,8 @@ namespace Flow { RelPermSutraIceEvaluator::RelPermSutraIceEvaluator(Teuchos::ParameterList& plist) : EvaluatorSecondaryMonotypeCV(plist), min_val_(0.) { - AMANZI_ASSERT(plist_.isSublist("WRM parameters")); - Teuchos::ParameterList sublist = plist_.sublist("WRM parameters"); + Teuchos::ParameterList& sublist = plist_.sublist("WRM parameters"); + sublist.remove("model type", false); wrms_ = createWRMPartition(sublist); InitializeFromPlist_(); } diff --git a/src/pks/flow/constitutive_relations/wrm/suction_head_evaluator.cc b/src/pks/flow/constitutive_relations/wrm/suction_head_evaluator.cc index 960570abc..4ac6230be 100644 --- a/src/pks/flow/constitutive_relations/wrm/suction_head_evaluator.cc +++ b/src/pks/flow/constitutive_relations/wrm/suction_head_evaluator.cc @@ -20,8 +20,8 @@ namespace Flow { SuctionHeadEvaluator::SuctionHeadEvaluator(Teuchos::ParameterList& plist) : EvaluatorSecondaryMonotypeCV(plist), min_val_(0.) { - AMANZI_ASSERT(plist_.isSublist("WRM parameters")); - Teuchos::ParameterList sublist = plist_.sublist("WRM parameters"); + Teuchos::ParameterList& sublist = plist_.sublist("WRM parameters"); + sublist.remove("model type", false); wrms_ = createWRMPartition(sublist); InitializeFromPlist_(); } diff --git a/src/pks/flow/constitutive_relations/wrm/wrm_evaluator.cc b/src/pks/flow/constitutive_relations/wrm/wrm_evaluator.cc index d402dc106..fdda233f1 100644 --- a/src/pks/flow/constitutive_relations/wrm/wrm_evaluator.cc +++ b/src/pks/flow/constitutive_relations/wrm/wrm_evaluator.cc @@ -22,8 +22,8 @@ namespace Flow { WRMEvaluator::WRMEvaluator(Teuchos::ParameterList& plist) : EvaluatorSecondaryMonotypeCV(plist), calc_other_sat_(true) { - AMANZI_ASSERT(plist_.isSublist("WRM parameters")); - Teuchos::ParameterList wrm_plist = plist_.sublist("WRM parameters"); + Teuchos::ParameterList& wrm_plist = plist_.sublist("WRM parameters"); + wrm_plist.remove("model type", false); wrms_ = createWRMPartition(wrm_plist); InitializeFromPlist_(); diff --git a/src/pks/flow/constitutive_relations/wrm/wrm_linear_relperm.cc b/src/pks/flow/constitutive_relations/wrm/wrm_linear_relperm.cc index acf9f8ba2..1932125bd 100644 --- a/src/pks/flow/constitutive_relations/wrm/wrm_linear_relperm.cc +++ b/src/pks/flow/constitutive_relations/wrm/wrm_linear_relperm.cc @@ -27,8 +27,8 @@ WRMLinearRelPerm::WRMLinearRelPerm(Teuchos::ParameterList& plist) : plist_(plist void WRMLinearRelPerm::InitializeFromPlist_() { - AMANZI_ASSERT(plist_.isSublist("WRM parameters")); - Teuchos::ParameterList sublist = plist_.sublist("WRM parameters"); + Teuchos::ParameterList& sublist = plist_.sublist("WRM parameters"); + sublist.remove("model type", false); WRMFactory fac; wrm_ = fac.createWRM(sublist); diff --git a/src/pks/flow/constitutive_relations/wrm/wrm_permafrost_evaluator.cc b/src/pks/flow/constitutive_relations/wrm/wrm_permafrost_evaluator.cc index 19ea2ceab..02f77a19e 100644 --- a/src/pks/flow/constitutive_relations/wrm/wrm_permafrost_evaluator.cc +++ b/src/pks/flow/constitutive_relations/wrm/wrm_permafrost_evaluator.cc @@ -26,8 +26,8 @@ WRMPermafrostEvaluator::WRMPermafrostEvaluator(Teuchos::ParameterList& plist) : EvaluatorSecondaryMonotypeCV(plist) { // get the WRMs - AMANZI_ASSERT(plist_.isSublist("WRM parameters")); - Teuchos::ParameterList wrm_plist = plist_.sublist("WRM parameters"); + Teuchos::ParameterList& wrm_plist = plist_.sublist("WRM parameters"); + wrm_plist.remove("model type", false); wrms_ = createWRMPartition(wrm_plist); // and the permafrost models diff --git a/src/pks/flow/constitutive_relations/wrm/wrm_zero_relperm.cc b/src/pks/flow/constitutive_relations/wrm/wrm_zero_relperm.cc index d6ea8b05b..9014543ea 100644 --- a/src/pks/flow/constitutive_relations/wrm/wrm_zero_relperm.cc +++ b/src/pks/flow/constitutive_relations/wrm/wrm_zero_relperm.cc @@ -27,8 +27,8 @@ WRMZeroRelPerm::WRMZeroRelPerm(Teuchos::ParameterList& plist) : plist_(plist) void WRMZeroRelPerm::InitializeFromPlist_() { - AMANZI_ASSERT(plist_.isSublist("WRM parameters")); - Teuchos::ParameterList sublist = plist_.sublist("WRM parameters"); + Teuchos::ParameterList& sublist = plist_.sublist("WRM parameters"); + sublist.remove("model type", false); WRMFactory fac; wrm_ = fac.createWRM(sublist); diff --git a/src/pks/surface_balance/CMakeLists.txt b/src/pks/surface_balance/CMakeLists.txt index 22f7afa5a..bc6c03868 100644 --- a/src/pks/surface_balance/CMakeLists.txt +++ b/src/pks/surface_balance/CMakeLists.txt @@ -14,7 +14,6 @@ set(ats_surface_balance_src_files constitutive_relations/land_cover/longwave_evaluator.cc constitutive_relations/land_cover/incident_shortwave_radiation_model.cc constitutive_relations/land_cover/incident_shortwave_radiation_evaluator.cc - constitutive_relations/land_cover/evaporation_downregulation_model.cc constitutive_relations/land_cover/evaporation_downregulation_evaluator.cc constitutive_relations/land_cover/LandCover.cc constitutive_relations/land_cover/drainage_evaluator.cc diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.cc b/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.cc index 5d99b335a..d1e4c6b39 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.cc +++ b/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.cc @@ -38,9 +38,6 @@ LandCover::LandCover(Teuchos::ParameterList& plist) beers_k_lw(plist.get("Beer's law extinction coefficient, longwave [-]", NAN)), snow_transition_depth(plist.get("snow transition depth [m]", NAN)), water_transition_depth(plist.get("water transition depth [m]", NAN)), - dessicated_zone_thickness(plist.get("dessicated zone thickness [m]", NAN)), - clapp_horn_b(plist.get("Clapp and Hornberger b [-]", NAN)), - rs_method(plist.get("soil resistance method", "sakagucki_zeng")), roughness_ground(plist.get("roughness length of bare ground [m]", NAN)), roughness_snow(plist.get("roughness length of snow [m]", NAN)), mannings_n(plist.get("Manning's n [?]", NAN)) @@ -133,10 +130,6 @@ checkValid(const std::string& region, const LandCover& lc, const std::string& pa throwInvalid(region, "snow transition depth [m]"); if (parname == "water_transition_depth" && std::isnan(lc.water_transition_depth)) throwInvalid(region, "water transition depth [m]"); - if (parname == "dessicated_zone_thickness" && std::isnan(lc.dessicated_zone_thickness)) - throwInvalid(region, "dessicated zone thickness [m]"); - if (parname == "clapp_horn_b" && std::isnan(lc.clapp_horn_b)) - throwInvalid(region, "Clapp and Hornberger b [-]"); if (parname == "roughness_ground" && std::isnan(lc.roughness_ground)) throwInvalid(region, "roughness length of bare ground [m]"); if (parname == "roughness_snow" && std::isnan(lc.roughness_snow)) diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.hh b/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.hh index d625230f6..3b16c2443 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.hh +++ b/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.hh @@ -177,10 +177,6 @@ struct LandCover { double water_transition_depth; // [m] // soil properties controlling evaporation - double dessicated_zone_thickness; // [m] Thickness over which vapor must diffuse - // when the soil is dry. - double clapp_horn_b; // [-] exponent of the WRM, Clapp & Hornberger eqn 1 - std::string rs_method; // soil resistance method {sakagucki_zeng, sellers, TBD} double roughness_ground; // [m] Fetch length for latent/sensible heat fluxes. double roughness_snow; // [m] Fetch length for latent/sensible heat fluxes. }; diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_evaluator.cc b/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_evaluator.cc index 7b7399f12..63d6b9156 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_evaluator.cc +++ b/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_evaluator.cc @@ -41,18 +41,12 @@ EvaporationDownregulationEvaluator::InitializeFromPlist_() domain_surf_ = Keys::getDomain(my_keys_.front().first); domain_sub_ = Keys::readDomainHint(plist_, domain_surf_, "surface", "domain"); - // sat gas, sat liquid, and porosity on subsurface - sat_gas_key_ = Keys::readKey(plist_, domain_sub_, "saturation gas", "saturation_gas"); - dependencies_.insert(KeyTag{ sat_gas_key_, tag }); - sat_liq_key_ = Keys::readKey(plist_, domain_sub_, "saturation liquid", "saturation_liquid"); - dependencies_.insert(KeyTag{ sat_liq_key_, tag }); - poro_key_ = Keys::readKey(plist_, domain_sub_, "porosity", "porosity"); - dependencies_.insert(KeyTag{ poro_key_, tag }); - // dependency: potential_evaporation on surface pot_evap_key_ = Keys::readKey(plist_, domain_surf_, "potential evaporation", "potential_evaporation"); dependencies_.insert(KeyTag{ pot_evap_key_, tag }); + rsoil_key_ = Keys::readKey(plist_, domain_surf_, "soil resistance", "soil_resistance"); + dependencies_.insert(KeyTag{ rsoil_key_, tag }); } @@ -61,12 +55,8 @@ EvaporationDownregulationEvaluator::Evaluate_(const State& S, const std::vector& result) { Tag tag = my_keys_.front().second; - const Epetra_MultiVector& sat_gas = - *S.Get(sat_gas_key_, tag).ViewComponent("cell", false); - const Epetra_MultiVector& sat_liq = - *S.Get(sat_liq_key_, tag).ViewComponent("cell", false); - const Epetra_MultiVector& poro = - *S.Get(poro_key_, tag).ViewComponent("cell", false); + const Epetra_MultiVector& rsoil = + *S.Get(rsoil_key_, tag).ViewComponent("cell", false); const Epetra_MultiVector& pot_evap = *S.Get(pot_evap_key_, tag).ViewComponent("cell", false); Epetra_MultiVector& surf_evap = *result[0]->ViewComponent("cell", false); @@ -81,7 +71,7 @@ EvaporationDownregulationEvaluator::Evaluate_(const State& S, for (AmanziMesh::Entity_ID sc : lc_ids) { auto c = sub_mesh.cells_of_column(sc)[0]; surf_evap[0][sc] = - region_model.second->Evaporation(sat_gas[0][c], poro[0][c], pot_evap[0][sc], sat_liq[0][c]); + region_model.second->Evaporation(pot_evap[0][sc], rsoil[0][sc]); } } } @@ -96,12 +86,8 @@ EvaporationDownregulationEvaluator::EvaluatePartialDerivative_( { Tag tag = my_keys_.front().second; if (wrt_key == pot_evap_key_) { - const Epetra_MultiVector& sat_gas = - *S.Get(sat_gas_key_, tag).ViewComponent("cell", false); - const Epetra_MultiVector& sat_liq = - *S.Get(sat_liq_key_, tag).ViewComponent("cell", false); - const Epetra_MultiVector& poro = - *S.Get(poro_key_, tag).ViewComponent("cell", false); + const Epetra_MultiVector& rsoil = + *S.Get(rsoil_key_, tag).ViewComponent("cell", false); const Epetra_MultiVector& pot_evap = *S.Get(pot_evap_key_, tag).ViewComponent("cell", false); Epetra_MultiVector& surf_evap = *result[0]->ViewComponent("cell", false); @@ -117,7 +103,7 @@ EvaporationDownregulationEvaluator::EvaluatePartialDerivative_( for (AmanziMesh::Entity_ID sc : lc_ids) { auto c = sub_mesh.cells_of_column(sc)[0]; surf_evap[0][sc] = region_model.second->DEvaporationDPotentialEvaporation( - sat_gas[0][c], poro[0][c], pot_evap[0][sc], sat_liq[0][c]); + pot_evap[0][sc], rsoil[0][sc]); } } } @@ -128,25 +114,14 @@ void EvaporationDownregulationEvaluator::EnsureCompatibility_ToDeps_(State& S) { if (!consistent_) { - land_cover_ = getLandCover(S.ICList().sublist("land cover types"), - { "dessicated_zone_thickness", "clapp_horn_b" }); - for (const auto& lc : land_cover_) { - models_[lc.first] = Teuchos::rcp(new EvaporationDownregulationModel(lc.second)); - } Tag tag = my_keys_.front().second; - S.Require(poro_key_, tag) - .SetMesh(S.GetMesh(domain_sub_)) - ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1); - S.Require(sat_gas_key_, tag) - .SetMesh(S.GetMesh(domain_sub_)) - ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1); - S.Require(sat_liq_key_, tag) - .SetMesh(S.GetMesh(domain_sub_)) - ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1); S.Require(pot_evap_key_, tag) .SetMesh(S.GetMesh(domain_surf_)) ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1); + S.Require(rsoil_key_, tag) + .SetMesh(S.GetMesh(domain_surf_)) + ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1); consistent_ = true; } diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_evaluator.hh b/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_evaluator.hh index 457b690f0..5ff6174d3 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_evaluator.hh +++ b/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_evaluator.hh @@ -60,9 +60,7 @@ class EvaporationDownregulationEvaluator : public EvaluatorSecondaryMonotypeCV { protected: void InitializeFromPlist_(); - Key sat_gas_key_; - Key sat_liq_key_; - Key poro_key_; + Key rsoil_key_; Key pot_evap_key_; Key domain_surf_; diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_model.cc b/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_model.cc deleted file mode 100644 index 8a3995ca9..000000000 --- a/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_model.cc +++ /dev/null @@ -1,95 +0,0 @@ -/* - Copyright 2010-202x held jointly by participating institutions. - ATS is released under the three-clause BSD License. - The terms of use and "as is" disclaimer for this license are - provided in the top-level COPYRIGHT file. - - Authors: Ethan Coon (ecoon@lanl.gov) -*/ - -//! Downregulates bare soil evaporation through a dessicated zone. -#include "Teuchos_ParameterList.hpp" -#include "dbc.hh" -#include "evaporation_downregulation_model.hh" -#include "seb_physics_funcs.hh" - -namespace Amanzi { -namespace SurfaceBalance { -namespace Relations { - -// Constructor from ParameterList -EvaporationDownregulationModel::EvaporationDownregulationModel(Teuchos::ParameterList& plist) -{ - dess_dz_ = plist.get("dessicated zone thickness [m]", 0.1); - Clapp_Horn_b_ = plist.get("Clapp and Hornberger b of surface soil [-]", 1.0); - rs_method_ = plist.get("Soil resistance method", "sakagucki_zeng"); -} - -// Constructor from LandCover -EvaporationDownregulationModel::EvaporationDownregulationModel(const LandCover& lc) -{ - dess_dz_ = lc.dessicated_zone_thickness; - Clapp_Horn_b_ = lc.clapp_horn_b; - rs_method_ = lc.rs_method; -} - -// main method -double -EvaporationDownregulationModel::Evaporation(double sg, - double poro, - double pot_evap, - double sl) const -{ - double rsoil; - if (rs_method_ == "sellers") { - rsoil = Relations::EvaporativeResistanceCoefSellers(sl); - } else { - rsoil = Relations::EvaporativeResistanceCoefSakaguckiZeng(sg, poro, dess_dz_, Clapp_Horn_b_); - } - return pot_evap / (1. + rsoil); -} - -double -EvaporationDownregulationModel::DEvaporationDSaturationGas(double sg, - double poro, - double pot_evap, - double sl) const -{ - return 0.; -} - -double -EvaporationDownregulationModel::DEvaporationDPorosity(double sg, - double poro, - double pot_evap, - double sl) const -{ - return 0.; -} - -double -EvaporationDownregulationModel::DEvaporationDSaturationLiquid(double sg, - double poro, - double pot_evap, - double sl) const -{ - return 0.; -} - -double -EvaporationDownregulationModel::DEvaporationDPotentialEvaporation(double sg, - double poro, - double pot_evap, - double sl) const -{ - if (rs_method_ == "sellers") { - return 1. / (1. + Relations::EvaporativeResistanceCoefSellers(sl)); - } else { - return 1. / (1. + Relations::EvaporativeResistanceCoefSakaguckiZeng( - sg, poro, dess_dz_, Clapp_Horn_b_)); - } -} - -} // namespace Relations -} // namespace SurfaceBalance -} // namespace Amanzi diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_model.hh b/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_model.hh index f241a4095..fe080a860 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_model.hh +++ b/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_model.hh @@ -32,9 +32,7 @@ These may be provided via parameter list or LandCover type. */ -#pragma once - -#include "LandCover.hh" +#include "Teuchos_ParameterList.hpp" namespace Amanzi { namespace SurfaceBalance { @@ -42,25 +40,23 @@ namespace Relations { class EvaporationDownregulationModel { public: - explicit EvaporationDownregulationModel(double dessicated_zone_thickness, double clapp_horn_b); - EvaporationDownregulationModel(Teuchos::ParameterList& plist); - EvaporationDownregulationModel(const LandCover& lc); + explicit EvaporationDownregulationModel(Teuchos::ParameterList& plist) {} - double Evaporation(double sg, double poro, double pot_evap, double sl) const; + double Evaporation(double pot_evap, double rsoil) + { + return pot_evap / (1. + rsoil); + } - double DEvaporationDSaturationGas(double sg, double poro, double pot_evap, double sl) const; - double DEvaporationDPorosity(double sg, double poro, double pot_evap, double sl) const; - double DEvaporationDSaturationLiquid(double sg, double poro, double pot_evap, double sl) const; - double - DEvaporationDPotentialEvaporation(double sg, double poro, double pot_evap, double sl) const; + double DEvaporationDPotentialEvaporation(double pot_evap, double rsoil) + { + return 1. / (1. + rsoil); + } - protected: - void InitializeFromPlist_(Teuchos::ParameterList& plist); + double DEvaporationDRsoil(double pot_evap, double rsoil) + { + return -pot_evap * std::pow(1. + rsoil, -2); + } - protected: - double dess_dz_; - double Clapp_Horn_b_; - std::string rs_method_; }; } // namespace Relations diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_defs.hh b/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_defs.hh index c7279b120..ea92d3846 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_defs.hh +++ b/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_defs.hh @@ -112,15 +112,10 @@ struct GroundProperties { double temp; // temperature [K] double pressure; // [Pa] double ponded_depth; // [m] - double porosity; // [-] double density_w; // density [kg/m^3] - double dz; // [m] - double clapp_horn_b; // [-] - std::string rs_method; // options: {sakagucki_zeng, sellers} double albedo; // [-] double emissivity; // [-] - double saturation_gas; // [-] - double saturation_liq; // [-] + double rsoil; // [s/m] soil resistance at top cell double roughness; // [m] surface roughness of a bare domain double snow_death_rate; // [kg/m^2/s] snow that must die this timestep, make it melt! double unfrozen_fraction; // [-] fraction of ground water that is unfrozen @@ -130,16 +125,11 @@ struct GroundProperties { GroundProperties() : temp(NaN), pressure(NaN), - porosity(NaN), density_w(NaN), - dz(NaN), albedo(NaN), emissivity(NaN), - saturation_gas(NaN), - saturation_liq(NaN), + rsoil(NaN), roughness(NaN), - clapp_horn_b(3), - rs_method("sakagucki_zeng"), snow_death_rate(0.), unfrozen_fraction(0.), water_transition_depth(0.01) diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_funcs.cc b/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_funcs.cc index c10d922c6..70942cab3 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_funcs.cc +++ b/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_funcs.cc @@ -175,61 +175,16 @@ VaporPressureGround(const GroundProperties& surf, const ModelParams& params) double EvaporativeResistanceGround(const GroundProperties& surf, const MetData& met, - const ModelParams& params, double vapor_pressure_ground) { // calculate evaporation prefactors if (met.vp_air > vapor_pressure_ground) { // condensation return 0.; } else { - if (surf.rs_method == "sakagucki_zeng") { - return EvaporativeResistanceCoefSakaguckiZeng( - surf.saturation_gas, surf.porosity, surf.dz, surf.clapp_horn_b); - } else if (surf.rs_method == "sellers") { - return EvaporativeResistanceCoefSellers(surf.saturation_liq); - } else { - throw("Currently support {sakagucki_zeng, sellers}"); - } - } -} - -double -EvaporativeResistanceCoefSakaguckiZeng(double saturation_gas, - double porosity, - double dessicated_zone_thickness, - double Clapp_Horn_b) -{ - double Rsoil; - if (saturation_gas == 0.) { - Rsoil = 0.; // ponded water - } else { - // Equation for reduced vapor diffusivity - // See Sakagucki and Zeng 2009 eqaution (9) and Moldrup et al., 2004. - // - // This really needs to be refactored independently of C&H. There are a - // lot of assumptions here of hard-coded parameters including C&H WRM, a - // residual water content of 0.0556 (not sure why this was chosen), and - // more. - // - // The result of using this with other WRMs requires some adaptation... - // also with arbitrary values. - double s_res = std::min(0.0556 / porosity, 0.4); - double vp_diffusion = - 0.000022 * std::pow(porosity, 2) * std::pow(1 - s_res, 2 + 3 * Clapp_Horn_b); - // Sakagucki and Zeng 2009 eqaution (10) - double L_Rsoil = - dessicated_zone_thickness * (std::exp(std::pow(saturation_gas, 5)) - 1) / (std::exp(1) - 1); - Rsoil = L_Rsoil / vp_diffusion; + return surf.rsoil; } - AMANZI_ASSERT(Rsoil >= 0); - return Rsoil; } -double -EvaporativeResistanceCoefSellers(double saturation_liq) -{ - return std::exp(8.206 - 4.255 * saturation_liq); -} double SensibleHeat(double resistance_coef, @@ -390,7 +345,7 @@ UpdateEnergyBalanceWithoutSnow(const GroundProperties& surf, (params.Da0_a * std::pow(Re0, params.Da0_b) - (params.Cd0_c * std::log(Re0) + params.Cd0_d)) * c_von_Karman; double Dhe_latent = WindFactor(met.Us, met.Z_Us, surf.roughness, KB); - double Rsoil = EvaporativeResistanceGround(surf, met, params, vapor_pressure_skin); + double Rsoil = EvaporativeResistanceGround(surf, met, vapor_pressure_skin); double coef = 1.0 / (Rsoil + 1.0 / (Dhe_latent * Sqig)); diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_funcs.hh b/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_funcs.hh index 3b9b4486a..baf4e9ff8 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_funcs.hh +++ b/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_funcs.hh @@ -122,17 +122,8 @@ VaporPressureGround(const GroundProperties& surf, const ModelParams& params); double EvaporativeResistanceGround(const GroundProperties& surf, const MetData& met, - const ModelParams& params, double vapor_pressure_ground); -double -EvaporativeResistanceCoefSakaguckiZeng(double saturation_gas, - double porosity, - double dessicated_zone_thickness, - double Clapp_Horn_b); - -double -EvaporativeResistanceCoefSellers(double saturation_liq); // // Basic sensible heat. diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/seb_threecomponent_evaluator.cc b/src/pks/surface_balance/constitutive_relations/land_cover/seb_threecomponent_evaluator.cc index 59de50e22..cc1682058 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/seb_threecomponent_evaluator.cc +++ b/src/pks/surface_balance/constitutive_relations/land_cover/seb_threecomponent_evaluator.cc @@ -152,14 +152,10 @@ SEBThreeComponentEvaluator::SEBThreeComponentEvaluator(Teuchos::ParameterList& p dependencies_.insert(KeyTag{ surf_temp_key_, tag }); surf_pres_key_ = Keys::readKey(plist, domain_, "pressure", "pressure"); dependencies_.insert(KeyTag{ surf_pres_key_, tag }); + rsoil_key_ = Keys::readKey(plist, domain_, "soil resistance", "soil_resistance"); + dependencies_.insert(KeyTag{ rsoil_key_, tag }); // -- subsurface properties for evaporating bare soil - sat_gas_key_ = Keys::readKey(plist, domain_ss_, "gas saturation", "saturation_gas"); - dependencies_.insert(KeyTag{ sat_gas_key_, tag }); - sat_liq_key_ = Keys::readKey(plist, domain_ss_, "liquid saturation", "saturation_liquid"); - dependencies_.insert(KeyTag{ sat_liq_key_, tag }); - poro_key_ = Keys::readKey(plist, domain_ss_, "porosity", "porosity"); - dependencies_.insert(KeyTag{ poro_key_, tag }); ss_pres_key_ = Keys::readKey(plist, domain_ss_, "subsurface pressure", "pressure"); dependencies_.insert(KeyTag{ ss_pres_key_, tag }); @@ -207,11 +203,9 @@ SEBThreeComponentEvaluator::Evaluate_(const State& S, const std::vector(area_frac_key_, tag).ViewComponent("cell", false); const auto& surf_pres = *S.Get(surf_pres_key_, tag).ViewComponent("cell", false); const auto& surf_temp = *S.Get(surf_temp_key_, tag).ViewComponent("cell", false); + const auto& rsoil = *S.Get(rsoil_key_, tag).ViewComponent("cell", false); // collect subsurface properties - const auto& sat_gas = *S.Get(sat_gas_key_, tag).ViewComponent("cell", false); - const auto& sat_liq = *S.Get(sat_liq_key_, tag).ViewComponent("cell", false); - const auto& poro = *S.Get(poro_key_, tag).ViewComponent("cell", false); const auto& ss_pres = *S.Get(ss_pres_key_, tag).ViewComponent("cell", false); // collect output vecs @@ -285,15 +279,10 @@ SEBThreeComponentEvaluator::Evaluate_(const State& S, const std::vector(surf_temp_key_, tag).ViewComponent("cell", false); // collect subsurface properties - const auto& sat_gas = *S.Get(sat_gas_key_, tag).ViewComponent("cell", false); - const auto& sat_liq = *S.Get(sat_liq_key_, tag).ViewComponent("cell", false); - const auto& poro = *S.Get(poro_key_, tag).ViewComponent("cell", false); const auto& ss_pres = *S.Get(ss_pres_key_, tag).ViewComponent("cell", false); + const auto& rsoil = *S.Get(rsoil_key_, tag).ViewComponent("cell", false); // collect output vecs auto& water_source = *results[0]->ViewComponent("cell", false); @@ -284,15 +278,11 @@ SEBTwoComponentEvaluator::Evaluate_(const State& S, const std::vector lc.second.water_transition_depth) { surf.pressure = surf_pres[0][c]; - surf.porosity = 1.; - surf.saturation_gas = 0.; - surf.saturation_liq = sat_liq[0][cells[0]]; + surf.rsoil = 0.; } else { double factor = std::max(ponded_depth[0][c], 0.) / lc.second.water_transition_depth; surf.pressure = factor * surf_pres[0][c] + (1 - factor) * ss_pres[0][cells[0]]; - surf.porosity = factor + (1 - factor) * poro[0][cells[0]]; - surf.saturation_gas = (1 - factor) * sat_gas[0][cells[0]]; - surf.saturation_liq = factor + (1 - factor) * sat_liq[0][cells[0]]; + surf.rsoil = (1 - factor) * rsoil[0][c]; } if (model_1p1_) surf.pressure = surf_pres[0][c]; surf.ponded_depth = ponded_depth[0][c]; @@ -302,9 +292,6 @@ SEBTwoComponentEvaluator::Evaluate_(const State& S, const std::vectorSetGhosted()->AddComponent("cell", AmanziMesh::CELL, 1); diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/seb_twocomponent_evaluator.hh b/src/pks/surface_balance/constitutive_relations/land_cover/seb_twocomponent_evaluator.hh index bf7c4475d..f1a45de58 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/seb_twocomponent_evaluator.hh +++ b/src/pks/surface_balance/constitutive_relations/land_cover/seb_twocomponent_evaluator.hh @@ -148,7 +148,7 @@ class SEBTwoComponentEvaluator : public EvaluatorSecondaryMonotypeCV { Key ponded_depth_key_, unfrozen_fraction_key_; Key sg_albedo_key_, sg_emissivity_key_, area_frac_key_; Key surf_temp_key_, surf_pres_key_; - Key sat_gas_key_, sat_liq_key_, poro_key_, ss_pres_key_; + Key ss_pres_key_, rsoil_key_; Key mol_dens_key_, mass_dens_key_; Key melt_key_, evap_key_;