chemistry changes -- ATS builds

src/pks/chem_pk_helpers.cc

@@ -17,56 +17,41 @@ namespace Amanzi {
 // -----------------------------------------------------------------------------
 void
 convertConcentrationToAmanzi(const Epetra_MultiVector& mol_dens,
-                             int num_aqueous,
                              const Epetra_MultiVector& tcc_ats,
                              Epetra_MultiVector& tcc_amanzi)
 {
   // convert from mole fraction [mol C / mol H20] to [mol C / L]
-  for (int k = 0; k != num_aqueous; ++k) {
-    for (int c = 0; c != tcc_ats.MyLength(); ++c) {
-      // 1.e-3 converts L to m^3
-      tcc_amanzi[k][c] = tcc_ats[k][c] * mol_dens[0][c] * 1.e-3;
-    }
-  }
+  tcc_amanzi.Multiply(1.e-3, mol_dens, tcc_ats, 0.);
 }
 
 void
 convertConcentrationToATS(const Epetra_MultiVector& mol_dens,
-                          int num_aqueous,
                           const Epetra_MultiVector& tcc_amanzi,
                           Epetra_MultiVector& tcc_ats)
 {
   // convert from [mol C / L] to mol fraction [mol C / mol H20]
-  for (int k = 0; k != num_aqueous; ++k) {
-    for (int c = 0; c != tcc_amanzi.MyLength(); ++c) {
-      tcc_ats[k][c] = tcc_amanzi[k][c] / (mol_dens[0][c] * 1.e-3);
-    }
-  }
+  tcc_ats.ReciprocalMultiply(1.e3, tcc_amanzi, mol_dens, 0.);
 }
 
 
 bool
-advanceChemistry(Teuchos::RCP<AmanziChemistry::Chemistry_PK> chem_pk,
+advanceChemistry(AmanziChemistry::Chemistry_PK& chem_pk,
                  double t_old,
                  double t_new,
                  bool reinit,
                  const Epetra_MultiVector& mol_dens,
-                 Teuchos::RCP<Epetra_MultiVector> tcc,
+                 Epetra_MultiVector& tcc,
                  Teuchos::Time& timer)
 {
   bool fail = false;
-  int num_aqueous = chem_pk->num_aqueous_components();
-  convertConcentrationToAmanzi(mol_dens, num_aqueous, *tcc, *tcc);
-  chem_pk->set_aqueous_components(tcc);
-
+  convertConcentrationToAmanzi(mol_dens, tcc, tcc);
   {
     auto monitor = Teuchos::rcp(new Teuchos::TimeMonitor(timer));
-    fail = chem_pk->AdvanceStep(t_old, t_new, reinit);
+    fail = chem_pk.AdvanceStep(t_old, t_new, reinit);
   }
   if (fail) return fail;
 
-  *tcc = *chem_pk->aqueous_components();
-  convertConcentrationToATS(mol_dens, num_aqueous, *tcc, *tcc);
+  convertConcentrationToATS(mol_dens, tcc, tcc);
   return fail;
 }
 

src/pks/chem_pk_helpers.hh

@@ -11,7 +11,7 @@
 #pragma once
 
 #include "Teuchos_TimeMonitor.hpp"
-#include "Epetra_MultiVector.hpp"
+#include "Epetra_MultiVector.h"
 
 namespace Amanzi {
 
@@ -24,23 +24,21 @@ class Chemistry_PK;
 // -----------------------------------------------------------------------------
 void
 convertConcentrationToAmanzi(const Epetra_MultiVector& mol_den,
-                             int num_aqueous,
                              const Epetra_MultiVector& tcc_ats,
                              Epetra_MultiVector& tcc_amanzi);
 
 void
 convertConcentrationToATS(const Epetra_MultiVector& mol_den,
-                          int num_aqueous,
                           const Epetra_MultiVector& tcc_ats,
                           Epetra_MultiVector& tcc_amanzi);
 
 bool
-advanceChemistry(Teuchos::RCP<AmanziChemistry::Chemistry_PK> chem_pk,
+advanceChemistry(AmanziChemistry::Chemistry_PK& chem_pk,
                  double t_old,
                  double t_new,
                  bool reinit,
                  const Epetra_MultiVector& mol_dens,
-                 Teuchos::RCP<Epetra_MultiVector> tcc,
+                 Epetra_MultiVector& tcc,
                  Teuchos::Time& timer);
 
 } // namespace Amanzi

src/pks/energy/energy_base_ti.cc

@@ -10,9 +10,11 @@
 #include "Debugger.hh"
 #include "BoundaryFunction.hh"
 #include "Evaluator.hh"
-#include "energy_base.hh"
+#include "pk_helpers.hh"
 #include "Op.hh"
 
+#include "energy_base.hh"
+
 namespace Amanzi {
 namespace Energy {
 

src/pks/flow/snow_distribution_ti.cc

@@ -8,6 +8,7 @@
 */
 
 #include "Op.hh"
+#include "pk_helpers.hh"
 #include "snow_distribution.hh"
 
 #define DEBUG_FLAG 1

src/pks/mpc/mpc_coupled_reactivetransport.cc

@@ -14,6 +14,7 @@
 */
 
 #include "mpc_coupled_reactivetransport.hh"
+#include "chem_pk_helpers.hh"
 #include "pk_helpers.hh"
 
 namespace Amanzi {
@@ -63,12 +64,8 @@ MPCCoupledReactiveTransport::parseParameterList()
 
   // communicate chemistry engine to transport.
 #ifdef ALQUIMIA_ENABLED
-  transport_pk_->setChemEngine(
-    Teuchos::rcp_static_cast<AmanziChemistry::Alquimia_PK>(chemistry_pk_),
-    chemistry_pk_->chem_engine());
-  transport_pk_surf_->setChemEngine(
-    Teuchos::rcp_static_cast<AmanziChemistry::Alquimia_PK>(chemistry_pk_surf_),
-    chemistry_pk_surf_->chem_engine());
+  transport_pk_->setChemEngine(Teuchos::rcp_static_cast<AmanziChemistry::Alquimia_PK>(chemistry_pk_));
+  transport_pk_surf_->setChemEngine(Teuchos::rcp_static_cast<AmanziChemistry::Alquimia_PK>(chemistry_pk_surf_));
 #endif
 
   coupled_transport_pk_->parseParameterList();
@@ -85,31 +82,6 @@ MPCCoupledReactiveTransport::Setup()
   // must Setup transport first to get alias for saturation, etc set up correctly
   coupled_transport_pk_->Setup();
   coupled_chemistry_pk_->Setup();
-
-  S_->Require<CompositeVector, CompositeVectorSpace>(tcc_key_, tag_next_, name_)
-    .SetMesh(S_->GetMesh(domain_))
-    ->SetGhosted()
-    ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, chemistry_pk_->num_aqueous_components());
-  S_->RequireEvaluator(tcc_key_, tag_next_);
-
-  S_->Require<CompositeVector, CompositeVectorSpace>(tcc_surf_key_, tag_next_, name_)
-    .SetMesh(S_->GetMesh(domain_surf_))
-    ->SetGhosted()
-    ->AddComponent(
-      "cell", AmanziMesh::Entity_kind::CELL, chemistry_pk_surf_->num_aqueous_components());
-  S_->RequireEvaluator(tcc_surf_key_, tag_next_);
-
-  S_->Require<CompositeVector, CompositeVectorSpace>(mol_dens_key_, tag_next_)
-    .SetMesh(S_->GetMesh(domain_))
-    ->SetGhosted()
-    ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1);
-  S_->RequireEvaluator(mol_dens_key_, tag_next_);
-
-  S_->Require<CompositeVector, CompositeVectorSpace>(mol_dens_surf_key_, tag_next_)
-    .SetMesh(S_->GetMesh(domain_surf_))
-    ->SetGhosted()
-    ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1);
-  S_->RequireEvaluator(mol_dens_surf_key_, tag_next_);
 }
 
 
@@ -163,20 +135,14 @@ MPCCoupledReactiveTransport::Initialize()
   Teuchos::RCP<const Epetra_MultiVector> mol_dens =
     S_->Get<CompositeVector>(mol_dens_key_, tag_next_).ViewComponent("cell", true);
 
-  int num_aqueous = chemistry_pk_surf_->num_aqueous_components();
-  convertConcentrationToAmanzi(*mol_dens_surf, num_aqueous, *tcc_surf, *tcc_surf);
-  convertConcentrationToAmanzi(*mol_dens, num_aqueous, *tcc, *tcc);
+  convertConcentrationToAmanzi(*mol_dens_surf, *tcc_surf, *tcc_surf);
+  convertConcentrationToAmanzi(*mol_dens, *tcc, *tcc);
 
-  chemistry_pk_surf_->set_aqueous_components(tcc_surf);
   chemistry_pk_surf_->Initialize();
-  //*tcc_surf = *chemistry_pk_surf_->aqueous_components();
-
-  chemistry_pk_->set_aqueous_components(tcc);
   chemistry_pk_->Initialize();
-  //*tcc = *chemistry_pk_->aqueous_components();
 
-  convertConcentrationToATS(*mol_dens_surf, num_aqueous, *tcc_surf, *tcc_surf);
-  convertConcentrationToATS(*mol_dens, num_aqueous, *tcc, *tcc);
+  convertConcentrationToATS(*mol_dens_surf, *tcc_surf, *tcc_surf);
+  convertConcentrationToATS(*mol_dens, *tcc, *tcc);
 
   transport_pk_surf_->Initialize();
   transport_pk_->Initialize();
@@ -223,8 +189,8 @@ MPCCoupledReactiveTransport::AdvanceStep(double t_old, double t_new, bool reinit
   S_->GetEvaluator(mol_dens_surf_key_, tag_next_).Update(*S_, name_);
   Teuchos::RCP<const Epetra_MultiVector> mol_dens_surf =
     S_->Get<CompositeVector>(mol_dens_surf_key_, tag_next_).ViewComponent("cell", true);
-  fail = advanceChemistry(
-    chemistry_pk_surf_, t_old, t_new, reinit, *mol_dens_surf, tcc_surf, *alquimia_surf_timer_);
+
+  fail = advanceChemistry(*chemistry_pk_surf_, t_old, t_new, reinit, *mol_dens_surf, *tcc_surf, *alquimia_surf_timer_);
   changedEvaluatorPrimary(tcc_surf_key_, tag_next_, *S_);
   if (fail) {
     if (vo_->os_OK(Teuchos::VERB_MEDIUM))
@@ -242,7 +208,7 @@ MPCCoupledReactiveTransport::AdvanceStep(double t_old, double t_new, bool reinit
   Teuchos::RCP<const Epetra_MultiVector> mol_dens =
     S_->Get<CompositeVector>(mol_dens_key_, tag_next_).ViewComponent("cell", true);
   try {
-    fail = advanceChemistry(chemistry_pk_, t_old, t_new, reinit, *mol_dens, tcc, *alquimia_timer_);
+    fail = advanceChemistry(*chemistry_pk_, t_old, t_new, reinit, *mol_dens, *tcc, *alquimia_timer_);
     changedEvaluatorPrimary(tcc_key_, tag_next_, *S_);
   } catch (const Errors::Message& chem_error) {
     fail = true;

src/pks/mpc/mpc_reactivetransport.cc

@@ -15,6 +15,7 @@
 
 #include "mpc_reactivetransport.hh"
 #include "pk_helpers.hh"
+#include "chem_pk_helpers.hh"
 
 namespace Amanzi {
 
@@ -52,9 +53,7 @@ MPCReactiveTransport::parseParameterList()
   chemistry_pk_->parseParameterList();
 
 #ifdef ALQUIMIA_ENABLED
-  transport_pk_->setChemEngine(
-    Teuchos::rcp_static_cast<AmanziChemistry::Alquimia_PK>(chemistry_pk_),
-    chemistry_pk_->chem_engine());
+  transport_pk_->setChemEngine(Teuchos::rcp_static_cast<AmanziChemistry::Alquimia_PK>(chemistry_pk_));
 #endif
 
   // now transport can be parse with knowledge of names
@@ -68,12 +67,6 @@ MPCReactiveTransport::Setup()
   transport_pk_->Setup();
   chemistry_pk_->Setup();
 
-  S_->Require<CompositeVector, CompositeVectorSpace>(tcc_key_, tag_next_, name_)
-    .SetMesh(S_->GetMesh(domain_))
-    ->SetGhosted()
-    ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, chemistry_pk_->num_aqueous_components());
-  S_->RequireEvaluator(tcc_key_, tag_next_);
-
   S_->Require<CompositeVector, CompositeVectorSpace>(mol_dens_key_, tag_next_)
     .SetMesh(S_->GetMesh(domain_))
     ->SetGhosted()
@@ -109,13 +102,9 @@ MPCReactiveTransport::Initialize()
   Teuchos::RCP<const Epetra_MultiVector> mol_dens =
     S_->Get<CompositeVector>(mol_dens_key_, tag_next_).ViewComponent("cell", true);
 
-  int num_aqueous = chemistry_pk_->num_aqueous_components();
-  convertConcentrationToAmanzi(*mol_dens, num_aqueous, *tcc, *tcc);
-  chemistry_pk_->set_aqueous_components(tcc);
+  convertConcentrationToAmanzi(*mol_dens, *tcc, *tcc);
   chemistry_pk_->Initialize();
-  //*tcc = *chemistry_pk_->aqueous_components();
-  convertConcentrationToATS(*mol_dens, num_aqueous, *tcc, *tcc);
-
+  convertConcentrationToATS(*mol_dens, *tcc, *tcc);
   transport_pk_->Initialize();
 }
 
@@ -153,15 +142,15 @@ MPCReactiveTransport::AdvanceStep(double t_old, double t_new, bool reinit)
   // Second, we do a chemistry step.
   int num_aq_componets = transport_pk_->get_num_aqueous_component();
 
-  Teuchos::RCP<Epetra_MultiVector> tcc_copy =
+  Teuchos::RCP<Epetra_MultiVector> tcc =
     S_->GetW<CompositeVector>(tcc_key_, tag_next_, name_).ViewComponent("cell", true);
 
   S_->GetEvaluator(mol_dens_key_, tag_next_).Update(*S_, name_);
   Teuchos::RCP<const Epetra_MultiVector> mol_dens =
     S_->Get<CompositeVector>(mol_dens_key_, tag_next_).ViewComponent("cell", true);
 
   fail |=
-    advanceChemistry(chemistry_pk_, t_old, t_new, reinit, *mol_dens, tcc_copy, *alquimia_timer_);
+    advanceChemistry(*chemistry_pk_, t_old, t_new, reinit, *mol_dens, *tcc, *alquimia_timer_);
   changedEvaluatorPrimary(tcc_key_, tag_next_, *S_);
   if (!fail) chem_step_succeeded_ = true;
 

src/pks/transport/transport_ats.hh

@@ -321,8 +321,7 @@ class Transport_ATS : public PK_Physical_Default {
 
   // coupling with chemistry
 #ifdef ALQUIMIA_ENABLED
-  void setChemEngine(Teuchos::RCP<AmanziChemistry::Alquimia_PK> chem_pk,
-                     Teuchos::RCP<AmanziChemistry::ChemistryEngine> chem_engine);
+  void setChemEngine(Teuchos::RCP<AmanziChemistry::Alquimia_PK> chem_pk);
 #endif
 
   virtual void Initialize() override;

src/pks/transport/transport_ats_pk.cc

@@ -182,11 +182,10 @@ Transport_ATS::readParameterMapByPhase(Teuchos::ParameterList& plist,
 ****************************************************************** */
 #ifdef ALQUIMIA_ENABLED
 void
-Transport_ATS::setChemEngine(Teuchos::RCP<AmanziChemistry::Alquimia_PK> chem_pk,
-                             Teuchos::RCP<AmanziChemistry::ChemistryEngine> chem_engine)
+Transport_ATS::setChemEngine(Teuchos::RCP<AmanziChemistry::Alquimia_PK> chem_pk)
 {
   chem_pk_ = chem_pk;
-  chem_engine_ = chem_engine;
+  chem_engine_ = chem_pk->getChemEngine();
 
   if (chem_engine_ != Teuchos::null) {
     // Retrieve the component names (primary and secondary) from the chemistry