WIP on chemistry

src/executables/coordinator.cc

@@ -303,16 +303,15 @@ Coordinator::initialize()
   int size = comm_->NumProc();
   int rank = comm_->MyPID();
 
-  S_->set_time(Amanzi::Tags::CURRENT, t0_);
-  S_->set_time(Amanzi::Tags::NEXT, t0_);
-  S_->set_cycle(cycle0_);
+  for (const auto& tag : S_->GetRecordSet("time"))
+    S_->set_time(tag.first, t0_);
 
   // Restart from checkpoint part 1:
   //  - get the time prior to initializing anything else
   if (restart_) {
     double t_restart = Amanzi::ReadCheckpointInitialTime(comm_, restart_filename_);
-    S_->set_time(Amanzi::Tags::CURRENT, t_restart);
-    S_->set_time(Amanzi::Tags::NEXT, t_restart);
+    for (const auto& tag : S_->GetRecordSet("time"))
+      S_->set_time(tag.first, t_restart);
     t0_ = t_restart;
   }
 
@@ -344,8 +343,8 @@ Coordinator::initialize()
     t0_ = S_->get_time();
     cycle0_ = S_->get_cycle();
 
-    S_->set_time(Amanzi::Tags::CURRENT, t0_);
-    S_->set_time(Amanzi::Tags::NEXT, t0_);
+    for (const auto& tag : S_->GetRecordSet("time"))
+      S_->set_time(tag.first, t0_);
 
     for (Amanzi::State::mesh_iterator mesh = S_->mesh_begin(); mesh != S_->mesh_end(); ++mesh) {
       if (S_->IsDeformableMesh(mesh->first) && !S_->IsAliasedMesh(mesh->first))

src/pks/mpc/mpc_coupled_reactivetransport.cc

@@ -63,12 +63,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_dynamic_cast<AmanziChemistry::Alquimia_PK>(chemistry_pk_));
+  transport_pk_surf_->setChemEngine(Teuchos::rcp_dynamic_cast<AmanziChemistry::Alquimia_PK>(chemistry_pk_surf_));
 #endif
 
   coupled_transport_pk_->parseParameterList();
@@ -151,32 +147,27 @@ MPCCoupledReactiveTransport::Initialize()
   // after the density of water can be evaluated.  This could be problematic
   // for, e.g., salinity intrusion problems where water density is a function
   // of concentration itself, but should work for all other problems?
-  Teuchos::RCP<Epetra_MultiVector> tcc_surf =
-    S_->GetW<CompositeVector>(tcc_surf_key_, tag_next_, name_).ViewComponent("cell", true);
+  Epetra_MultiVector& tcc_surf =
+    *S_->GetW<CompositeVector>(tcc_surf_key_, tag_next_, name_).ViewComponent("cell", true);
   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);
+  const Epetra_MultiVector& mol_dens_surf =
+    *S_->Get<CompositeVector>(mol_dens_surf_key_, tag_next_).ViewComponent("cell", true);
 
-  Teuchos::RCP<Epetra_MultiVector> tcc =
-    S_->GetW<CompositeVector>(tcc_key_, tag_next_, name_).ViewComponent("cell", true);
+  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);
+  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, num_aqueous, tcc_surf, tcc_surf);
+  convertConcentrationToAmanzi(mol_dens, num_aqueous, 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, num_aqueous, tcc_surf, tcc_surf);
+  convertConcentrationToATS(mol_dens, num_aqueous, tcc, tcc);
 
   transport_pk_surf_->Initialize();
   transport_pk_->Initialize();

src/pks/mpc/mpc_reactivetransport.cc

@@ -52,9 +52,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_dynamic_cast<AmanziChemistry::Alquimia_PK>(chemistry_pk_));
 #endif
 
   // now transport can be parse with knowledge of names
@@ -103,18 +101,16 @@ MPCReactiveTransport::Initialize()
   // after the density of water can be evaluated.  This could be problematic
   // for, e.g., salinity intrusion problems where water density is a function
   // of concentration itself, but should work for all other problems?
-  Teuchos::RCP<Epetra_MultiVector> tcc =
-    S_->GetW<CompositeVector>(tcc_key_, tag_next_, name_).ViewComponent("cell", true);
+  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);
+  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, num_aqueous, tcc, tcc);
   chemistry_pk_->Initialize();
-  //*tcc = *chemistry_pk_->aqueous_components();
-  convertConcentrationToATS(*mol_dens, num_aqueous, *tcc, *tcc);
+  convertConcentrationToATS(mol_dens, num_aqueous, tcc, tcc);
 
   transport_pk_->Initialize();
 }

src/pks/pk_helpers.cc

@@ -267,22 +267,20 @@ advanceChemistry(Teuchos::RCP<AmanziChemistry::Chemistry_PK> chem_pk,
                  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, num_aqueous, tcc, tcc);
 
   {
     auto monitor = Teuchos::rcp(new Teuchos::TimeMonitor(timer));
     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, num_aqueous, tcc, tcc);
   return fail;
 }
 

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(const Teuchos::RCP<AmanziChemistry::Alquimia_PK>& chem_pk);
 #endif
 
   virtual void Initialize() override;

src/pks/transport/transport_ats_pk.cc

@@ -182,13 +182,12 @@ 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(const Teuchos::RCP<AmanziChemistry::Alquimia_PK>& chem_pk)
 {
   chem_pk_ = chem_pk;
-  chem_engine_ = chem_engine;
+  chem_engine_ = chem_pk->chem_engine();
 
-  if (chem_engine_ != Teuchos::null) {
+  if (chem_pk_ != Teuchos::null) {
     // Retrieve the component names (primary and secondary) from the chemistry
     // engine.
     chem_engine_->GetPrimarySpeciesNames(component_names_);