From 89e43d46f4cb2132666348eea411ef8f545f3c99 Mon Sep 17 00:00:00 2001 From: gaobhub Date: Fri, 12 Jun 2026 17:53:42 -0400 Subject: [PATCH 1/6] add third gas diffusive transport pk aliased evaluators, add phong's q-c divide evaluator --- src/pks/mpc/mpc_flow_transport.cc | 56 +++++++++++ .../sources/qc_relation_divide_evaluator.cc | 64 +++++++++++++ .../sources/qc_relation_divide_evaluator.hh | 93 +++++++++++++++++++ .../qc_relation_divide_evaluator_reg.hh | 21 +++++ 4 files changed, 234 insertions(+) create mode 100644 src/pks/transport/constitutive_relations/sources/qc_relation_divide_evaluator.cc create mode 100644 src/pks/transport/constitutive_relations/sources/qc_relation_divide_evaluator.hh create mode 100644 src/pks/transport/constitutive_relations/sources/qc_relation_divide_evaluator_reg.hh diff --git a/src/pks/mpc/mpc_flow_transport.cc b/src/pks/mpc/mpc_flow_transport.cc index a67d51b0d..0388f7718 100644 --- a/src/pks/mpc/mpc_flow_transport.cc +++ b/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("evaluator type", "alias"); + flux_list.set("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("evaluator type", "alias"); + velo_list.set("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("evaluator type", "temporal interpolation"); + lwc_list_next.set("current tag", flow_current_tag.get()); + lwc_list_next.set("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("evaluator type", "temporal interpolation"); + poro_list_next.set("current tag", flow_current_tag.get()); + poro_list_next.set("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_); + + } } diff --git a/src/pks/transport/constitutive_relations/sources/qc_relation_divide_evaluator.cc b/src/pks/transport/constitutive_relations/sources/qc_relation_divide_evaluator.cc new file mode 100644 index 000000000..4c407fc1a --- /dev/null +++ b/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 +#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& result) +{ + + Tag tag = my_keys_.front().second; + + const auto& cv = *S.Get(cv_key_, tag).ViewComponent("cell", false); + const auto& molar_den = + *S.Get(molar_density_key_, tag).ViewComponent("cell", false); + const auto& first_source = + *S.Get(first_src_key_, tag).ViewComponent("cell", false); + const auto& second_source = + *S.Get(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 \ No newline at end of file diff --git a/src/pks/transport/constitutive_relations/sources/qc_relation_divide_evaluator.hh b/src/pks/transport/constitutive_relations/sources/qc_relation_divide_evaluator.hh new file mode 100644 index 000000000..14a3f4650 --- /dev/null +++ b/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 + + + + + + + + + + + + + + +*/ + +#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 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& result) override; + virtual void EvaluatePartialDerivative_(const State& S, + const Key& wrt_key, + const Tag& wrt_tag, + const std::vector& result) override{}; + + protected: + Key domain_; + Key cv_key_; + Key molar_density_key_; + Key first_src_key_; + Key second_src_key_; + + private: + static Utils::RegisteredFactory reg_; +}; + +} // namespace Relations +} // namespace Flow +} // namespace Amanzi + diff --git a/src/pks/transport/constitutive_relations/sources/qc_relation_divide_evaluator_reg.hh b/src/pks/transport/constitutive_relations/sources/qc_relation_divide_evaluator_reg.hh new file mode 100644 index 000000000..0b6643144 --- /dev/null +++ b/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 + QCRelationDivideEvaluator::reg_("q-c divide"); + +} // namespace Relations +} // namespace Flow +} // namespace Amanzi From 01366c0c80ba34d6654a913bb84084c929ed341e Mon Sep 17 00:00:00 2001 From: gaobhub Date: Sat, 13 Jun 2026 00:27:17 -0400 Subject: [PATCH 2/6] require primary_free_ion_concentration from alquimia --- src/pks/transport/transport_ats_pk.cc | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/pks/transport/transport_ats_pk.cc b/src/pks/transport/transport_ats_pk.cc index fdb27fd34..2d41cf6f1 100644 --- a/src/pks/transport/transport_ats_pk.cc +++ b/src/pks/transport/transport_ats_pk.cc @@ -580,6 +580,16 @@ Transport_ATS::SetupPhysicalEvaluators_() // ->AddComponent("face", AmanziMesh::Entity_kind::FACE, num_components_); S_->GetRecordSetW(key_).set_subfieldnames(component_names_); + if (chem_engine_ != Teuchos::null) { + Key primary_free_ion_concentration_key_ = + Keys::readKey(*plist_, domain_, "primary free ion concentration", "primary_free_ion_concentration"); + + S_->Require(primary_free_ion_concentration_key_, tag_next_, passwd_) + .SetMesh(mesh_) + ->SetGhosted(true) + ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, num_components_); + } + // -- water flux requireEvaluatorAtNext(water_flux_key_, tag_next_, *S_) .SetMesh(mesh_) From 2e2598cb32ca69d7d7bf003ca730b404082d8d7c Mon Sep 17 00:00:00 2001 From: gaobhub Date: Sat, 13 Jun 2026 18:52:06 -0400 Subject: [PATCH 3/6] try to require Amanzi primary free ion for source calculation in ats transport --- .../generic_evaluators/MultiplicativeEvaluator.cc | 3 +++ src/pks/transport/transport_ats.hh | 1 + src/pks/transport/transport_ats_pk.cc | 11 ++++++----- 3 files changed, 10 insertions(+), 5 deletions(-) diff --git a/src/constitutive_relations/generic_evaluators/MultiplicativeEvaluator.cc b/src/constitutive_relations/generic_evaluators/MultiplicativeEvaluator.cc index 6952f074c..f58e66ade 100644 --- a/src/constitutive_relations/generic_evaluators/MultiplicativeEvaluator.cc +++ b/src/constitutive_relations/generic_evaluators/MultiplicativeEvaluator.cc @@ -62,6 +62,9 @@ MultiplicativeEvaluator::Evaluate_(const State& S, const std::vectorViewComponent(lcv_name, false)); int i = 0; for (const auto& key_tag : dependencies_) { + std::cout << "DEBUG MultiplicativeEvaluator::Evaluate_ for " + << my_keys_[0].first << std::endl; + const auto& dep_v = *(*S.Get(key_tag.first, key_tag.second).ViewComponent(lcv_name, false))( dofs_[i]); diff --git a/src/pks/transport/transport_ats.hh b/src/pks/transport/transport_ats.hh index 78aef4d41..e3160795e 100644 --- a/src/pks/transport/transport_ats.hh +++ b/src/pks/transport/transport_ats.hh @@ -319,6 +319,7 @@ class Transport_ATS : public PK_Physical_Default { bool is_source_term_; Key cv_key_; Key molec_diff_key_; + Key primary_free_ion_concentration_key_; // workspace Key solid_residue_mass_key_; // residue -- mass that was left behind by water diff --git a/src/pks/transport/transport_ats_pk.cc b/src/pks/transport/transport_ats_pk.cc index 2d41cf6f1..68ece8fe1 100644 --- a/src/pks/transport/transport_ats_pk.cc +++ b/src/pks/transport/transport_ats_pk.cc @@ -150,6 +150,10 @@ Transport_ATS::parseParameterList() // needed by geochemical bcs molar_dens_key_ = Keys::readKey(*plist_, domain_, "molar density liquid", "molar_density_liquid"); + + primary_free_ion_concentration_key_ = + Keys::readKey(*plist_, domain_, "primary free ion concentration", + "primary_free_ion_concentration"); } // dispersion coefficient tensor @@ -581,12 +585,9 @@ Transport_ATS::SetupPhysicalEvaluators_() S_->GetRecordSetW(key_).set_subfieldnames(component_names_); if (chem_engine_ != Teuchos::null) { - Key primary_free_ion_concentration_key_ = - Keys::readKey(*plist_, domain_, "primary free ion concentration", "primary_free_ion_concentration"); - - S_->Require(primary_free_ion_concentration_key_, tag_next_, passwd_) + requireEvaluatorAtNext(primary_free_ion_concentration_key_, tag_next_, *S_) .SetMesh(mesh_) - ->SetGhosted(true) + ->SetGhosted(false) ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, num_components_); } From 3bc7b801c63bd80f0619423eebc701d1fe1187b1 Mon Sep 17 00:00:00 2001 From: gaobhub Date: Sun, 14 Jun 2026 05:38:49 -0400 Subject: [PATCH 4/6] clean code, use primary_free_ion in subsurface only --- .../generic_evaluators/MultiplicativeEvaluator.cc | 3 --- src/pks/transport/transport_ats_pk.cc | 10 ++++++---- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/src/constitutive_relations/generic_evaluators/MultiplicativeEvaluator.cc b/src/constitutive_relations/generic_evaluators/MultiplicativeEvaluator.cc index f58e66ade..6952f074c 100644 --- a/src/constitutive_relations/generic_evaluators/MultiplicativeEvaluator.cc +++ b/src/constitutive_relations/generic_evaluators/MultiplicativeEvaluator.cc @@ -62,9 +62,6 @@ MultiplicativeEvaluator::Evaluate_(const State& S, const std::vectorViewComponent(lcv_name, false)); int i = 0; for (const auto& key_tag : dependencies_) { - std::cout << "DEBUG MultiplicativeEvaluator::Evaluate_ for " - << my_keys_[0].first << std::endl; - const auto& dep_v = *(*S.Get(key_tag.first, key_tag.second).ViewComponent(lcv_name, false))( dofs_[i]); diff --git a/src/pks/transport/transport_ats_pk.cc b/src/pks/transport/transport_ats_pk.cc index 68ece8fe1..0cdc83665 100644 --- a/src/pks/transport/transport_ats_pk.cc +++ b/src/pks/transport/transport_ats_pk.cc @@ -151,9 +151,11 @@ Transport_ATS::parseParameterList() molar_dens_key_ = Keys::readKey(*plist_, domain_, "molar density liquid", "molar_density_liquid"); - primary_free_ion_concentration_key_ = - Keys::readKey(*plist_, domain_, "primary free ion concentration", - "primary_free_ion_concentration"); + if (domain_ == "domain") { + primary_free_ion_concentration_key_ = + Keys::readKey(*plist_, domain_, "primary free ion concentration", + "primary_free_ion_concentration"); + } } // dispersion coefficient tensor @@ -584,7 +586,7 @@ Transport_ATS::SetupPhysicalEvaluators_() // ->AddComponent("face", AmanziMesh::Entity_kind::FACE, num_components_); S_->GetRecordSetW(key_).set_subfieldnames(component_names_); - if (chem_engine_ != Teuchos::null) { + if (chem_engine_ != Teuchos::null && domain_ == "domain") { requireEvaluatorAtNext(primary_free_ion_concentration_key_, tag_next_, *S_) .SetMesh(mesh_) ->SetGhosted(false) From e7934fd1ba932a7eeff81c583a2df38cb05e6c52 Mon Sep 17 00:00:00 2001 From: gaobhub Date: Sun, 14 Jun 2026 16:33:38 -0400 Subject: [PATCH 5/6] clean code, primary_free_ion has been required and managed in Alquimia --- src/pks/transport/transport_ats.hh | 1 - src/pks/transport/transport_ats_pk.cc | 13 ------------- 2 files changed, 14 deletions(-) diff --git a/src/pks/transport/transport_ats.hh b/src/pks/transport/transport_ats.hh index e3160795e..78aef4d41 100644 --- a/src/pks/transport/transport_ats.hh +++ b/src/pks/transport/transport_ats.hh @@ -319,7 +319,6 @@ class Transport_ATS : public PK_Physical_Default { bool is_source_term_; Key cv_key_; Key molec_diff_key_; - Key primary_free_ion_concentration_key_; // workspace Key solid_residue_mass_key_; // residue -- mass that was left behind by water diff --git a/src/pks/transport/transport_ats_pk.cc b/src/pks/transport/transport_ats_pk.cc index 0cdc83665..618ec8377 100644 --- a/src/pks/transport/transport_ats_pk.cc +++ b/src/pks/transport/transport_ats_pk.cc @@ -150,13 +150,7 @@ Transport_ATS::parseParameterList() // needed by geochemical bcs molar_dens_key_ = Keys::readKey(*plist_, domain_, "molar density liquid", "molar_density_liquid"); - - if (domain_ == "domain") { - primary_free_ion_concentration_key_ = - Keys::readKey(*plist_, domain_, "primary free ion concentration", - "primary_free_ion_concentration"); } - } // dispersion coefficient tensor dispersion_tensor_key_ = @@ -586,13 +580,6 @@ Transport_ATS::SetupPhysicalEvaluators_() // ->AddComponent("face", AmanziMesh::Entity_kind::FACE, num_components_); S_->GetRecordSetW(key_).set_subfieldnames(component_names_); - if (chem_engine_ != Teuchos::null && domain_ == "domain") { - requireEvaluatorAtNext(primary_free_ion_concentration_key_, tag_next_, *S_) - .SetMesh(mesh_) - ->SetGhosted(false) - ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, num_components_); - } - // -- water flux requireEvaluatorAtNext(water_flux_key_, tag_next_, *S_) .SetMesh(mesh_) From d36428b38f6ba1f1b4917be2e0a262f558f59123 Mon Sep 17 00:00:00 2001 From: gaobhub Date: Sun, 14 Jun 2026 16:34:46 -0400 Subject: [PATCH 6/6] clean code, primary_free_ion has been required and managed in Alquimia --- src/pks/transport/transport_ats_pk.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pks/transport/transport_ats_pk.cc b/src/pks/transport/transport_ats_pk.cc index 618ec8377..fdb27fd34 100644 --- a/src/pks/transport/transport_ats_pk.cc +++ b/src/pks/transport/transport_ats_pk.cc @@ -150,7 +150,7 @@ Transport_ATS::parseParameterList() // needed by geochemical bcs molar_dens_key_ = Keys::readKey(*plist_, domain_, "molar density liquid", "molar_density_liquid"); - } + } // dispersion coefficient tensor dispersion_tensor_key_ =