work on coupled flow and transport

src/pks/flow/richards_physics.cc

@@ -226,6 +226,8 @@ Richards::UpdateVelocity_(const Tag& tag)
 
     for (int i = 0; i != d; ++i) velocity[i][c] = rhs[i] / nliq_c[0][c];
   }
+
+  changedEvaluatorPrimary(velocity_key_, tag, *S_);
 }
 
 

src/pks/flow/richards_pk.cc

@@ -90,7 +90,10 @@ Richards::parseParameterList()
   requireAtNext(flux_key_, tag_next_, *S_, name_);
 
   flux_dir_key_ = Keys::readKey(*plist_, domain_, "darcy flux direction", "water_flux_direction");
+
   velocity_key_ = Keys::readKey(*plist_, domain_, "darcy velocity", "darcy_velocity");
+  requireAtNext(velocity_key_, Tags::NEXT, *S_, name_);
+
   sat_key_ = Keys::readKey(*plist_, domain_, "saturation", "saturation_liquid");
   sat_gas_key_ = Keys::readKey(*plist_, domain_, "saturation gas", "saturation_gas");
   sat_ice_key_ = Keys::readKey(*plist_, domain_, "saturation ice", "saturation_ice");
@@ -451,9 +454,6 @@ Richards::SetupPhysicalEvaluators_()
   S_->RequireDerivative<CompositeVector, CompositeVectorSpace>(
     conserved_key_, tag_next_, key_, tag_next_);
 
-  //    and at the current time, where it is a copy evaluator
-  requireAtCurrent(conserved_key_, tag_current_, *S_, name_);
-
   // -- Water retention evaluators
   // -- saturation
   requireAtNext(sat_key_, tag_next_, *S_, true)

src/pks/mpc/mpc_flow_transport.cc

@@ -7,6 +7,7 @@
   Authors: Ethan Coon
 */
 
+#include "pk_helpers.hh"
 #include "mpc_flow_transport.hh"
 
 namespace Amanzi {
@@ -30,10 +31,10 @@ MPCFlowTransport::parseParameterList()
     subsurface = true;
   }
 
-  if (subsurface) {
-    auto [flow_current_tag, flow_next_tag] = tags_[0];
-    auto [transport_current_tag, transport_next_tag] = tags_[1];
+  auto [flow_current_tag, flow_next_tag] = tags_[0];
+  auto [transport_current_tag, transport_next_tag] = tags_[1];
 
+  if (subsurface) {
     if (transport_next_tag != flow_next_tag) {
       // set the flow field evaluator as the flow's NEXT tag
       //
@@ -45,27 +46,26 @@ MPCFlowTransport::parseParameterList()
       flux_list.set<std::string>("evaluator type", "alias");
       flux_list.set<std::string>("target", Keys::getKey("water_flux", flow_next_tag, true));
 
-      // set the saturation_liquid as an interpolated field
-      Teuchos::ParameterList& sat_list_current = S_->GetEvaluatorList(Keys::getKey("saturation_liquid", transport_current_tag));
-      sat_list_current.set<std::string>("evaluator type", "temporal interpolation");
-      sat_list_current.set<std::string>("current tag", flow_current_tag.get());
-      sat_list_current.set<std::string>("next tag", flow_next_tag.get());
+      // velocity for dispersivity
+      Teuchos::ParameterList& velo_list = S_->GetEvaluatorList(Keys::getKey("velocity", transport_next_tag));
+      velo_list.set<std::string>("evaluator type", "alias");
+      velo_list.set<std::string>("target", Keys::getKey("velocity", flow_next_tag, true));
 
-      Teuchos::ParameterList& sat_list_next = S_->GetEvaluatorList(Keys::getKey("saturation_liquid", transport_next_tag));
-      sat_list_next.set<std::string>("evaluator type", "temporal interpolation");
-      sat_list_next.set<std::string>("current tag", flow_current_tag.get());
-      sat_list_next.set<std::string>("next tag", flow_next_tag.get());
+      // now set the liquid water content as an interpolated field at next
+      // note that current copy is kept by transport PK
+      Teuchos::ParameterList& lwc_list_next = S_->GetEvaluatorList(Keys::getKey("liquid_water_content", transport_next_tag));
+      lwc_list_next.set<std::string>("evaluator type", "temporal interpolation");
+      lwc_list_next.set<std::string>("current tag", flow_current_tag.get());
+      lwc_list_next.set<std::string>("next tag", flow_next_tag.get());
 
+      // porosity used with velocity to compute particle velocity
       Teuchos::ParameterList& poro_list_next = S_->GetEvaluatorList(Keys::getKey("porosity", transport_next_tag));
       poro_list_next.set<std::string>("evaluator type", "alias");
       poro_list_next.set<std::string>("target", Keys::getKey("porosity", flow_next_tag, true));
     }
   }
 
   if (surface) {
-    auto [flow_current_tag, flow_next_tag] = tags_[0];
-    auto [transport_current_tag, transport_next_tag] = tags_[1];
-
     if (transport_next_tag != flow_next_tag) {
       // set the flow field evaluator as the flow's NEXT tag
       //
@@ -77,20 +77,40 @@ MPCFlowTransport::parseParameterList()
       flux_list.set<std::string>("evaluator type", "alias");
       flux_list.set<std::string>("target", Keys::getKey("surface-water_flux", flow_next_tag, true));
 
-      // set the surface-ponded_depth as an interpolated field
-      Teuchos::ParameterList& pd_list_current = S_->GetEvaluatorList(Keys::getKey("surface-ponded_depth", transport_current_tag));
-      pd_list_current.set<std::string>("evaluator type", "temporal interpolation");
-      pd_list_current.set<std::string>("current tag", flow_current_tag.get());
-      pd_list_current.set<std::string>("next tag", flow_next_tag.get());
-
-      Teuchos::ParameterList& pd_list_next = S_->GetEvaluatorList(Keys::getKey("surface-ponded_depth", transport_next_tag));
-      pd_list_next.set<std::string>("evaluator type", "temporal interpolation");
-      pd_list_next.set<std::string>("current tag", flow_current_tag.get());
-      pd_list_next.set<std::string>("next tag", flow_next_tag.get());
+      // velocity for dispersivity
+      Teuchos::ParameterList& velo_list = S_->GetEvaluatorList(Keys::getKey("surface-velocity", transport_next_tag));
+      velo_list.set<std::string>("evaluator type", "alias");
+      velo_list.set<std::string>("target", Keys::getKey("velocity", flow_next_tag, true));
+
+      // set the liquid water content as an interpolated field
+      // note that current copy is kept by transport PK
+      Teuchos::ParameterList& lwc_list_next = S_->GetEvaluatorList(Keys::getKey("surface-liquid_water_content", transport_next_tag));
+      lwc_list_next.set<std::string>("evaluator type", "temporal interpolation");
+      lwc_list_next.set<std::string>("current tag", flow_current_tag.get());
+      lwc_list_next.set<std::string>("next tag", flow_next_tag.get());
     }
   }
 
   MPCSubcycled::parseParameterList();
+
+  if (subsurface) {
+    // first require lwc@NEXT -- otherwise it will get called to be an alias
+    // lwc@TRANSPORT_NEXT, which is an interpolator that depends upon this.
+    requireAtNext("liquid_water_content", Tags::NEXT, *S_);
+
+    // now call at lwc@transport_next, which will be the interpolant
+    requireAtNext("liquid_water_content", transport_next_tag, *S_);
+  }
+
+  if (surface) {
+    // first require lwc@NEXT -- otherwise it will get called to be an alias
+    // lwc@TRANSPORT_NEXT, which is an interpolator that depends upon this.
+    requireAtNext("surface-liquid_water_content", Tags::NEXT, *S_);
+
+    // now call at lwc@transport_next, which will be the interpolant
+    requireAtNext("surface-liquid_water_content", transport_next_tag, *S_);
+  }
+
 }
 
 

src/pks/mpc/mpc_subcycled.cc

@@ -223,7 +223,9 @@ MPCSubcycled::CommitStep(double t_old, double t_new, const Tag& tag)
     // initial commit, also do the substep commits
     int i = 0;
     for (auto& pk : sub_pks_) {
-      if (subcycling_[i]) { pk->CommitStep(t_old, t_new, tags_[i].second); }
+      if (subcycling_[i]) {
+        pk->CommitStep(t_old, t_new, tags_[i].second);
+      }
       ++i;
     }
   }

src/pks/pk_bdf_default.cc

@@ -128,21 +128,6 @@ PK_BDF_Default::AdvanceStep(double t_old, double t_new, bool reinit)
   State_to_Solution(Tags::NEXT, *solution_);
 
   // take a bdf timestep
-  //
-  // Three dts:
-  // --  dt is the requested timestep size.  It must be less than or equal to...
-  // --  dt_internal is the max valid dt, and is set by physics/solvers
-  // --  dt_solver is what the solver wants to do
-  double dt_internal = S_->Get<double>(Keys::cleanName(name_, true) + "_dt_internal", Tags::DEFAULT);
-
-  // roundoff calculation follows that of Amanzi::Utils::isNearEqual() in
-  // Event.hh.  Unfortunately, still cannot have this assertion for manually
-  // overridden runs (where a file controls the timestep). --ETC
-  // if (dt > dt_internal) {
-  //   AMANZI_ASSERT(dt <= dt_internal + std::max(1., std::max(dt, dt_internal)) * 1.e4 * Utils::Event_EPS<double>::value);
-  // }
-
-  double dt_solver = -1;
   bool fail = false;
   try {
     fail = time_stepper_->AdvanceStep(dt, dt_next_, solution_);

src/pks/pk_physical_bdf_default.cc

@@ -28,6 +28,8 @@ PK_PhysicalBDF_Default::parseParameterList()
   max_valid_change_ = plist_->get<double>("max valid change", -1.0);
 
   conserved_key_ = Keys::readKey(*plist_, domain_, "conserved quantity");
+  requireAtCurrent(conserved_key_, tag_current_, *S_, name_);
+
   atol_ = plist_->get<double>("absolute error tolerance", 1.0);
   rtol_ = plist_->get<double>("relative error tolerance", 1.0);
   fluxtol_ = plist_->get<double>("flux error tolerance", 1.0);

src/pks/transport/transport_ats_pk.cc

@@ -98,15 +98,23 @@ Transport_ATS::parseParameterList()
   has_dispersion_ = plist_->get<bool>("has dispersion", false);
 
   // keys, dependencies, etc
+  // -- flux -- only needed at new time, evaluator controlled elsewhere
   flux_key_ = Keys::readKey(*plist_, domain_, "water flux", "water_flux");
+
+  // -- liquid water content - need at new time, copy at current time
   lwc_key_ = Keys::readKey(*plist_, domain_, "liquid water content", "liquid_water_content");
+  requireAtCurrent(lwc_key_, tag_current_, *S_, name_);
+
   water_src_key_ = Keys::readKey(*plist_, domain_, "water source", "water_source");
   cv_key_ = Keys::readKey(*plist_, domain_, "cell volume", "cell_volume");
 
-  // workspace
+  // workspace, no evaluator
   conserve_qty_key_ =
     Keys::readKey(*plist_, domain_, "conserved quantity", "total_component_quantity");
+  requireAtNext(conserve_qty_key_, tag_next_, *S_, name_);
+
   solid_residue_mass_key_ = Keys::readKey(*plist_, domain_, "solid residue", "solid_residue_mass");
+
   geochem_src_factor_key_ =
     Keys::readKey(*plist_, domain_, "geochem source factor", "geochem_src_factor");
 
@@ -491,7 +499,7 @@ Transport_ATS::SetupPhysicalEvaluators_()
     .SetMesh(mesh_)
     ->SetGhosted(true)
     ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1);
-  requireAtCurrent(lwc_key_, tag_current_, *S_);
+  requireAtCurrent(lwc_key_, tag_current_, *S_, name_);
 
   requireAtNext(molar_dens_key_, tag_next_, *S_)
     .SetMesh(mesh_)
@@ -747,14 +755,19 @@ Transport_ATS::AdvanceStep(double t_old, double t_new, bool reinit)
                << "----------------------------------------------------------------" << std::endl;
   AMANZI_ASSERT(std::abs(S_->get_time(tag_current_) - t_old) < 1.e-4);
   AMANZI_ASSERT(std::abs(S_->get_time(tag_next_) - t_new) < 1.e-4);
+  db_->WriteCellInfo(true);
 
   // check the stable step size again -- flow can now be computed at the
   // correct, new time, and may have changed, resulting in a smaller dt.
   //
   // By failing the step, it will get repeated with the smaller dt, again
   // potentially recomputing a flow field, but should be closer.
   dt_stable_ = ComputeStableTimeStep_();
-  if (dt > dt_stable_ + 1.e-4) return true;
+  if (dt > dt_stable_ + 1.e-4) {
+    if (vo_->os_OK(Teuchos::VERB_LOW))
+      *vo_->os() << "Failed step: requested dt = " << dt << " > stable dt = " << dt_stable_ << std::endl;
+    return true;
+  }
 
   if (vo_->os_OK(Teuchos::VERB_HIGH))
     *vo_->os() << "Water state:" << std::endl;
@@ -935,9 +948,8 @@ Transport_ATS::CommitStep(double t_old, double t_new, const Tag& tag_next)
 
   PK_Physical_Default::CommitStep(t_old, t_new, tag_next);
 
-  if (tag_next == Tags::NEXT) {
-    assign(lwc_key_, Tags::CURRENT, tag_next, *S_);
-  }
+  Tag tag_current = tag_next == tag_next_ ? tag_current_ : Tags::CURRENT;
+  assign(lwc_key_, tag_current, tag_next, *S_);
 }
 
 

src/pks/transport/transport_ats_ti.cc

