major cleanup of PKs.

* moves multiple evaluators out of PKs and into the more typical evaluator list
* cleans up use of common geometric/mesh shortcut functions
* standardizes the use of keys in PKs

closes #49
closes #39
closes #38
closes #35
closes #33
closes #3

src/operators/upwinding/UpwindFluxFactory.cc

@@ -9,6 +9,8 @@
 
 #include "errors.hh"
 
+#include "Teuchos_ParameterList.hpp"
+
 #include "UpwindFluxFactory.hh"
 #include "upwind_total_flux.hh"
 #include "upwind_flux_harmonic_mean.hh"
@@ -19,23 +21,15 @@
 namespace Amanzi {
 namespace Operators {
 
-Teuchos::RCP<Upwinding> UpwindFluxFactory::Create(Teuchos::ParameterList& oplist,
-                                                  std::string pkname,
-                                                  std::string cell_coef,
-                                                  std::string face_coef,
-                                                  std::string flux) {
-  //  AMANZI_ASSERT(oplist.isSublist("overland conductivity model"));
-  Teuchos::ParameterList sublist;
-  if (oplist.isSublist("overland conductivity model")) {
-    sublist = oplist.sublist("overland conductivity model");
-  }
-  else {
-    sublist = oplist.sublist("overland conductivity subgrid model");
-  }
-
-  std::string model_type = sublist.get<std::string>("overland conductivity type", "manning");
-
-  double flux_eps = sublist.get<double>("upwind flux epsilon", 1.e-8);
+Teuchos::RCP<Upwinding>
+UpwindFluxFactory::Create(Teuchos::ParameterList& oplist,
+                          std::string pkname,
+                          std::string cell_coef,
+                          std::string face_coef,
+                          std::string flux)
+{
+  std::string model_type = oplist.get<std::string>("upwind type", "manning");
+  double flux_eps = oplist.get<double>("upwind flux epsilon", 1.e-8);
 
   if (model_type == "manning") {
     return Teuchos::rcp(new UpwindTotalFlux(pkname, cell_coef, face_coef, flux, flux_eps));
@@ -60,15 +54,15 @@ Teuchos::RCP<Upwinding> UpwindFluxFactory::Create(Teuchos::ParameterList& oplist
   } else if (model_type == "manning split denominator") {
     std::string slope = oplist.get<std::string>("slope key", "slope_magnitude");
     std::string manning_coef = oplist.get<std::string>("coefficient key", "manning_coefficient");
-    double slope_regularization = sublist.get<double>("slope regularization epsilon", 1.e-8);
+    double slope_regularization = oplist.get<double>("slope regularization epsilon", 1.e-8);
     std::string ponded_depth = oplist.get<std::string>("height key", "ponded_depth");
     return Teuchos::rcp(new UpwindFluxSplitDenominator(pkname, cell_coef, face_coef, flux, flux_eps, slope, manning_coef, slope_regularization, ponded_depth));
 
   } else if (model_type == "manning ponded depth passthrough") {
     std::string slope = oplist.get<std::string>("slope key", "slope_magnitude");
     std::string manning_coef = oplist.get<std::string>("coefficient key", "manning_coefficient");
-    double slope_regularization = sublist.get<double>("slope regularization epsilon", 1.e-8);
-    double manning_exp = sublist.get<double>("Manning exponent");
+    double slope_regularization = oplist.get<double>("slope regularization epsilon", 1.e-8);
+    double manning_exp = oplist.get<double>("Manning exponent");
     return Teuchos::rcp(new UpwindFluxFOCont(pkname, cell_coef, face_coef, flux, slope, manning_coef, "elevation", slope_regularization, manning_exp));
 
   } else if (model_type == "manning cell centered") {

src/pks/CMakeLists.txt

@@ -6,6 +6,7 @@
 #
 
 set(ats_pks_src_files
+  pk_helpers.cc
   pk_bdf_default.cc
   pk_physical_default.cc
   pk_physical_bdf_default.cc

src/pks/energy/constitutive_relations/thermal_conductivity/thermal_conductivity_surface_evaluator.cc

@@ -29,8 +29,8 @@ ThermalConductivitySurfaceEvaluator::ThermalConductivitySurfaceEvaluator(
 
   AMANZI_ASSERT(plist_.isSublist("thermal conductivity parameters"));
   Teuchos::ParameterList sublist = plist_.sublist("thermal conductivity parameters");
-  K_liq_ = sublist.get<double>("thermal conductivity of water [W/(m-K)]", 0.58);
-  K_ice_ = sublist.get<double>("thermal conductivity of ice [W/(m-K)]", 2.18);
+  K_liq_ = sublist.get<double>("thermal conductivity of water [W m^-1 K^-1]", 0.58);
+  K_ice_ = sublist.get<double>("thermal conductivity of ice [W m^-1 K^-1]", 2.18);
   min_K_ = sublist.get<double>("minimum thermal conductivity", 1.e-14);
 }
 

src/pks/energy/constitutive_relations/thermal_conductivity/thermal_conductivity_threephase_peterslidard.cc

@@ -37,10 +37,10 @@ void ThermalConductivityThreePhasePetersLidard::InitializeFromPlist_() {
   eps_ = plist_.get<double>("epsilon [-]", 1.e-10);
   alpha_u_ = plist_.get<double>("unsaturated alpha unfrozen [-]");
   alpha_f_ = plist_.get<double>("unsaturated alpha frozen [-]");
-  k_soil_ = plist_.get<double>("thermal conductivity of soil [W/(m-K)]");
-  k_ice_ = plist_.get<double>("thermal conductivity of ice [W/(m-K)]");
-  k_liquid_ = plist_.get<double>("thermal conductivity of liquid [W/(m-K)]");
-  k_gas_ = plist_.get<double>("thermal conductivity of gas [W/(m-K)]");
+  k_soil_ = plist_.get<double>("thermal conductivity of soil [W m^-1 K^-1]");
+  k_ice_ = plist_.get<double>("thermal conductivity of ice [W m^-1 K^-1]");
+  k_liquid_ = plist_.get<double>("thermal conductivity of liquid [W m^-1 K^-1]");
+  k_gas_ = plist_.get<double>("thermal conductivity of gas [W m^-1 K^-1]");
 };
 
 } // namespace Relations

src/pks/energy/constitutive_relations/thermal_conductivity/thermal_conductivity_threephase_sutra_hacked.cc

@@ -31,11 +31,10 @@ double ThermalConductivityThreePhaseSutraHacked::ThermalConductivity(double poro
 };
 
 void ThermalConductivityThreePhaseSutraHacked::InitializeFromPlist_() {
-  k_frozen_ = plist_.get<double>("thermal conductivity of frozen zone [W/(m-K)]");
-  k_unfrozen_ = plist_.get<double>("thermal conductivity of unfrozen zone [W/(m-K)]");
-  k_mushy_ = plist_.get<double>("thermal conductivity of mushy zone [W/(m-K)]");
+  k_frozen_ = plist_.get<double>("thermal conductivity of frozen zone [W m^-1 K^-1]");
+  k_unfrozen_ = plist_.get<double>("thermal conductivity of unfrozen zone [W m^-1 K^-1]");
+  k_mushy_ = plist_.get<double>("thermal conductivity of mushy zone [W m^-1 K^-1]");
   sr_ = plist_.get<double>("residual saturation [-]");
-
 };
 
 } // namespace Relations

src/pks/energy/constitutive_relations/thermal_conductivity/thermal_conductivity_threephase_volume_averaged.cc

@@ -29,10 +29,10 @@ double ThermalConductivityThreePhaseVolumeAveraged::ThermalConductivity(double p
 };
 
 void ThermalConductivityThreePhaseVolumeAveraged::InitializeFromPlist_() {
-  k_soil_ = plist_.get<double>("thermal conductivity of soil [W/(m-K)]");
-  k_ice_ = plist_.get<double>("thermal conductivity of ice [W/(m-K)]");
-  k_liquid_ = plist_.get<double>("thermal conductivity of liquid [W/(m-K)]");
-  k_gas_ = plist_.get<double>("thermal conductivity of gas [W/(m-K)]");
+  k_soil_ = plist_.get<double>("thermal conductivity of soil [W m^-1 K^-1]");
+  k_ice_ = plist_.get<double>("thermal conductivity of ice [W m^-1 K^-1]");
+  k_liquid_ = plist_.get<double>("thermal conductivity of liquid [W m^-1 K^-1]");
+  k_gas_ = plist_.get<double>("thermal conductivity of gas [W m^-1 K^-1]");
 };
 
 } // namespace Relations

src/pks/energy/constitutive_relations/thermal_conductivity/thermal_conductivity_threephase_wetdry.cc

@@ -75,8 +75,8 @@ void ThermalConductivityThreePhaseWetDry::InitializeFromPlist_() {
   eps_ = plist_.get<double>("epsilon [-]", 1.e-10);
   alpha_u_ = plist_.get<double>("unsaturated alpha unfrozen [-]");
   alpha_f_ = plist_.get<double>("unsaturated alpha frozen [-]");
-  k_dry_ = plist_.get<double>("thermal conductivity, dry [W/(m-K)]");
-  k_sat_u_ = plist_.get<double>("thermal conductivity, saturated (unfrozen) [W/(m-K)]");
+  k_dry_ = plist_.get<double>("thermal conductivity, dry [W m^-1 K^-1]");
+  k_sat_u_ = plist_.get<double>("thermal conductivity, saturated (unfrozen) [W m^-1 K^-1]");
   beta_sat_f_ = plist_.get<double>("saturated beta frozen [-]",1.0);
 };
 

src/pks/energy/constitutive_relations/thermal_conductivity/thermal_conductivity_twophase_wetdry.cc

@@ -38,8 +38,8 @@ double ThermalConductivityTwoPhaseWetDry::ThermalConductivity(double poro,
 void ThermalConductivityTwoPhaseWetDry::InitializeFromPlist_() {
   eps_ = plist_.get<double>("epsilon [-]", 1.e-10);
   alpha_ = plist_.get<double>("unsaturated alpha [-]", 1.0);
-  k_dry_ = plist_.get<double>("thermal conductivity, dry [W/(m-K)]");
-  k_wet_ = plist_.get<double>("thermal conductivity, wet [W/(m-K)]");
+  k_dry_ = plist_.get<double>("thermal conductivity, dry [W m^-1 K^-1]");
+  k_wet_ = plist_.get<double>("thermal conductivity, wet [W m^-1 K^-1]");
 };
 
 } // namespace Relations

src/pks/energy/energy_base.hh

@@ -1,6 +1,6 @@
 /*
-  ATS is released under the three-clause BSD License. 
-  The terms of use and "as is" disclaimer for this license are 
+  ATS is released under the three-clause BSD License.
+  The terms of use and "as is" disclaimer for this license are
   provided in the top-level COPYRIGHT file.
 
   Authors: Ethan Coon (ecoon@lanl.gov)
@@ -15,7 +15,7 @@ Solves an advection-diffusion equation for energy:
     \frac{\partial E}{\partial t} - \nabla \cdot \kappa \nabla T + \nabla \cdot \mathbf{q} e(T) = Q_w e(T) + Q_e
 
 .. todo:: Document the energy error norm!
-    
+
 .. _energy-pk-spec:
 .. admonition:: energy-pk-spec
 
@@ -27,24 +27,24 @@ Solves an advection-diffusion equation for energy:
 
     * `"boundary conditions`" ``[energy-bc-spec]`` Defaults to 0 diffusive flux
       boundary condition.  See `Energy-specific Boundary Conditions`_
-      
+
     * `"thermal conductivity evaluator`"
       ``[thermal-conductivity-evaluator-spec]`` The thermal conductivity.  This
       needs to go away, and should get moved to State.
 
     * `"absolute error tolerance`" ``[double]`` **76.e-6** A small amount of
       energy, see error norm. `[MJ]`
-      
+
     * `"upwind conductivity method`" ``[string]`` **arithmetic mean** Method of
       moving cell-based thermal conductivities onto faces.  One of:
 
       - `"arithmetic mean`" the default, average of neighboring cells
       - `"cell centered`" harmonic mean
 
     IF
-      
+
     * `"explicit advection`" ``[bool]`` **false** Treat the advection term implicitly.
-    
+
     ELSE
 
     * `"supress advective terms in preconditioner`" ``[bool]`` **false**
@@ -57,19 +57,19 @@ Solves an advection-diffusion equation for energy:
       Typically defaults are correct.
 
     END
-      
+
     * `"diffusion`" ``[pde-diffusion-spec]`` See PDE_Diffusion_, the diffusion operator.
 
     * `"diffusion preconditioner`" ``[pde-diffusion-spec]`` See
       PDE_Diffusion_, the inverse operator.  Typically only adds Jacobian
       terms, as all the rest default to those values from `"diffusion`".
 
     IF
-    
+
     * `"source term`" ``[bool]`` **false** Is there a source term?
-    
+
     THEN
-    
+
     * `"source key`" ``[string]`` **DOMAIN-total_energy_source** Typically
       not set, as the default is good. ``[MJ s^-1]``
 
@@ -84,11 +84,11 @@ Solves an advection-diffusion equation for energy:
     EVALUATORS:
 
     - `"source term`"
-    
+
     END
 
     Globalization:
-    
+
     * `"modify predictor with consistent faces`" ``[bool]`` **false** In a
       face+cell diffusion discretization, this modifies the predictor to make
       sure that faces, which are a DAE, are consistent with the predicted cells
@@ -101,7 +101,7 @@ Solves an advection-diffusion equation for energy:
       0, stops nonlinear updates from being too big through clipping.
 
     The following are rarely set by the user, as the defaults are typically right.
-    
+
     Variable names:
 
     * `"conserved quantity key`" ``[string]`` **DOMAIN-energy** The total energy :math:`E` `[MJ]`
@@ -113,24 +113,24 @@ Solves an advection-diffusion equation for energy:
     * `"advected energy flux`" ``[string]`` **DOMAIN-advected_energy_flux** :math:`\mathbf{q_e^{adv}} = q e` `[MJ s^-1]`
     * `"thermal conductivity`" ``[string]`` **DOMAIN-thermal_conductivity** Thermal conductivity on cells `[W m^-1 K^-1]`
     * `"upwinded thermal conductivity`" ``[string]`` **DOMAIN-upwinded_thermal_conductivity** Thermal conductivity on faces `[W m^-1 K^-1]`
-    
+
     * `"advection`" ``[pde-advection-spec]`` **optional** The PDE_Advection_ spec.  Only one current implementation, so defaults are typically fine.
-    
+
     * `"accumulation preconditioner`" ``[pde-accumulation-spec]`` **optional**
       The inverse of the accumulation operator.  See PDE_Accumulation_.
       Typically not provided by users, as defaults are correct.
 
     IF
-    
+
     * `"coupled to surface via flux`" ``[bool]`` **false** If true, apply
       surface boundary conditions from an exchange flux.  Note, if this is a
       coupled problem, it is probably set by the MPC.  No need for a user to
       set it.
-      
+
     THEN
 
     * `"surface-subsurface energy flux key`" ``[string]`` **DOMAIN-surface_subsurface_energy_flux**
-    
+
     END
 
     * `"coupled to surface via temperature`" ``[bool]`` **false** If true, apply
@@ -144,7 +144,7 @@ Solves an advection-diffusion equation for energy:
     - `"thermal conductivity`"
     - `"conserved quantity`"
     - `"energy`"
-    
+
  */
 
 
@@ -178,7 +178,7 @@ public:
              const Teuchos::RCP<Teuchos::ParameterList>& plist,
              const Teuchos::RCP<State>& S,
              const Teuchos::RCP<TreeVector>& solution);
-  
+
   // Virtual destructor
   virtual ~EnergyBase() {}
 
@@ -215,7 +215,7 @@ public:
 
   virtual bool ModifyPredictor(double h, Teuchos::RCP<const TreeVector> u0,
           Teuchos::RCP<TreeVector> u) override;
-    
+
   // evaluating consistent faces for given BCs and cell values
   virtual void CalculateConsistentFaces(const Teuchos::Ptr<CompositeVector>& u);
 
@@ -233,7 +233,6 @@ public:
   //    will get replaced by a better system when we get maps on the boundary
   //    faces.
   virtual void ApplyDirichletBCsToEnthalpy_(const Teuchos::Ptr<State>& S);
-  virtual void ApplyDirichletBCsToBoundaryFace_(const Teuchos::Ptr<CompositeVector>& temp);
 
   // -- Add any source terms into the residual.
   virtual void AddSources_(const Teuchos::Ptr<State>& S,
@@ -263,9 +262,6 @@ public:
   virtual void ApplyDiffusion_(const Teuchos::Ptr<State>& S,
           const Teuchos::Ptr<CompositeVector>& f);
 
-  virtual int BoundaryFaceGetCell(int f) const;
-
-
  protected:
   int niter_;
 
@@ -306,14 +302,14 @@ public:
   double T_limit_;
   double mass_atol_;
   double soil_atol_;
-  
+
   bool coupled_to_subsurface_via_temp_;
   bool coupled_to_subsurface_via_flux_;
   bool coupled_to_surface_via_temp_;
   bool coupled_to_surface_via_flux_;
 
   // Keys
-  Key energy_key_;
+  Key ss_primary_key_;
   Key wc_key_;
   Key enthalpy_key_;
   Key flux_key_;

src/pks/energy/energy_base_physics.cc

@@ -26,12 +26,12 @@ void EnergyBase::AddAccumulation_(const Teuchos::Ptr<CompositeVector>& g) {
   double dt = S_next_->time() - S_inter_->time();
 
   // update the energy at both the old and new times.
-  S_next_->GetFieldEvaluator(energy_key_)->HasFieldChanged(S_next_.ptr(), name_);
-  S_inter_->GetFieldEvaluator(energy_key_)->HasFieldChanged(S_inter_.ptr(), name_);
+  S_next_->GetFieldEvaluator(conserved_key_)->HasFieldChanged(S_next_.ptr(), name_);
+  S_inter_->GetFieldEvaluator(conserved_key_)->HasFieldChanged(S_inter_.ptr(), name_);
 
   // get the energy at each time
-  Teuchos::RCP<const CompositeVector> e1 = S_next_->GetFieldData(energy_key_);
-  Teuchos::RCP<const CompositeVector> e0 = S_inter_->GetFieldData(energy_key_);
+  Teuchos::RCP<const CompositeVector> e1 = S_next_->GetFieldData(conserved_key_);
+  Teuchos::RCP<const CompositeVector> e0 = S_inter_->GetFieldData(conserved_key_);
 
   // Update the residual with the accumulation of energy over the
   // timestep, on cells.