improves memory management, initialization of reconstruction objects

src/pks/transport/transport_ats.hh

@@ -456,7 +456,6 @@ class Transport_ATS : public PK_Physical_Default {
   Teuchos::RCP<Epetra_IntVector> upwind_cell_, downwind_cell_;
   Teuchos::RCP<Operators::ReconstructionCellLinear> lifting_;
   Teuchos::RCP<Operators::LimiterCell> limiter_;
-  int limiter_model_;
 
   Teuchos::RCP<Operators::BCs> adv_bcs_;
   std::vector<Teuchos::RCP<TransportDomainFunction>> bcs_;

src/pks/transport/transport_ats_pk.cc

@@ -95,7 +95,6 @@ Transport_ATS::parseParameterList()
 
     tortuosity_ = readParameterMapByPhase(plist_->sublist("tortuosity [-]"), 1.);
   }
-  has_dispersion_ = plist_->get<bool>("has dispersion", false);
 
   // keys, dependencies, etc
   // -- flux -- only needed at new time, evaluator controlled elsewhere
@@ -123,6 +122,7 @@ Transport_ATS::parseParameterList()
 
   // dispersion coefficient tensor
   dispersion_tensor_key_ = Keys::readKey(*plist_, domain_, "dispersion coefficient", "dispersion_coefficient");
+  has_dispersion_ = S_->HasEvaluatorList(dispersion_tensor_key_);
 
   // other parameters
   // -- a small amount of water, used to define when we are going to completely dry out a grid cell
@@ -221,8 +221,19 @@ Transport_ATS::SetupTransport_()
 
   if (adv_spatial_disc_order_ == 2) {
     // reconstruction initialization
-    limiter_ = Teuchos::rcp(new Operators::LimiterCell(mesh_));
+    Teuchos::ParameterList& recon_list = plist_->sublist("reconstruction");
+
+    // check and set defaults
+    if (!recon_list.isParameter("limiter extention for transport"))
+      recon_list.set<bool>("limiter extention for transport", true);
+    if (!recon_list.isParameter("limiter"))
+      recon_list.set<std::string>("limiter", "tensorial");
+
     lifting_ = Teuchos::rcp(new Operators::ReconstructionCellLinear(mesh_));
+    lifting_->Init(recon_list);
+
+    limiter_ = Teuchos::rcp(new Operators::LimiterCell(mesh_));
+    limiter_->Init(recon_list);
   }
 
   adv_bcs_ = Teuchos::rcp(
@@ -901,8 +912,6 @@ Transport_ATS ::AdvanceDispersionDiffusion_(double t_old, double t_new)
 
     // -- build the matrices if needed
     if (changed_tensor || i == 0) {
-      std::cout << "D = " << (*D_)[0] << std::endl;
-
       // update mass, stiffness matrices of diffusion operator
       diff_global_op_->Init();
       diff_op_->SetTensorCoefficient(Teuchos::rcpFromRef(D_->data));

src/pks/transport/transport_ats_ti.cc

@@ -160,24 +160,18 @@ Transport_ATS::AddAdvection_SecondOrderUpwind_(double t_old,
   int nfaces_all =
     mesh_->getNumEntities(AmanziMesh::Entity_kind::FACE, AmanziMesh::Parallel_kind::ALL);
 
-  // should Init be moved to Setup?  Or does it also zero things out? --ETC
-  Teuchos::ParameterList recon_list = plist_->sublist("reconstruction");
-  lifting_->Init(recon_list);
-
-  // ugly API funkiness
-  auto tcc_tmp = Teuchos::rcp(new Epetra_Vector(tcc));
-  *tcc_tmp = tcc;
-  lifting_->Compute(tcc_tmp, 0);
+  // API funkiness
+  Teuchos::RCP<const Epetra_MultiVector> tcc_rcp = Teuchos::rcpFromRef<const Epetra_MultiVector>(tcc);
+  lifting_->Compute(tcc_rcp, 0);
 
   // populate boundary conditions for the current component
   PopulateBoundaryData_(component, *adv_bcs_);
   auto& bc_model = adv_bcs_->bc_model();
   auto& bc_value = adv_bcs_->bc_value();
 
   auto flux = S_->Get<CompositeVector>(flux_key_, tag_next_).ViewComponent("face", true);
-  limiter_->Init(recon_list, flux);
-  limiter_->ApplyLimiter(tcc_tmp, 0, lifting_, bc_model, bc_value);
-  lifting_->data()->ScatterMasterToGhosted("cell");
+  limiter_->SetFlux(flux);
+  limiter_->ApplyLimiter(tcc_rcp, 0, lifting_, bc_model, bc_value);
 
   // ADVECTIVE FLUXES
   // We assume that limiters made their job up to round-off errors.
@@ -188,14 +182,14 @@ Transport_ATS::AddAdvection_SecondOrderUpwind_(double t_old,
 
     double u1, u2, umin, umax;
     if (c1 >= 0 && c2 >= 0) {
-      u1 = (*tcc_tmp)[c1];
-      u2 = (*tcc_tmp)[c2];
+      u1 = tcc[c1];
+      u2 = tcc[c2];
       umin = std::min(u1, u2);
       umax = std::max(u1, u2);
     } else if (c1 >= 0) {
-      u1 = u2 = umin = umax = (*tcc_tmp)[c1];
+      u1 = u2 = umin = umax = tcc[c1];
     } else if (c2 >= 0) {
-      u1 = u2 = umin = umax = (*tcc_tmp)[c2];
+      u1 = u2 = umin = umax = tcc[c2];
     }
 
     double u = std::abs((*flux)[0][f]);
@@ -220,7 +214,7 @@ Transport_ATS::AddAdvection_SecondOrderUpwind_(double t_old,
       cq_flux[c1] -= tcc_flux;
 
     } else if (c1 >= 0 && c1 < cq_flux.MyLength() && (c2 < 0)) {
-      upwind_tcc = (*tcc_tmp)[c1];
+      upwind_tcc = tcc[c1];
       upwind_tcc = std::max(upwind_tcc, umin);
       upwind_tcc = std::min(upwind_tcc, umax);
 
@@ -290,10 +284,6 @@ Transport_ATS::InvertTccNew_(const Epetra_MultiVector& conserve_qty,
         // at the new time + stuff leaving through the domain coupling, divided
         // by water of both
         tcc[i][c] = conserve_qty[i][c] / water_total;
-        if (i == 0 && (c == 0 || c == 99)) {
-          std::cout << "inverting: conserve_qty = " << conserve_qty[i][c] << " / water_sink = "
-                    << water_sink << " + water_new = " << water_new << " = water_total = " << water_total << " --> tcc = " << tcc[i][c] << std::endl;
-        }
         if (include_current_water_mass) conserve_qty[num_aqueous_][c] = water_total;
       } else if (water_sink > water_tolerance_ / cv[0][c] && conserve_qty[i][c] > 0) {
         // there is water and stuff leaving through the domain coupling, but it

tools/input_converters/xml-1.5-master.py

@@ -151,8 +151,6 @@ def fixTransportPK(pk, evals_list):
     if pk.isElement("material properties"):
         mat_props = pk.pop("material properties").sublist("domain")
         if mat_props.isElement("model"):
-            pk.setParameter("has dispersion", "bool", True)
-
             if domain == "domain":
                 disp_key = "dispersion_coefficient"
             else: