Move calculation of soil resistance from landcover to sources; move W…

…RM related parameters from evaluators to model parameters

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<SoilResistanceModelPartition>
+createSoilResistanceModelPartition(Teuchos::ParameterList& plist)
+{
+  SoilResistanceSakaguckiZengModelList rs_list;
+  std::vector<std::string> 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<std::string>("region"));
+      rs_list.push_back(Teuchos::rcp(new SoilResistanceSakaguckiZengModel(sublist)));
+    } else {
+      AMANZI_ASSERT(0);
+    }
+  }
+
+  Teuchos::RCP<Functions::MeshPartition> part =
+    Teuchos::rcp(new Functions::MeshPartition(AmanziMesh::CELL, region_list));
+
+  return Teuchos::rcp(new SoilResistanceModelPartition(part, rs_list));
+}
+
+} // namespace Flow
+} // namespace Amanzi

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<Teuchos::RCP<SoilResistanceSakaguckiZengModel>> 
+  SoilResistanceSakaguckiZengModelList;
+typedef std::pair<Teuchos::RCP<Functions::MeshPartition>, SoilResistanceSakaguckiZengModelList> 
+  SoilResistanceModelPartition;
+
+// Non-member factory
+Teuchos::RCP<SoilResistanceModelPartition>
+createSoilResistanceModelPartition(Teuchos::ParameterList& plist);
+
+} // namespace Flow
+} // namespace Amanzi
+
+#endif

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<Evaluator>
+SoilResistanceSakaguckiZengEvaluator::Clone() const
+{
+  return Teuchos::rcp(new SoilResistanceSakaguckiZengEvaluator(*this));
+}
+
+
+// Required methods from EvaluatorSecondaryMonotypeCV
+void
+SoilResistanceSakaguckiZengEvaluator::Evaluate_(const State& S,
+                   const std::vector<CompositeVector*>& 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<const CompositeVector> sat_gas = S.GetPtr<CompositeVector>(sat_gas_key_, tag);
+  Teuchos::RCP<const CompositeVector> poro = S.GetPtr<CompositeVector>(poro_key_, tag);
+  Teuchos::RCP<const AmanziMesh::Mesh> mesh = sat_gas->Mesh();
+  Teuchos::RCP<const AmanziMesh::Mesh> 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<CompositeVector*>& 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<const CompositeVector> sat_gas = S.GetPtr<CompositeVector>(sat_gas_key_, tag);
+  Teuchos::RCP<const CompositeVector> poro = S.GetPtr<CompositeVector>(poro_key_, tag);
+  Teuchos::RCP<const AmanziMesh::Mesh> mesh = sat_gas->Mesh();
+  Teuchos::RCP<const AmanziMesh::Mesh> 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<CompositeVector, CompositeVectorSpace>(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<CompositeVector, CompositeVectorSpace>(dep.first, dep.second).Update(fac);
+      } else {
+        S.Require<CompositeVector, CompositeVectorSpace>(dep.first, dep.second).Update(dep_fac);
+      }
+    }
+  }
+}
+
+
+} // namespace Flow
+} // namespace Amanzi

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<Evaluator> Clone() const override;
+
+  Teuchos::RCP<SoilResistanceModelPartition> get_Models() { return models_; }
+
+ protected:
+  // Required methods from EvaluatorSecondaryMonotypeCV
+  virtual void Evaluate_(const State& S, const std::vector<CompositeVector*>& result) override;
+  virtual void EvaluatePartialDerivative_(const State& S,
+                                          const Key& wrt_key,
+                                          const Tag& wrt_tag,
+                                          const std::vector<CompositeVector*>& result) override;
+
+  virtual void EnsureCompatibility_ToDeps_(State& S) override;
+
+ protected:
+  Key sat_gas_key_;
+  Key poro_key_;
+
+  Teuchos::RCP<SoilResistanceModelPartition> models_;
+
+ private:
+  static Utils::RegisteredFactory<Evaluator, SoilResistanceSakaguckiZengEvaluator> fac_;
+};
+
+} // namespace Flow
+} // namespace Amanzi
+
+#endif

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<Evaluator, SoilResistanceSakaguckiZengEvaluator>
+  SoilResistanceSakaguckiZengEvaluator::fac_("sakagucki-zeng soil resistance");
+
+} // namespace Flow
+} // namespace Amanzi