@@ -270,7 +270,7 @@ Transport_ATS::InvertTccNew_(const Epetra_MultiVector& conserve_qty,
   // calculate the new conc based on advected term
   for (int c = 0; c < conserve_qty.MyLength(); c++) {
     double water_new = lwc_new[0][c];
-    double water_sink = conserve_qty[num_components_+1][c];
+    double water_sink = conserve_qty[num_aqueous_][c];
     double water_total = water_sink + water_new;
 
     for (int i = 0; i != num_components_; ++i) {

tools/input_converters/xml-1.5-master.py

@@ -111,11 +111,125 @@ def tensorPerm_(perm):
         else:
             perm.setParameter("tensor type", "string", "scalar")
 
+
+def fixTransportPK(pk, evals_list):
+    if pk.isElement("domain name"):
+        domain = pk.getElement("domain name").getValue()
+    else:
+        domain = "domain"
+
+
+    if not pk.isElement("primary variable key") and not pk.isElement("primary variable key suffix"):
+        pk.setParameter("primary variable key suffix", "string", "total_component_concentration")
+
+    if not pk.isElement("advection spatial discretization order") and pk.isElement("spatial discretization order"):
+        order = pk.getElement("spatial discretization order")
+        order.setName("advection spatial discretization order")
+
+    if pk.isElement("molecular diffusion") and not pk.isElement("molecular diffusivity [m^2 s^-1]"):
+        md = pk.pop("molecular diffusion")
+        names = md.getElement("aqueous names").getValue()
+        vals = md.getElement("aqueous values").getValue()
+        print("found MD for names:", names)
+        if len(names) > 0 and names[0] != '':
+            diff = pk.sublist("molecular diffusivity [m^2 s^-1]")
+            for name, val in zip(names, vals):
+                diff.setParameter(name, "string", val)
+
+    if pk.isElement("inverse") and pk.isElement("diffusion"):
+        inv = pk.pop("inverse")
+        pk.sublist("diffusion").append(inv)
+    elif pk.isElement("inverse"):
+        pk.pop("inverse")
+
+    if pk.isElement("transport subcycling"):
+        pk.pop("transport subcycling")
+
+    if pk.isElement("runtime diagnostics: regions"):
+        pk.pop("runtime diagnostics: regions")
+
+    if pk.isElement("material properties"):
+        mat_props = pk.pop("material properties").sublist("domain")
+        if mat_props.isElement("model"):
+            pk.setParameter("has dispersivity", "bool", True)
+
+            if domain == "domain":
+                disp_key = "dispersion_coefficient"
+            else:
+                disp_key = domain + "-dispersion_coefficient"
+
+            if not evals_list.isElement(disp_key):
+                disp_list = evals_list.sublist(disp_key)
+                disp_list.setParameter("evaluator type", "string", "dispersion tensor")
+                disp_list2 = disp_list.sublist("mechanical dispersion parameters").sublist(domain)
+                if domain == "domain":
+                    disp_list2.setParameter("region", "string", "computational domain")
+                else:
+                    disp_list2.setParameter("region", "string", domain + " domain")
+
+                disp_type = mat_props.getElement("model").getValue()
+                if disp_type == "scalar":
+                    disp_list2.setParameter("mechanical dispersion type", "string", "isotropic")
+                    params_list = mat_props.sublist(f"parameters for {disp_type}")
+                    params_list.setName("isotropic parameters")
+                    disp_list2.append(params_list)
+                else:
+                    disp_list2.setParameter("mechanical dispersion type", "string", disp_type)
+                    params_list = mat_props.sublist(f"parameters for {disp_type}")
+                    params_list.setName(f"{disp_type} parameters")
+                    disp_list2.append(params_list)
+
+        if mat_props.isElement("aqueous tortuosity"):
+            pk.sublist("tortuosity [-]").setParameter("aqueous", "double", mat_props.getElement("aqueous tortuosity").getValue())
+        if mat_props.isElement("gaseous tortuosity"):
+            pk.sublist("tortuosity [-]").setParameter("gaseous", "double", mat_props.getElement("gaseous tortuosity").getValue())
+
+    if pk.isElement("component molar masses"):
+        mms = pk.pop("component molar masses")
+        names = pk.getElement("component names").getValue()
+        for name, mm in zip(names, mms):
+            pk.sublist("molar mass [kg mol^-1]").setParameter(name, "double", mm)
+
+    if domain == "domain":
+        lwc_key = "liquid_water_content"
+    else:
+        lwc_key = f"{domain}-liquid_water_content"
+
+    print(f'searching for LWC = {lwc_key}')
+
+    if not evals_list.isElement(lwc_key):
+        print(' ... is not eval')
+        lwc_list = evals_list.sublist(lwc_key)
+        lwc_list.setParameter("evaluator type", "string", "multiplicative evaluator")
+
+        if "surface" in domain:
+            print('adding surface quantities for LWC')
+            lwc_list.setParameter("dependencies", "Array(string)", [f"{domain}-ponded_depth", f"{domain}-molar_density_liquid", f"{domain}-cell_volume"])
+        else:            
+            print('adding subsurface quantities for LWC')
+            lwc_list.setParameter("dependencies", "Array(string)", ["saturation_liquid", "molar_density_liquid", "porosity", "cell_volume"])
+
+
+
+
+def fixTransportPKs(xml):
+    evals_list = asearch.find_path(xml, ["state", "evaluators"], True)
+    pk_list = asearch.child_by_name(xml, "PKs")
+    for pk in pk_list:
+        pk_type = pk.getElement("PK type")
+        print(pk_type)
+
+        if pk_type.getValue() == "transport ATS":
+            print("fixing transport pk")
+            fixTransportPK(pk, evals_list)
+
+
 def update(xml):
     #enforceDtHistory(xml)
     timeStep(xml)
     tensorPerm(xml)
-
+    fixTransportPKs(xml)
+
 
 if __name__ == "__main__":
     import argparse