Added third pk to mpc_flow_transport

src/pks/mpc/mpc_flow_transport.cc

@@ -252,6 +252,62 @@ MPCFlowTransport::parseParameterList()
       requireEvaluatorAtCurrent("surface-porosity", transport_current_tag, *S_);
     }
   }
+
+  if (sub_pks_.size() == 3) {
+    Key gas_sub_lwc_key = Keys::readKey(*getSubPKPlist_(2), "subsurface", "liquid water content", "water_content");
+    auto [gas_transport_current_tag, gas_transport_next_tag] = tags_[2];
+    if (gas_transport_next_tag != flow_next_tag) {
+      // set the flow field evaluator as the flow's NEXT tag
+      //
+      // Note, we could be more careful here and readKey() the flow field's name
+      // from the flow PK's sublist (which may be nested two deep).  Instead we
+      // hard-code this as the default.  If this breaks in the future it can be
+      // fixed. --ETC
+      Teuchos::ParameterList& flux_list =
+        S_->GetEvaluatorList(Keys::getKey("water_flux", gas_transport_next_tag));
+      if (!flux_list.isParameter("evaluator type")) {
+        flux_list.set<std::string>("evaluator type", "alias");
+        flux_list.set<std::string>("target", Keys::getKey("water_flux", flow_next_tag, true));
+      }
+
+      // velocity for dispersivity
+      Teuchos::ParameterList& velo_list =
+        S_->GetEvaluatorList(Keys::getKey("darcy_velocity", gas_transport_next_tag));
+      if (!velo_list.isParameter("evaluator type")) {
+        velo_list.set<std::string>("evaluator type", "alias");
+        velo_list.set<std::string>("target", Keys::getKey("darcy_velocity", flow_next_tag, true));
+      }
+
+      // now set the liquid water content as an interpolated field at next
+      // note that flow_current copy is kept by flow PK, and transport_current copy is kept by transport PK
+      Teuchos::ParameterList& lwc_list_next =
+        S_->GetEvaluatorList(Keys::getKey(gas_sub_lwc_key, gas_transport_next_tag));
+      if (!lwc_list_next.isParameter("evaluator type")) {
+        lwc_list_next.set<std::string>("evaluator type", "temporal interpolation");
+        lwc_list_next.set<std::string>("current tag", flow_current_tag.get());
+        lwc_list_next.set<std::string>("next tag", flow_next_tag.get());
+      }
+
+      // porosity used with velocity to compute particle velocity when dispersion is on
+      // -- and an interpolation at transport's next
+      Teuchos::ParameterList& poro_list_next =
+        S_->GetEvaluatorList(Keys::getKey("porosity", gas_transport_next_tag));
+      if (!poro_list_next.isParameter("evaluator type")) {
+        poro_list_next.set<std::string>("evaluator type", "temporal interpolation");
+        poro_list_next.set<std::string>("current tag", flow_current_tag.get());
+        poro_list_next.set<std::string>("next tag", flow_next_tag.get());
+      }
+    }
+
+    requireEvaluatorAtNext(gas_sub_lwc_key, flow_next_tag, *S_);
+    requireEvaluatorAtNext("porosity", flow_next_tag, *S_);
+    requireEvaluatorAtCurrent("porosity", flow_current_tag, *S_, name_);
+
+    // now require key@transport_next, which will be the interpolant
+    requireEvaluatorAtNext(gas_sub_lwc_key, gas_transport_next_tag, *S_);
+    requireEvaluatorAtNext("porosity", gas_transport_next_tag, *S_);
+
+  }
 }
 
 

src/pks/transport/constitutive_relations/sources/qc_relation_divide_evaluator.cc

@@ -0,0 +1,64 @@
+/*
+  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: Phong V.V. Le (lepv@ornl.gov)
+*/
+#include <algorithm>
+#include "Key.hh"
+#include "Factory.hh"
+#include "Function.hh"
+#include "qc_relation_divide_evaluator.hh"
+#include "FunctionFactory.hh"
+
+namespace Amanzi {
+namespace Flow {
+namespace Relations {
+
+
+QCRelationDivideEvaluator::QCRelationDivideEvaluator(Teuchos::ParameterList& plist) : EvaluatorSecondaryMonotypeCV(plist)
+{
+  domain_ = Keys::getDomain(my_keys_.front().first);
+  auto tag = my_keys_.front().second;
+
+  cv_key_ = Keys::readKey(plist, domain_, "cell volume", "cell_volume");
+  dependencies_.insert(KeyTag{ cv_key_, tag });
+  molar_density_key_ = Keys::readKey(plist, domain_, "molar density liquid", "molar_density_liquid");
+  dependencies_.insert(KeyTag{ molar_density_key_, tag });
+  first_src_key_ = Keys::readKey(plist, domain_, "first source", "first_source");
+  dependencies_.insert(KeyTag{ first_src_key_, tag });
+  second_src_key_ = Keys::readKey(plist, domain_, "second source", "second_source");
+  dependencies_.insert(KeyTag{ second_src_key_, tag });
+}
+
+// Required methods from SecondaryVariableDivideEvaluator
+void
+QCRelationDivideEvaluator::Evaluate_(const State& S, const std::vector<CompositeVector*>& result) 
+{
+
+  Tag tag = my_keys_.front().second;
+
+  const auto& cv = *S.Get<CompositeVector>(cv_key_, tag).ViewComponent("cell", false);
+  const auto& molar_den =
+    *S.Get<CompositeVector>(molar_density_key_, tag).ViewComponent("cell", false);
+  const auto& first_source =
+    *S.Get<CompositeVector>(first_src_key_, tag).ViewComponent("cell", false);
+  const auto& second_source =
+    *S.Get<CompositeVector>(second_src_key_, tag).ViewComponent("cell", false);    
+
+  auto& surf_src = *result[0]->ViewComponent("cell"); // not being reference
+  const AmanziMesh::Mesh& mesh = *result[0]->Mesh();
+
+  // Loop through each cell
+  AmanziMesh::Entity_ID ncells = cv.MyLength();
+  for (AmanziMesh::Entity_ID c = 0; c != ncells; ++c) {  
+    // transport source as a function of discharge and another variable
+    surf_src[0][c] = first_source[0][c] / second_source[0][c];
+  }
+}
+
+} // namespace Relations
+} // namespace Flow
+} // namespace Amanzi 

