Merge branch 'master' into phongle/optimize_transport

amanzi · Mar 21, 2024 · c4702b6 · c4702b6
2 parents d10a81c + e790ec8
commit c4702b6
Show file tree

Hide file tree

Showing 14 changed files with 126 additions and 99 deletions.
diff --git a/src/constitutive_relations/column_integrators/ColumnSumEvaluator.cc b/src/constitutive_relations/column_integrators/ColumnSumEvaluator.cc
@@ -62,12 +62,12 @@ IntegratorColumnSum::IntegratorColumnSum(Teuchos::ParameterList& plist,
   divide_by_density_ = plist.get<bool>("divide by density");
 
   if (volume_factor_) {
-    AMANZI_ASSERT(deps.size() > 3);
+    AMANZI_ASSERT(deps.size() >= 3);
     cv_ = deps[i_dep++];
     surf_cv_ = deps[i_dep++];
   }
   if (volume_average_) {
-    AMANZI_ASSERT(deps.size() > 2);
+    AMANZI_ASSERT(deps.size() >= 2);
     cv_ = deps[i_dep++];
   }
 

diff --git a/src/pks/flow/constitutive_relations/overland_conductivity/unfrozen_fraction_model.hh b/src/pks/flow/constitutive_relations/overland_conductivity/unfrozen_fraction_model.hh
@@ -22,7 +22,7 @@
    * `"freezing point [K]`" ``[double]`` **273.15** Center of the transition,
      at this point unfrozen fraction is 0.5.
 
-   * `"minimum unfrozen fraction [-]`` ``[double]`` **0** Sets a minimum value.
+   * `"minimum unfrozen fraction [-]`" ``[double]`` **0** Sets a minimum value.
 
 */
 

diff --git a/src/pks/flow/constitutive_relations/sources/soil_resistance_sellers_evaluator.hh b/src/pks/flow/constitutive_relations/sources/soil_resistance_sellers_evaluator.hh
@@ -14,8 +14,8 @@ referred to Sellers et al. (1992).
 
 `"evaluator type`" = `"soil resistance, Sellers`"
 
-.. _compressible-porosity-evaluator-spec
-.. admonition:: compressible-porosity-evaluator-spec
+.. _soil-resistance-sellers-evaluator-spec
+.. admonition:: soil-resistance-sellers-evaluator-spec
 
    KEYS:
 

diff --git a/src/pks/flow/constitutive_relations/wrm/wrm.hh b/src/pks/flow/constitutive_relations/wrm/wrm.hh
@@ -21,10 +21,9 @@ WRMs are available in a collection of types:
 .. admonition:: wrm-typedinline-spec
 
    * `"region`" ``[string]`` Region to which this applies
-   * `"wrm type`" ``[string]`` Type of the WRM.  One of:
-
-     - `"van Genuchten`"
-     - `"linear system`" saturation a linear function of pressure
+   * `"wrm type`" ``[string]`` Type of the WRM.  One of `"van Genuchten`",
+     `"Brooks-Corey`", or `"linear system`", a saturation given as a linear
+     function of pressure
 
 */
 

diff --git a/src/pks/flow/constitutive_relations/wrm/wrm_brooks_corey.hh b/src/pks/flow/constitutive_relations/wrm/wrm_brooks_corey.hh
@@ -18,15 +18,12 @@ to Brooks-Corey parameters.
 .. _WRM-Brooks-Corey-spec
 .. admonition:: WRM-Brooks-Corey-spec
 
-    * `"region`" ``[string]`` Region to which this applies
-    * `"Brooks-Corey lambda [-]`" ``[double]`` Brooks-Corey's lambda
-    * `"Brooks-Corey saturated matric suction [Pa]`" ``[double]`` Brooks-Corey
-    saturated matric suction in Pa
-
-    END
-
-    * `"residual saturation [-]`" ``[double]`` **0.0**
-    * `"smoothing interval width [saturation]`" ``[double]`` **0.0**
+   * `"region`" ``[string]`` Region to which this applies
+   * `"Brooks-Corey lambda [-]`" ``[double]`` Brooks-Corey's lambda
+   * `"Brooks-Corey saturated matric suction [Pa]`" ``[double]`` Brooks-Corey
+     saturated matric suction in Pa
+   * `"residual saturation [-]`" ``[double]`` **0.0**
+   * `"smoothing interval width [saturation]`" ``[double]`` **0.0**
 
 Example:
 

diff --git a/src/pks/flow/constitutive_relations/wrm/wrm_evaluator.hh b/src/pks/flow/constitutive_relations/wrm/wrm_evaluator.hh
@@ -19,7 +19,7 @@ commonly used in practice is the van Genuchten model, but others are available d
 .. _wrm-evaluator-spec:
 .. admonition:: wrm-evaluator-spec
 
-   * `"model parameters`" ``[string]``  **WRM parameters** ``[WRM-typedinline-spec-list]``
+   * `"model parameters`" ``[string]``
      List (by region) of WRM specs. This will copy `"WRM parameters`" given in `"model parameters`"
      under state here to evaluate WRM.
 

diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/albedo_threecomponent_evaluator.cc b/src/pks/surface_balance/constitutive_relations/land_cover/albedo_threecomponent_evaluator.cc
@@ -20,6 +20,16 @@ namespace Relations {
 AlbedoThreeComponentEvaluator::AlbedoThreeComponentEvaluator(Teuchos::ParameterList& plist)
   : EvaluatorSecondaryMonotypeCV(plist)
 {
+  // parameters
+  a_ice_ = plist_.get<double>("albedo ice [-]", 0.44);
+  a_water_ = plist_.get<double>("albedo water [-]", 0.1168);
+  is_constant_snow_albedo_ = plist_.isParameter("albedo snow [-]");
+  if (is_constant_snow_albedo_) { a_snow_ = plist_.get<double>("albedo snow [-]"); }
+
+  e_ice_ = plist_.get<double>("emissivity ice [-]", 0.98);
+  e_water_ = plist_.get<double>("emissivity water [-]", 0.995);
+  e_snow_ = plist_.get<double>("emissivity snow [-]", 0.98);
+
   // determine the domain
   Key akey = my_keys_.front().first;
   domain_ = Keys::getDomain(akey);
@@ -40,22 +50,14 @@ AlbedoThreeComponentEvaluator::AlbedoThreeComponentEvaluator(Teuchos::ParameterL
 
   // dependencies
   // -- snow properties
-  snow_dens_key_ = Keys::readKey(plist, domain_snow_, "snow density", "density");
-  dependencies_.insert(KeyTag{ snow_dens_key_, tag });
+  if (!is_constant_snow_albedo_) {
+    snow_dens_key_ = Keys::readKey(plist, domain_snow_, "snow density", "density");
+    dependencies_.insert(KeyTag{ snow_dens_key_, tag });
+  }
 
   // -- skin properties
   unfrozen_fraction_key_ = Keys::readKey(plist, domain_, "unfrozen fraction", "unfrozen_fraction");
   dependencies_.insert(KeyTag{ unfrozen_fraction_key_, tag });
-
-  // parameters
-  a_ice_ = plist_.get<double>("albedo ice [-]", 0.44);
-  a_water_ = plist_.get<double>("albedo water [-]", 0.1168);
-  is_constant_snow_albedo_ = plist_.isParameter("albedo snow [-]");
-  if (is_constant_snow_albedo_) { a_snow_ = plist_.get<double>("albedo snow [-]"); }
-
-  e_ice_ = plist_.get<double>("emissivity ice [-]", 0.98);
-  e_water_ = plist_.get<double>("emissivity water [-]", 0.995);
-  e_snow_ = plist_.get<double>("emissivity snow [-]", 0.98);
 }
 
 // Required methods from EvaluatorSecondaryMonotypeCV
@@ -113,7 +115,10 @@ AlbedoThreeComponentEvaluator::EvaluatePartialDerivative_(
 void
 AlbedoThreeComponentEvaluator::EnsureCompatibility_Structure_(State& S)
 {
-  S.GetRecordSetW(my_keys_.front().first).set_subfieldnames({ "bare", "water", "snow" });
+  EnsureCompatibility_StructureSame_(S);
+  for (const auto& key : my_keys_) {
+    S.GetRecordSetW(key.first).set_subfieldnames({ "bare", "water", "snow" });
+  }
 }
 
 

diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/albedo_twocomponent_evaluator.cc b/src/pks/surface_balance/constitutive_relations/land_cover/albedo_twocomponent_evaluator.cc
@@ -20,6 +20,16 @@ namespace Relations {
 AlbedoTwoComponentEvaluator::AlbedoTwoComponentEvaluator(Teuchos::ParameterList& plist)
   : EvaluatorSecondaryMonotypeCV(plist)
 {
+  // parameters
+  a_ice_ = plist_.get<double>("albedo ice [-]", 0.44);
+  a_water_ = plist_.get<double>("albedo water [-]", 0.1168);
+  is_constant_snow_albedo_ = plist_.isParameter("albedo snow [-]");
+  if (is_constant_snow_albedo_) { a_snow_ = plist_.get<double>("albedo snow [-]"); }
+
+  e_ice_ = plist_.get<double>("emissivity ice [-]", 0.98);
+  e_water_ = plist_.get<double>("emissivity water [-]", 0.995);
+  e_snow_ = plist_.get<double>("emissivity snow [-]", 0.98);
+
   // determine the domain
   Key akey = my_keys_.front().first;
   domain_ = Keys::getDomain(akey);
@@ -40,8 +50,10 @@ AlbedoTwoComponentEvaluator::AlbedoTwoComponentEvaluator(Teuchos::ParameterList&
 
   // dependencies
   // -- snow properties
-  snow_dens_key_ = Keys::readKey(plist, domain_snow_, "snow density", "density");
-  dependencies_.insert(KeyTag{ snow_dens_key_, tag });
+  if (!is_constant_snow_albedo_) {
+    snow_dens_key_ = Keys::readKey(plist, domain_snow_, "snow density", "density");
+    dependencies_.insert(KeyTag{ snow_dens_key_, tag });
+  }
 
   // -- skin properties
   unfrozen_fraction_key_ = Keys::readKey(plist, domain_, "unfrozen fraction", "unfrozen_fraction");
@@ -50,16 +62,6 @@ AlbedoTwoComponentEvaluator::AlbedoTwoComponentEvaluator(Teuchos::ParameterList&
   // -- skin properties
   ponded_depth_key_ = Keys::readKey(plist, domain_, "ponded depth", "ponded_depth");
   dependencies_.insert(KeyTag{ ponded_depth_key_, tag });
-
-  // parameters
-  a_ice_ = plist_.get<double>("albedo ice [-]", 0.44);
-  a_water_ = plist_.get<double>("albedo water [-]", 0.1168);
-  is_constant_snow_albedo_ = plist_.isParameter("albedo snow [-]");
-  if (is_constant_snow_albedo_) { a_snow_ = plist_.get<double>("albedo snow [-]"); }
-
-  e_ice_ = plist_.get<double>("emissivity ice [-]", 0.98);
-  e_water_ = plist_.get<double>("emissivity water [-]", 0.995);
-  e_snow_ = plist_.get<double>("emissivity snow [-]", 0.98);
 }
 
 // Required methods from EvaluatorSecondaryMonotypeCV
@@ -127,7 +129,10 @@ AlbedoTwoComponentEvaluator::EvaluatePartialDerivative_(
 void
 AlbedoTwoComponentEvaluator::EnsureCompatibility_Structure_(State& S)
 {
-  S.GetRecordSetW(my_keys_.front().first).set_subfieldnames({ "bare_or_water", "snow" });
+  EnsureCompatibility_StructureSame_(S);
+  for (const auto& key : my_keys_) {
+    S.GetRecordSetW(key.first).set_subfieldnames({ "bare_or_water", "snow" });
+  }
 }
 
 void

diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/canopy_radiation_evaluator.hh b/src/pks/surface_balance/constitutive_relations/land_cover/canopy_radiation_evaluator.hh
@@ -85,6 +85,9 @@ class CanopyRadiationEvaluator : public EvaluatorSecondaryMonotypeCV {
 
  protected:
   virtual void EnsureCompatibility_ToDeps_(State& S) override;
+  virtual void EnsureCompatibility_Structure_(State& S) {
+    EnsureCompatibility_StructureSame_(S);
+  }
 
   // Required methods from EvaluatorSecondaryMonotypeCV
   virtual void Evaluate_(const State& S, const std::vector<CompositeVector*>& results) override;

diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/drainage_evaluator.hh b/src/pks/surface_balance/constitutive_relations/land_cover/drainage_evaluator.hh
@@ -72,6 +72,10 @@ class DrainageEvaluator : public EvaluatorSecondaryMonotypeCV {
                                           const Tag& wrt_tag,
                                           const std::vector<CompositeVector*>& results) override;
 
+  virtual void EnsureCompatibility_Structure_(State& S) {
+    EnsureCompatibility_StructureSame_(S);
+  }
+
  protected:
   Key drainage_key_;
   Key fracwet_key_;

diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/interception_fraction_evaluator.hh b/src/pks/surface_balance/constitutive_relations/land_cover/interception_fraction_evaluator.hh
@@ -74,6 +74,9 @@ class InterceptionFractionEvaluator : public EvaluatorSecondaryMonotypeCV {
                                           const Key& wrt_key,
                                           const Tag& wrt_tag,
                                           const std::vector<CompositeVector*>& results) override;
+  virtual void EnsureCompatibility_Structure_(State& S) {
+    EnsureCompatibility_StructureSame_(S);
+  }
 
   void InitializeFromPlist_();
 

diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/radiation_balance_evaluator.hh b/src/pks/surface_balance/constitutive_relations/land_cover/radiation_balance_evaluator.hh
@@ -14,6 +14,9 @@ Here the net radiation is positive for energy inputs to the layer.  Note that
 ground is based on the two-channel (land + snow) while canopy is assumed to be
 a simple, single layer.
 
+This evaluator requires that the surface temperature, snow temperature, and
+canopy temperature are known, or at least being solved for.
+
 Requires the use of LandCover types, for albedo and Beer's law coefficients.
 
 This is combination of CLM v4.5 Tech Note and Beer's law for attenuation of
@@ -90,6 +93,9 @@ class RadiationBalanceEvaluator : public EvaluatorSecondaryMonotypeCV {
 
  protected:
   virtual void EnsureCompatibility_ToDeps_(State& S) override;
+  virtual void EnsureCompatibility_Structure_(State& S) {
+    EnsureCompatibility_StructureSame_(S);
+  }
 
   // Required methods from EvaluatorSecondaryMonotypeCV
   virtual void Evaluate_(const State& S, const std::vector<CompositeVector*>& results) override;

diff --git a/src/pks/surface_balance/surface_balance_base.cc b/src/pks/surface_balance/surface_balance_base.cc
@@ -96,7 +96,13 @@ SurfaceBalanceBase::Setup()
       .SetMesh(mesh_)
       ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1);
     S_->RequireDerivative<CompositeVector, CompositeVectorSpace>(
-      conserved_key_, tag_next_, key_, tag_next_);
+      conserved_key_, tag_next_, key_, tag_next_)
+          .SetMesh(mesh_)
+          ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1);
+        // NOTE, remove SetMesh/AddComponent lines after fixing amanzi/ats#167.
+        // The mesh should get set by the evaluator, but when
+        // the evaluator isn't actually differentiable, it
+        // doesn't get done.
 
     //    and at the current time, where it is a copy evaluator
     requireAtCurrent(conserved_key_, tag_current_, *S_, name_);