From fc5fa4f4bfec011ddfc001110ec7ca062a61d8e0 Mon Sep 17 00:00:00 2001 From: Daniil Svyatsky Date: Tue, 16 Sep 2025 22:53:38 -0600 Subject: [PATCH 1/6] Setting Dirichlet BC for diffusion solver for each component --- src/pks/transport/transport_ats_pk.cc | 58 ++++++++++++++++----------- 1 file changed, 34 insertions(+), 24 deletions(-) diff --git a/src/pks/transport/transport_ats_pk.cc b/src/pks/transport/transport_ats_pk.cc index 057bd2d09..c5642e26d 100644 --- a/src/pks/transport/transport_ats_pk.cc +++ b/src/pks/transport/transport_ats_pk.cc @@ -286,26 +286,6 @@ Transport_ATS::SetupTransport_() S_->RequireEvaluator(dispersion_tensor_key_, tag_next_); } - // operator and boundary conditions for diffusion/dispersion solve - if (has_dispersion_ || has_diffusion_) { - // default boundary conditions (none inside domain and Neumann on its boundary) - diff_bcs_ = Teuchos::rcp( - new Operators::BCs(mesh_, AmanziMesh::Entity_kind::FACE, WhetStone::DOF_Type::SCALAR)); - PopulateBoundaryData_(-1, *diff_bcs_); - - // diffusion operator - Operators::PDE_DiffusionFactory opfactory; - Teuchos::ParameterList& op_list = plist_->sublist("diffusion"); - diff_op_ = opfactory.Create(op_list, mesh_, diff_bcs_); - diff_global_op_ = diff_op_->global_operator(); - diff_acc_op_ = - Teuchos::rcp(new Operators::PDE_Accumulation(AmanziMesh::Entity_kind::CELL, diff_global_op_)); - - // diffusion workspace - const CompositeVectorSpace& cvs = diff_global_op_->DomainMap(); - diff_sol_ = Teuchos::rcp(new CompositeVector(cvs)); - } - // source term setup // -------------------------------------------------------------------------------- if (plist_->isSublist("source terms")) { @@ -449,6 +429,7 @@ Transport_ATS::SetupTransport_() auto bcs_list = Teuchos::sublist(plist_, "boundary conditions"); auto conc_bcs_list = Teuchos::sublist(bcs_list, "mole fraction"); + int m = 0; for (const auto& it : *conc_bcs_list) { std::string name = it.first; if (conc_bcs_list->isSublist(name)) { @@ -503,13 +484,40 @@ Transport_ATS::SetupTransport_() // set the component indicies for (const auto& n : bc->tcc_names()) { - bc->tcc_index().push_back(FindComponentNumber_(n)); + std::cout<<"n "<tcc_index().push_back(FindComponentNumber_(n)); } bcs_.push_back(bc); + for (auto it = bc->begin(); it != bc->end(); ++it) { + std::cout << it->first <<"\n"; + } } } } + // operator and boundary conditions for diffusion/dispersion solve + if (has_dispersion_ || has_diffusion_) { + // default boundary conditions (none inside domain and Neumann on its boundary) + diff_bcs_ = Teuchos::rcp( + new Operators::BCs(mesh_, AmanziMesh::Entity_kind::FACE, WhetStone::DOF_Type::SCALAR)); + + PopulateBoundaryData_(-1, *diff_bcs_); + + // diffusion operator + Operators::PDE_DiffusionFactory opfactory; + Teuchos::ParameterList& op_list = plist_->sublist("diffusion"); + diff_op_ = opfactory.Create(op_list, mesh_, diff_bcs_); + diff_global_op_ = diff_op_->global_operator(); + diff_acc_op_ = + Teuchos::rcp(new Operators::PDE_Accumulation(AmanziMesh::Entity_kind::CELL, diff_global_op_)); + + // diffusion workspace + const CompositeVectorSpace& cvs = diff_global_op_->DomainMap(); + diff_sol_ = Teuchos::rcp(new CompositeVector(cvs)); + } + + + #ifdef ALQUIMIA_ENABLED // -- try geochemical conditions auto geochem_list = Teuchos::sublist(bcs_list, "geochemical"); @@ -950,7 +958,7 @@ Transport_ATS ::AdvanceDispersionDiffusion_(double t_old, double t_new) { if (!has_diffusion_ && !has_dispersion_) return; double dt = t_new - t_old; - + Epetra_MultiVector& tcc_new = *S_->GetW(key_, tag_next_, passwd_).ViewComponent("cell", false); @@ -983,6 +991,9 @@ Transport_ATS ::AdvanceDispersionDiffusion_(double t_old, double t_new) double md_old(0.0); for (int i = 0; i != num_aqueous_; ++i) { + + // Set Dirichlet_BC for + PopulateBoundaryData_(i, *diff_bcs_); // add molecular diffusion to the dispersion tensor bool changed_tensor(false); if (has_diffusion_) { @@ -1346,8 +1357,7 @@ Transport_ATS::PopulateBoundaryData_(int component, Operators::BCs& bc) if (component >= 0) { for (int m = 0; m < bcs_.size(); m++) { std::vector& tcc_index = bcs_[m]->tcc_index(); - int ncomp = tcc_index.size(); - + int ncomp = tcc_index.size(); for (auto it = bcs_[m]->begin(); it != bcs_[m]->end(); ++it) { int f = it->first; std::vector& values = it->second; From bb3075a92e055c29b924e09bace03380a439df80 Mon Sep 17 00:00:00 2001 From: Daniil Svyatsky Date: Wed, 26 Nov 2025 13:13:01 -0700 Subject: [PATCH 2/6] Remove debug output --- src/pks/transport/transport_ats_pk.cc | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/pks/transport/transport_ats_pk.cc b/src/pks/transport/transport_ats_pk.cc index c5642e26d..df0183c06 100644 --- a/src/pks/transport/transport_ats_pk.cc +++ b/src/pks/transport/transport_ats_pk.cc @@ -484,13 +484,9 @@ Transport_ATS::SetupTransport_() // set the component indicies for (const auto& n : bc->tcc_names()) { - std::cout<<"n "<tcc_index().push_back(FindComponentNumber_(n)); } bcs_.push_back(bc); - for (auto it = bc->begin(); it != bc->end(); ++it) { - std::cout << it->first <<"\n"; - } } } } From 733441491ac3c5ae2ee109e2c8f1f211cee4be31 Mon Sep 17 00:00:00 2001 From: Daniil Svyatsky Date: Fri, 12 Dec 2025 17:10:01 -0700 Subject: [PATCH 3/6] Add as an option enforcement of Dirichlet BC in transport diffusion --- src/pks/transport/transport_ats.hh | 2 +- src/pks/transport/transport_ats_pk.cc | 6 ++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/pks/transport/transport_ats.hh b/src/pks/transport/transport_ats.hh index f42c2380f..7e5218ff0 100644 --- a/src/pks/transport/transport_ats.hh +++ b/src/pks/transport/transport_ats.hh @@ -372,7 +372,7 @@ class Transport_ATS : public PK_Physical_Default { std::vector> bcs_; // operators for dispersion/diffusion - bool has_diffusion_, has_dispersion_; + bool has_diffusion_, has_dispersion_, enforce_bc_; Teuchos::RCP D_; // workspace, disp + diff Teuchos::RCP diff_bcs_; Teuchos::RCP diff_global_op_; diff --git a/src/pks/transport/transport_ats_pk.cc b/src/pks/transport/transport_ats_pk.cc index df0183c06..4b3a6ccaa 100644 --- a/src/pks/transport/transport_ats_pk.cc +++ b/src/pks/transport/transport_ats_pk.cc @@ -63,7 +63,8 @@ Transport_ATS::Transport_ATS(Teuchos::ParameterList& pk_tree, dt_stable_(-1.0), dt_max_(-1.0), has_diffusion_(false), - has_dispersion_(false) + has_dispersion_(false), + enforce_bc_(false) { // initialize io units_.Init(global_plist->sublist("units")); @@ -502,6 +503,7 @@ Transport_ATS::SetupTransport_() // diffusion operator Operators::PDE_DiffusionFactory opfactory; Teuchos::ParameterList& op_list = plist_->sublist("diffusion"); + enforce_bc_ = op_list.get("enforce boundary conditions", false); diff_op_ = opfactory.Create(op_list, mesh_, diff_bcs_); diff_global_op_ = diff_op_->global_operator(); diff_acc_op_ = @@ -989,7 +991,7 @@ Transport_ATS ::AdvanceDispersionDiffusion_(double t_old, double t_new) for (int i = 0; i != num_aqueous_; ++i) { // Set Dirichlet_BC for - PopulateBoundaryData_(i, *diff_bcs_); + if (enforce_bc_) PopulateBoundaryData_(i, *diff_bcs_); // add molecular diffusion to the dispersion tensor bool changed_tensor(false); if (has_diffusion_) { From 56551184fc3ed9a6100725a6a95df292bb942975 Mon Sep 17 00:00:00 2001 From: Daniil Svyatsky Date: Mon, 5 Jan 2026 18:50:06 -0700 Subject: [PATCH 4/6] For EOS evaluator modified EnsureCompatibility_ToDeps_ to support multicomponent dependencies for mass and molar densities/ --- src/constitutive_relations/eos/eos_evaluator.cc | 13 +++++++++++++ src/constitutive_relations/eos/eos_evaluator.hh | 2 ++ 2 files changed, 15 insertions(+) diff --git a/src/constitutive_relations/eos/eos_evaluator.cc b/src/constitutive_relations/eos/eos_evaluator.cc index 96e00de78..f7bc7df51 100644 --- a/src/constitutive_relations/eos/eos_evaluator.cc +++ b/src/constitutive_relations/eos/eos_evaluator.cc @@ -348,6 +348,19 @@ EOSEvaluator::EvaluatePartialDerivative_(const State& S, } } +void +EOSEvaluator::EnsureCompatibility_ToDeps_(State& S){ + + auto akeytag = my_keys_.front(); + const auto& fac = + S.Require(akeytag.first, akeytag.second); + + for (const auto& dep : dependencies_) { + auto& dep_fac = S.Require(dep.first, dep.second); + dep_fac.UpdateComponentsSameNumDofs(fac); + } + +} } // namespace Relations } // namespace Amanzi diff --git a/src/constitutive_relations/eos/eos_evaluator.hh b/src/constitutive_relations/eos/eos_evaluator.hh index a6b640a41..d6e4d68d2 100644 --- a/src/constitutive_relations/eos/eos_evaluator.hh +++ b/src/constitutive_relations/eos/eos_evaluator.hh @@ -61,6 +61,8 @@ class EOSEvaluator : public EvaluatorSecondaryMonotypeCV { } + virtual void EnsureCompatibility_ToDeps_(State& S) override; + // Required methods from EvaluatorSecondaryMonotypeCV virtual void Evaluate_(const State& S, const std::vector& results) override; virtual void EvaluatePartialDerivative_(const State& S, From 7d2f5f8a43663fd43e3f0f8cb6c36165e3af2464 Mon Sep 17 00:00:00 2001 From: Daniil Svyatsky Date: Thu, 8 Jan 2026 18:08:20 -0700 Subject: [PATCH 5/6] Initialization of a mesh during EnsureCompatibility_ToDeps_ --- src/constitutive_relations/eos/eos_evaluator.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/constitutive_relations/eos/eos_evaluator.cc b/src/constitutive_relations/eos/eos_evaluator.cc index f7bc7df51..064b545f2 100644 --- a/src/constitutive_relations/eos/eos_evaluator.cc +++ b/src/constitutive_relations/eos/eos_evaluator.cc @@ -357,7 +357,7 @@ EOSEvaluator::EnsureCompatibility_ToDeps_(State& S){ for (const auto& dep : dependencies_) { auto& dep_fac = S.Require(dep.first, dep.second); - dep_fac.UpdateComponentsSameNumDofs(fac); + dep_fac.UpdateSameNumDofs(fac); } } From e6194cb0d61730587971ed158bd8abff31d628bd Mon Sep 17 00:00:00 2001 From: Daniil Svyatsky Date: Mon, 19 Jan 2026 22:38:38 -0700 Subject: [PATCH 6/6] Fix the bug related to determenation of the presense of dispertion and/or diffusion --- src/pks/transport/transport_ats_pk.cc | 47 +++++++++++++-------------- 1 file changed, 23 insertions(+), 24 deletions(-) diff --git a/src/pks/transport/transport_ats_pk.cc b/src/pks/transport/transport_ats_pk.cc index 4b3a6ccaa..6e0a6d99f 100644 --- a/src/pks/transport/transport_ats_pk.cc +++ b/src/pks/transport/transport_ats_pk.cc @@ -287,6 +287,29 @@ Transport_ATS::SetupTransport_() S_->RequireEvaluator(dispersion_tensor_key_, tag_next_); } + // operator and boundary conditions for diffusion/dispersion solve + if (has_dispersion_ || has_diffusion_) { + // default boundary conditions (none inside domain and Neumann on its boundary) + diff_bcs_ = Teuchos::rcp( + new Operators::BCs(mesh_, AmanziMesh::Entity_kind::FACE, WhetStone::DOF_Type::SCALAR)); + + PopulateBoundaryData_(-1, *diff_bcs_); + + // diffusion operator + Operators::PDE_DiffusionFactory opfactory; + Teuchos::ParameterList& op_list = plist_->sublist("diffusion"); + enforce_bc_ = op_list.get("enforce boundary conditions", false); + diff_op_ = opfactory.Create(op_list, mesh_, diff_bcs_); + diff_global_op_ = diff_op_->global_operator(); + diff_acc_op_ = + Teuchos::rcp(new Operators::PDE_Accumulation(AmanziMesh::Entity_kind::CELL, diff_global_op_)); + + // diffusion workspace + const CompositeVectorSpace& cvs = diff_global_op_->DomainMap(); + diff_sol_ = Teuchos::rcp(new CompositeVector(cvs)); + } + + // source term setup // -------------------------------------------------------------------------------- if (plist_->isSublist("source terms")) { @@ -491,30 +514,6 @@ Transport_ATS::SetupTransport_() } } } - - // operator and boundary conditions for diffusion/dispersion solve - if (has_dispersion_ || has_diffusion_) { - // default boundary conditions (none inside domain and Neumann on its boundary) - diff_bcs_ = Teuchos::rcp( - new Operators::BCs(mesh_, AmanziMesh::Entity_kind::FACE, WhetStone::DOF_Type::SCALAR)); - - PopulateBoundaryData_(-1, *diff_bcs_); - - // diffusion operator - Operators::PDE_DiffusionFactory opfactory; - Teuchos::ParameterList& op_list = plist_->sublist("diffusion"); - enforce_bc_ = op_list.get("enforce boundary conditions", false); - diff_op_ = opfactory.Create(op_list, mesh_, diff_bcs_); - diff_global_op_ = diff_op_->global_operator(); - diff_acc_op_ = - Teuchos::rcp(new Operators::PDE_Accumulation(AmanziMesh::Entity_kind::CELL, diff_global_op_)); - - // diffusion workspace - const CompositeVectorSpace& cvs = diff_global_op_->DomainMap(); - diff_sol_ = Teuchos::rcp(new CompositeVector(cvs)); - } - - #ifdef ALQUIMIA_ENABLED // -- try geochemical conditions