src/pks/transport/constitutive_relations/sources/qc_relation_divide_evaluator.hh

@@ -0,0 +1,93 @@
+/*
+  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: Phong V.V. Le (lepv@ornl.gov)
+*/
+
+//! Evaluates transport source (mass) from a ponded depth
+/*!
+
+Mass sources into stream/river from a ponded depth
+
+.. _qc_relation_depth_evaluator-spec:
+.. admonition:: qc_relation_depth_evaluator-spec
+
+   `"function`" ``[function-spec]`` Function describing the relationship between ponded depth (e.g. tile, groundwater) and solute mass running into river/stream
+
+   KEYS:
+   - `"cell volume`" **DOMAIN-cell_volume**
+   - `"molar density liquid`" **DOMAIN-molar_density_liquid** 
+   - `"depth source`" **DOMAIN-depth_source** source
+   - `"extensive`" ``[bool]`` checks if source is extensive. Default value is *false*.
+
+Example
+
+.. code-block:: xml
+
+  <ParameterList name="RIVER_DOMAIN-depth_sources">
+    <Parameter name="evaluator type" type="string" value="q-c depth" />
+    <Parameter name="depth source key" type="string" value="RIVER_DOMAIN-water_source_depth" />
+    <Parameter name="extensive" type="bool" value="false" />
+    <ParameterList name="function" type="ParameterList">
+      <ParameterList name="function-tabular" type="ParameterList">
+        <Parameter name="x values" type="Array(double)" value="{0, 0.1, 1, 10}" />
+        <Parameter name="y values" type="Array(double)" value="{0, 0.08, 0.1, 0.2}" /> 
+        <Parameter name="forms" type="Array(string)" value="{linear, linear, linear}" />
+      </ParameterList>
+    </ParameterList>
+  </ParameterList>
+
+*/
+
+#pragma once
+
+#include "Factory.hh"
+#include "EvaluatorSecondaryMonotype.hh"
+#include "FunctionFactory.hh"
+
+namespace Amanzi {
+namespace Flow {
+namespace Relations {
+
+class QCRelationDivideEvaluator : public EvaluatorSecondaryMonotypeCV {
+ public:
+  explicit QCRelationDivideEvaluator(Teuchos::ParameterList& plist);
+  QCRelationDivideEvaluator(const QCRelationDivideEvaluator& other) = default;
+
+  virtual Teuchos::RCP<Evaluator> Clone() const override
+  {
+    return Teuchos::rcp(new QCRelationDivideEvaluator(*this));
+  }
+
+  // virtual void EnsureCompatibility(State& S) override;
+  virtual bool
+  IsDifferentiableWRT(const State& S, const Key& wrt_key, const Tag& wrt_tag) const override
+  {
+    return false;
+  }
+
+ protected:
+  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{};
+
+ protected:
+  Key domain_;
+  Key cv_key_;
+  Key molar_density_key_;
+  Key first_src_key_;
+  Key second_src_key_;
+
+ private:
+  static Utils::RegisteredFactory<Evaluator, QCRelationDivideEvaluator> reg_;
+};
+
+} // namespace Relations
+} // namespace Flow
+} // namespace Amanzi
+

src/pks/transport/constitutive_relations/sources/qc_relation_divide_evaluator_reg.hh

@@ -0,0 +1,21 @@
+/*
+  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: Phong V.V. Le (lepv@ornl.gov)
+*/
+
+#include "qc_relation_divide_evaluator.hh"
+
+namespace Amanzi {
+namespace Flow {
+namespace Relations {
+
+Utils::RegisteredFactory<Evaluator, QCRelationDivideEvaluator>
+  QCRelationDivideEvaluator::reg_("q-c divide");
+
+} // namespace Relations
+} // namespace Flow
+} // namespace Amanzi