[WIP] Refactor of ELM-ATS API

src/constitutive_relations/CMakeLists.txt

@@ -44,6 +44,11 @@ register_evaluator_with_factory(
   LISTNAME   ATS_RELATIONS_REG
   )
 
+register_evaluator_with_factory(
+  HEADERFILE surface_subsurface_fluxes/flux_divergence_evaluator_reg.hh
+  LISTNAME   ATS_RELATIONS_REG
+  )
+
 
 register_evaluator_with_factory(
   HEADERFILE generic_evaluators/MultiplicativeEvaluator_reg.hh
@@ -80,6 +85,11 @@ register_evaluator_with_factory(
   LISTNAME ATS_RELATIONS_REG
   )
 
+register_evaluator_with_factory(
+  HEADERFILE generic_evaluators/EvaluatorTimeAccumulated_reg.hh
+  LISTNAME ATS_RELATIONS_REG
+  )
+
 register_evaluator_with_factory(
   HEADERFILE generic_evaluators/ExtractionEvaluator_reg.hh
   LISTNAME ATS_RELATIONS_REG

src/constitutive_relations/generic_evaluators/CMakeLists.txt

@@ -12,6 +12,7 @@ set(ats_generic_evals_src_files
     TimeMaxEvaluator.cc
     ExtractionEvaluator.cc
     InitialTimeEvaluator.cc
+    EvaluatorTimeAccumulated.cc
    )
 
 file(GLOB ats_generic_evals_inc_files "*.hh")

src/constitutive_relations/generic_evaluators/EvaluatorTimeAccumulated.cc

@@ -0,0 +1,165 @@
+/*
+  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 ([email protected])
+*/
+
+#include "EvaluatorTimeAccumulated.hh"
+
+namespace Amanzi {
+namespace Relations {
+
+EvaluatorTimeAccumulated::EvaluatorTimeAccumulated(Teuchos::ParameterList& plist)
+  : EvaluatorSecondary(plist)
+{
+  // my_keys_[0] was populated by EvaluatorSecondary from the plist name/tag.
+  const auto& my_key = my_keys_.front().first;
+  const auto& my_tag = my_keys_.front().second;
+
+  accumulation_type_ = plist_.get<std::string>("accumulation type", "integral");
+  if (accumulation_type_ != "integral" && accumulation_type_ != "min" &&
+      accumulation_type_ != "max") {
+    Errors::Message msg;
+    msg << "EvaluatorTimeAccumulated for \"" << my_key << "\": invalid \"accumulation type\" \""
+        << accumulation_type_ << "\", must be one of: integral, min, max";
+    Exceptions::amanzi_throw(msg);
+  }
+
+  auto domain = Keys::getDomain(my_key);
+  accumulated_key_ = Keys::readKey(plist_, domain, "accumulated", Keys::getVarName(my_key));
+
+  if (!plist_.isParameter("accumulated tag")) {
+    Errors::Message msg;
+    msg << "EvaluatorTimeAccumulated for \"" << my_key
+        << "\": missing required parameter \"accumulated tag\"";
+    Exceptions::amanzi_throw(msg);
+  }
+  accumulated_tag_ = Tag(plist_.get<std::string>("accumulated tag"));
+
+  // my_keys_[1]: the accumulated-dt scalar, same tag as the CV result
+  Key acc_dt_key = my_key + "_accumulated_dt";
+  my_keys_.emplace_back(KeyTag{ acc_dt_key, my_tag });
+
+  // accumulated_key@accumulated_tag changes each inner step — triggers Update_()
+  dependencies_.insert(KeyTag(accumulated_key_, accumulated_tag_));
+
+  // no evaluator for dt, so can't add it to dependencies
+  // dt@accumulated_tag provides the inner timestep size for weighting
+  //dependencies_.insert(KeyTag("dt", accumulated_tag_));
+}
+
+
+Teuchos::RCP<Evaluator>
+EvaluatorTimeAccumulated::Clone() const
+{
+  return Teuchos::rcp(new EvaluatorTimeAccumulated(*this));
+}
+
+
+void
+EvaluatorTimeAccumulated::EnsureCompatibility(State& S)
+{
+  // my_keys_[0]: CV result — claim ownership; structure flows from client requirements
+  const auto& [cv_key, cv_tag] = my_keys_[0];
+  auto& my_fac = S.Require<CompositeVector, CompositeVectorSpace>(cv_key, cv_tag, cv_key);
+
+  // my_keys_[1]: accumulated-dt scalar — claim ownership
+  const auto& [dt_key, dt_tag] = my_keys_[1];
+  S.Require<double>(dt_key, dt_tag, dt_key);
+
+  // wire dependency evaluators into the graph
+  if (my_fac.Mesh() != Teuchos::null) {
+    for (const auto& dep : dependencies_)
+      // dependency fac includes my fac
+      S.Require<CompositeVector, CompositeVectorSpace>(dep.first, dep.second)
+        .Update(my_fac);
+  }
+
+  EvaluatorSecondary::EnsureCompatibility_DepEnsureCompatibility_(S);
+  EvaluatorSecondary::EnsureCompatibility_Flags_(S);
+}
+
+
+void
+EvaluatorTimeAccumulated::Update_(State& S)
+{
+  const auto& [cv_key, cv_tag] = my_keys_[0];
+  const auto& [dt_key, dt_tag] = my_keys_[1];
+
+  double dt_inner = S.Get<double>("dt", accumulated_tag_);
+  const auto& source = S.Get<CompositeVector>(accumulated_key_, accumulated_tag_);
+  double acc_dt = S.Get<double>(dt_key, dt_tag);
+  if (dt_inner > 0) {
+    // cv is read and overwritten in place — valid because we own it
+    auto& cv = S.GetW<CompositeVector>(cv_key, cv_tag, cv_key);
+
+    if (accumulation_type_ == "integral") {
+      double total_dt = acc_dt + dt_inner;
+      for (const auto& comp : cv) {
+        auto& res = *cv.ViewComponent(comp, false);
+        const auto& src = *source.ViewComponent(comp, false);
+        for (int j = 0; j != res.NumVectors(); ++j)
+          for (int i = 0; i != res.MyLength(); ++i)
+            res[j][i] = (res[j][i] * acc_dt + src[j][i] * dt_inner) / total_dt;
+      }
+    } else if (accumulation_type_ == "max") {
+      for (const auto& comp : cv) {
+        auto& res = *cv.ViewComponent(comp, false);
+        const auto& src = *source.ViewComponent(comp, false);
+        for (int j = 0; j != res.NumVectors(); ++j)
+          for (int i = 0; i != res.MyLength(); ++i)
+            res[j][i] = std::max(res[j][i], src[j][i]);
+      }
+    } else { // min
+      for (const auto& comp : cv) {
+        auto& res = *cv.ViewComponent(comp, false);
+        const auto& src = *source.ViewComponent(comp, false);
+        for (int j = 0; j != res.NumVectors(); ++j)
+          for (int i = 0; i != res.MyLength(); ++i)
+            res[j][i] = std::min(res[j][i], src[j][i]);
+      }
+    }
+
+    S.GetW<double>(dt_key, dt_tag, dt_key) = acc_dt + dt_inner;
+
+    if (vo_.os_OK(Teuchos::VERB_EXTREME)) {
+      Teuchos::OSTab tab = vo_.getOSTab();
+      *vo_.os() << "Updated time accumulation time from " << acc_dt << " to " << acc_dt + dt_inner << std::endl;
+      cv.Print(std::cout);
+    }
+  }
+}
+
+
+void
+EvaluatorTimeAccumulated::Reset(State& S)
+{
+  const auto& [cv_key, cv_tag] = my_keys_[0];
+  const auto& [dt_key, dt_tag] = my_keys_[1];
+
+  if (vo_.os_OK(Teuchos::VERB_EXTREME)) {
+    Teuchos::OSTab tab = vo_.getOSTab();
+    *vo_.os() << "Resetting time accumulation" << std::endl;
+  }
+
+  auto& cv = S.GetW<CompositeVector>(cv_key, cv_tag, cv_key);
+  if (accumulation_type_ == "max") {
+    cv.PutScalar(-1.e16);
+  } else if (accumulation_type_ == "min") {
+    cv.PutScalar(1.e16);
+  } else {
+    cv.PutScalar(0.);
+  }
+
+  S.GetW<double>(dt_key, dt_tag, dt_key) = 0.;
+
+  // invalidate lazy cache so the next Update() call always runs Update_()
+  requests_.clear();
+}
+
+
+} // namespace Relations
+} // namespace Amanzi

src/constitutive_relations/generic_evaluators/EvaluatorTimeAccumulated.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: Ethan Coon ([email protected])
+*/
+
+/*!
+
+Accumulates a field pointwise over time, computing either a time-integrated
+average, a running maximum, or a running minimum.
+
+This evaluator is designed for use with subcycled MPCs.  It lives at
+``tag`` (typically ``Tags::NEXT``) and depends on the source field at
+``accumulated_tag`` (the subcycled inner tag, e.g. ``flow_next``).  The
+``TimeAdvancer`` driving the inner subcycle loop is responsible for calling
+``Reset()`` at the start of each outer timestep and ``Update()`` after each
+successful inner step.
+
+For ``integral`` mode the result is the time-weighted mean:
+
+  new_val = (old_val * acc_dt + source * dt) / (acc_dt + dt)
+
+For ``min`` / ``max`` modes the result is the running pointwise
+minimum / maximum; accumulated time is still tracked so that ``Reset()``
+can initialise correctly.
+
+Both the result ``CompositeVector`` and the ``KEY_accumulated_dt`` scalar are
+checkpointed so that the accumulation survives restarts mid outer-step.
+
+`"evaluator type`" = `"time accumulated`"
+
+.. _evaluator-time-accumulated-spec:
+.. admonition:: evaluator-time-accumulated-spec
+
+   * `"accumulation type`" ``[string]`` **"integral"** One of ``integral``,
+     ``min``, or ``max``.
+
+   * `"accumulated key`" ``[string]`` **my_key** Key of the field to
+     accumulate.  Defaults to the same key as this evaluator.
+
+   * `"accumulated tag`" ``[string]`` **required** Tag at which the source
+     field is evaluated (the subcycled inner tag).
+
+*/
+
+#pragma once
+
+#include "Factory.hh"
+#include "EvaluatorSecondary.hh"
+
+namespace Amanzi {
+namespace Relations {
+
+class EvaluatorTimeAccumulated : public EvaluatorSecondary {
+ public:
+  explicit EvaluatorTimeAccumulated(Teuchos::ParameterList& plist);
+  EvaluatorTimeAccumulated(const EvaluatorTimeAccumulated& other) = default;
+  Teuchos::RCP<Evaluator> Clone() const override;
+
+  bool IsDifferentiableWRT(const State& S,
+                           const Key& wrt_key,
+                           const Tag& wrt_tag) const override
+  {
+    return false;
+  }
+
+  void EnsureCompatibility(State& S) override;
+
+  // Zeros the accumulator and clears the lazy-evaluation cache.
+  // Called by TimeAdvancer before each outer timestep's inner loop.
+  void Reset(State& S);
+
+ protected:
+  void Update_(State& S) override;
+  void UpdateDerivative_(State& S, const Key& wrt_key, const Tag& wrt_tag) override {}
+
+ protected:
+  Key accumulated_key_;
+  Tag accumulated_tag_;
+  std::string accumulation_type_;  // "integral", "min", "max"
+
+  // my_keys_[0] = {key, tag}         — the CV result
+  // my_keys_[1] = {acc_dt_key, tag}  — the accumulated-dt scalar
+
+ private:
+  static Utils::RegisteredFactory<Evaluator, EvaluatorTimeAccumulated> reg_;
+};
+
+} // namespace Relations
+} // namespace Amanzi

src/constitutive_relations/generic_evaluators/EvaluatorTimeAccumulated_reg.hh

@@ -0,0 +1,19 @@
+/*
+  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 ([email protected])
+*/
+
+#include "EvaluatorTimeAccumulated.hh"
+
+namespace Amanzi {
+namespace Relations {
+
+Utils::RegisteredFactory<Evaluator, EvaluatorTimeAccumulated>
+  EvaluatorTimeAccumulated::reg_("time accumulated");
+
+} // namespace Relations
+} // namespace Amanzi

src/constitutive_relations/surface_subsurface_fluxes/CMakeLists.txt

@@ -10,6 +10,7 @@ set(ats_surf_subsurf_src_files
   surface_top_cells_evaluator.cc
   top_cells_surface_evaluator.cc
   volumetric_darcy_flux_evaluator.cc
+  flux_divergence_evaluator.cc
  )
 
 file(GLOB ats_surf_subsurf_inc_files "*.hh")

src/constitutive_relations/surface_subsurface_fluxes/flux_divergence_evaluator.cc

@@ -0,0 +1,65 @@
+/*
+  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  ([email protected])
+*/
+
+#include "flux_divergence_evaluator.hh"
+
+namespace Amanzi {
+namespace Relations {
+
+FluxDivergenceEvaluator::FluxDivergenceEvaluator(Teuchos::ParameterList& plist)
+  : EvaluatorSecondaryMonotypeCV(plist)
+{
+  // determine the domain
+  Key akey = my_keys_.front().first;
+  auto tag = my_keys_.front().second;
+  Key domain = Keys::getDomain(akey);
+
+  flux_key_ = Keys::readKey(plist, domain, "flux", "water_flux");
+  dependencies_.insert(KeyTag{ flux_key_, tag });
+}
+
+Teuchos::RCP<Evaluator>
+FluxDivergenceEvaluator::Clone() const
+{
+  return Teuchos::rcp(new FluxDivergenceEvaluator(*this));
+}
+
+void
+FluxDivergenceEvaluator::EnsureCompatibility_ToDeps_(State& S)
+{
+  const auto& my_fac = S.Require<CompositeVector, CompositeVectorSpace>(my_keys_.front().first, my_keys_.front().second);
+  if (my_fac.Mesh() != Teuchos::null) {
+    for (auto dep : dependencies_) {
+      auto& fac = S.Require<CompositeVector, CompositeVectorSpace>(dep.first, dep.second);
+      fac.SetMesh(my_fac.Mesh())->AddComponent("face", AmanziMesh::Entity_kind::FACE, 1);
+    }
+  }
+}
+
+
+void
+FluxDivergenceEvaluator::Evaluate_(const State& S, const std::vector<CompositeVector*>& result)
+{
+  auto tag = my_keys_.front().second;
+  const auto& flux_f = *S.Get<CompositeVector>(flux_key_, tag).ViewComponent("face", true);
+  const AmanziMesh::Mesh& m = *result[0]->Mesh();
+  auto& div_flux = *result[0]->ViewComponent("cell", false);
+  div_flux.PutScalar(0.);
+
+  for (AmanziMesh::Entity_ID c = 0; c != div_flux.MyLength(); ++c) {
+    auto [cfaces, cfdirs] = m.getCellFacesAndDirections(c);
+    for (int i = 0; i != cfaces.size(); ++i) {
+      div_flux[0][c] += flux_f[0][cfaces[i]] * cfdirs[i];
+    }
+  }
+}
+
+
+} // namespace Relations
+} // namespace Amanzi