Merge branch 'master' into jb/elm_api

.github/workflows/ats-ci.yml

@@ -21,22 +21,24 @@ jobs:
         repository: amanzi/amanzi
         ref: master
         submodules: true
-    - name: Extract the Amanzi branch name
-      id: amanzi_branch
-      run:
-        echo "::set-output name=AMANZI_BRANCH::$(echo master)"
     - name: Extract the ATS branch name
       id: branch
       working-directory: Docker
       run:
-        echo "::set-output name=ATS_BRANCH::$(echo $GITHUB_REF | sed -e 's/refs\/heads\///')"
+        echo "::set-output name=ATS_BRANCH::$GITHUB_REF_NAME"
+    - name: Does Amanzi branch of the same name exist?
+      id: amanzi_branch
+      run:
+        echo "::set-output name=AMANZI_BRANCH::$(git ls-remote --heads origin ${GITHUB_REF_NAME} | sed 's/.*refs\/heads\///')"
     - name: Filter the branch name to generate a tag for Docker
       id: tag
       run:
         echo "::set-output name=ATS_BRANCH_TAG::$(echo ${{steps.branch.outputs.ATS_BRANCH}} | sed -e 's/\//--/g')"
     - name: Output the branch name
       run:
-        echo "Branch or Tag reference = ${{steps.branch.outputs.ATS_BRANCH}}"
+        echo "Amanzi branch = ${{steps.amanzi_branch.outputs.AMANZI_BRANCH}}";
+        echo "ATS branch = ${{steps.branch.outputs.ATS_BRANCH}}";
+        echo "Tag reference = ${{steps.tag.outputs.ATS_BRANCH_TAG}}"
     - name: Get TPLs version
       id: version 
       working-directory: Docker

src/constitutive_relations/column_integrators/EvaluatorColumnIntegrator.hh

@@ -138,10 +138,12 @@ EvaluatorColumnIntegrator<Parser,Integrator>::Evaluate_(const State& S,
     }
 
     // val[1] is typically e.g. cell volume, but can be 0 to indicate no
-    // denominator.
-    if (val[1] > 0.) res[0][col] = val[0] / val[1];
-    else res[0][col] = val[0];
+    // denominator.  Coefficient provides a hook for column-wide multiples
+    // (e.g. 1/surface area).
+    if (val[1] > 0.) res[0][col] = integrator.coefficient(col) * val[0] / val[1];
+    else res[0][col] = integrator.coefficient(col) * val[0];
   }
+
 }
 
 

src/constitutive_relations/column_integrators/thaw_depth_evaluator.cc

@@ -23,10 +23,10 @@ ParserThawDepth::ParserThawDepth(Teuchos::ParameterList& plist,
   Key temp_key = Keys::readKey(plist, domain_ss, "temperature", "temperature");
   dependencies.insert(KeyTag{temp_key, tag});
 
-  Key cv_key = Keys::readKey(plist, domain_ss, "cell volume", "cell_volume");
+  Key cv_key = Keys::readKey(plist, domain_ss, "subsurface cell volume", "cell_volume");
   dependencies.insert(KeyTag{cv_key, tag});
 
-  Key cv_surf_key = Keys::readKey(plist, domain, "cell volume", "cell_volume");
+  Key cv_surf_key = Keys::readKey(plist, domain, "surface cell volume", "cell_volume");
   dependencies.insert(KeyTag{cv_surf_key, tag});
 }
 

src/constitutive_relations/eos/eos_constant.cc

@@ -22,16 +22,16 @@ EOSConstant::EOSConstant(Teuchos::ParameterList& eos_plist) :
 
 void EOSConstant::InitializeFromPlist_() {
   // defaults to water
-  if (eos_plist_.isParameter("molar mass [kg/mol]")) {
-    M_ = eos_plist_.get<double>("molar mass [kg/mol]");
+  if (eos_plist_.isParameter("molar mass [kg mol^-1]")) {
+    M_ = eos_plist_.get<double>("molar mass [kg mol^-1]");
   } else {
-    M_ = eos_plist_.get<double>("molar mass [g/mol]", 18.0153) * 1.e-3;
+    M_ = eos_plist_.get<double>("molar mass [g mol^-1]", 18.0153) * 1.e-3;
   }
 
-  if (eos_plist_.isParameter("density [mol/m^3]")) {
-    rho_ = eos_plist_.get<double>("density [mol/m^3]") * M_;
+  if (eos_plist_.isParameter("density [mol m^-3]")) {
+    rho_ = eos_plist_.get<double>("density [mol m^-3]") * M_;
   } else {
-    rho_ = eos_plist_.get<double>("density [kg/m^3]");
+    rho_ = eos_plist_.get<double>("density [kg m^-3]");
   }
 
 };

src/constitutive_relations/generic_evaluators/AdditiveEvaluator.cc

@@ -20,16 +20,21 @@ AdditiveEvaluator::AdditiveEvaluator(Teuchos::ParameterList& plist) :
   }
 
   for (const auto& dep : dependencies_) {
-    Key varname = Keys::getKey(dep.first, dep.second);
-    if (plist.isParameter(varname+" coefficient")) {
-      coefs_[varname] = plist.get<double>(varname+" coefficient");
-    } else if (plist.isParameter(dep.first+" coefficient")) {
-      coefs_[varname] = plist.get<double>(dep.first+" coefficient");
+    Key variable = dep.first;
+    Key variable_tag = Keys::getKey(dep.first, dep.second);
+    Key varname = Keys::getVarName(variable);
+    if (plist.isParameter(variable_tag+" coefficient")) {
+      coefs_[variable_tag] = plist.get<double>(variable_tag+" coefficient");
+    } else if (plist.isParameter(variable+" coefficient")) {
+      coefs_[variable_tag] = plist.get<double>(variable+" coefficient");
+    } else if (plist.isParameter(varname+" coefficient")) {
+      coefs_[variable_tag] = plist.get<double>(varname+" coefficient");
     } else {
-      coefs_[varname] = 1.0;
+      coefs_[variable_tag] = 1.0;
     }
   }
   shift_ = plist.get<double>("constant shift", 0.);
+  positive_ = plist.get<bool>("enforce positivity", false);
 }
 
 
@@ -51,13 +56,33 @@ AdditiveEvaluator::Evaluate_(const State& S,
     double coef = coefs_[Keys::getKey(key_tag.first, key_tag.second)];
     result[0]->Update(coef, dep, 1.0);
   }
+
+  if (positive_) {
+    for (const auto& name : *result[0]) {
+      auto& res = *result[0]->ViewComponent(name, false);
+      for (int i=0; i!=res.MyLength(); ++i) res[0][i] = std::max(res[0][i], 0.);
+    }
+  }
 }
 
 void
 AdditiveEvaluator::EvaluatePartialDerivative_(const State& S,
         const Key& wrt_key, const Tag& wrt_tag, const std::vector<CompositeVector*>& result)
 {
   result[0]->PutScalar(coefs_[Keys::getKey(wrt_key, wrt_tag)]);
+
+  if (positive_) {
+    const auto& value = S.Get<CompositeVector>(my_keys_.front().first, my_keys_.front().second);
+    for (const auto& name : *result[0]) {
+      auto& res = *result[0]->ViewComponent(name, false);
+      const auto& value_v = *value.ViewComponent(name, false);
+      for (int i=0; i!=res.MyLength(); ++i) {
+        if (value_v[0][i] == 0.0) {
+          res[0][i] = 0.;
+        }
+      }
+    }
+  }
 }
 
 

src/constitutive_relations/generic_evaluators/AdditiveEvaluator.hh

@@ -13,6 +13,8 @@
 .. admonition:: additive-evaluator-spec
    * `"constant shift`" ``[double]`` **0** A constant value to add to the sum.
 
+   * `"enforce positivity`" ``[bool]`` **false** If true, max the result with 0.
+     
    * `"DEPENDENCY coefficient`" ``[double]`` **1.0** A multiple for each dependency in
      the list below.
 
@@ -52,6 +54,7 @@ class AdditiveEvaluator : public EvaluatorSecondaryMonotypeCV {
  protected:
   std::map<Key, double> coefs_;
   double shift_;
+  bool positive_;
 
  private:
   static Utils::RegisteredFactory<Evaluator,AdditiveEvaluator> factory_;

src/constitutive_relations/generic_evaluators/MultiplicativeEvaluator.cc

@@ -96,8 +96,12 @@ MultiplicativeEvaluator::EvaluatePartialDerivative_(const State& S,
       }
     }
     if (positive_) {
+      const auto& value_c = *S.Get<CompositeVector>(my_keys_.front().first, my_keys_.front().second)
+        .ViewComponent(lcv_name, false);
       for (int c=0; c!=res_c.MyLength(); ++c) {
-        res_c[c] = std::max(res_c[c], 0.);
+        if (value_c[0][c] == 0) {
+          res_c[c] = 0.;
+        }
       }
     }
   }

src/constitutive_relations/generic_evaluators/SubgridAggregateEvaluator.cc

@@ -70,11 +70,6 @@ SubgridAggregateEvaluator::EnsureEvaluators(State& S)
       msg << "SubgridAggregateEvaluator: DomainSet \"" << source_domain_ << "\" does not have a referencing parent but must have one to aggregate.";
       Exceptions::amanzi_throw(msg);
     }
-    if (S.GetMesh(domain_) != ds->get_referencing_parent()) {
-      Errors::Message msg;
-      msg << "SubgridAggregateEvaluator: DomainSet \"" << source_domain_ << "\" has a referencing parent, but it does not match the aggregate vector's domain, \"" << domain_ << "\"";
-      Exceptions::amanzi_throw(msg);
-    }
 
     for (const auto& subdomain : *ds) {
       dependencies_.insert(KeyTag{Keys::getKey(subdomain, var_key_), dep_tag});
@@ -85,11 +80,34 @@ SubgridAggregateEvaluator::EnsureEvaluators(State& S)
 }
 
 
+// Make sure that this vector is set on the referencing parent mesh of the
+// domain set.
+void
+SubgridAggregateEvaluator::EnsureCompatibility_Structure_(State& S)
+{
+  auto ds = S.GetDomainSet(source_domain_);
+  auto& dep_fac = S.Require<CompositeVector,CompositeVectorSpace>(dependencies_.front().first, dependencies_.front().second);
+  if (dep_fac.HasComponent("cell")) {
+    S.Require<CompositeVector,CompositeVectorSpace>(my_keys_.front().first, my_keys_.front().second)
+      .SetMesh(ds->get_referencing_parent())
+      ->AddComponent("cell", AmanziMesh::CELL, dep_fac.NumVectors("cell"));
+  }
+
+  if (S.GetRecordSet(dependencies_.front().first).subfieldnames()) {
+    S.GetRecordSetW(my_keys_.front().first).set_subfieldnames(*S.GetRecordSet(dependencies_.front().first).subfieldnames());
+  }
+}
+
+
 // Implements custom EC to use dependencies from subgrid
 void
 SubgridAggregateEvaluator::EnsureCompatibility_ToDeps_(State& S)
 {
-  EvaluatorSecondaryMonotypeCV::EnsureCompatibility_ToDeps_(S, {"cell",}, {AmanziMesh::CELL,}, {1,});
+  auto& fac = S.Require<CompositeVector,CompositeVectorSpace>(my_keys_.front().first, my_keys_.front().second);
+  if (fac.HasComponent("cell")) {
+    int num_vectors = fac.NumVectors("cell");
+    EvaluatorSecondaryMonotypeCV::EnsureCompatibility_ToDeps_(S, {"cell",}, {AmanziMesh::CELL,}, {num_vectors,});
+  }
 }
 
 

src/constitutive_relations/generic_evaluators/SubgridAggregateEvaluator.hh

@@ -45,6 +45,10 @@ class SubgridAggregateEvaluator : public EvaluatorSecondaryMonotypeCV {
   void EnsureEvaluators(State& S) override;
 
  protected:
+  // custom EC required to make sure that this vector is generated on the
+  // referencing mesh
+  void EnsureCompatibility_Structure_(State& S) override;
+
   // custom EC required to depend on cells of the subgrid mesh
   void EnsureCompatibility_ToDeps_(State& S) override;
 

src/constitutive_relations/generic_evaluators/SubgridDisaggregateEvaluator.cc

@@ -58,10 +58,14 @@ SubgridDisaggregateEvaluator::EvaluatePartialDerivative_(const State& S,
 void
 SubgridDisaggregateEvaluator::EnsureCompatibility_ToDeps_(State& S)
 {
-  CompositeVectorSpace fac;
-  fac.SetMesh(S.GetMesh(source_domain_))
-    ->AddComponent("cell", AmanziMesh::CELL, 1);
-  EvaluatorSecondaryMonotypeCV::EnsureCompatibility_ToDeps_(S, fac);
+  auto& my_fac = S.Require<CompositeVector,CompositeVectorSpace>(my_keys_.front().first, my_keys_.front().second);
+  if (my_fac.HasComponent("cell")) {
+    int num_vectors = my_fac.NumVectors("cell");
+    CompositeVectorSpace fac;
+    fac.SetMesh(S.GetMesh(source_domain_))
+      ->AddComponent("cell", AmanziMesh::CELL, num_vectors);
+    EvaluatorSecondaryMonotypeCV::EnsureCompatibility_ToDeps_(S, fac);
+  }
 }
 
 

src/executables/ats_mesh_factory.cc

@@ -72,7 +72,7 @@ createMeshFromFile(const std::string& mesh_name,
     if (mesh_plist.isParameter("build columns from set")) {
       std::string regionname = mesh_plist.get<std::string>("build columns from set");
       mesh->build_columns(regionname);
-    } else if (mesh_plist.get("build columns", false)) {
+    } else if (mesh_plist.get("build columns", true)) {
       mesh->build_columns();
     }
 

src/executables/ats_mesh_factory.hh

@@ -85,9 +85,9 @@ Example:
      <ParameterList name="domain">
        <Parameter name="mesh type" type="string" value="generate mesh"/>
        <ParameterList name="generate mesh parameters"/>
-         <Parameter name="number of cells" type="Array(int)" value="{{100, 1, 100}}"/>
-         <Parameter name="domain low coordinate" type="Array(double)" value="{{0.0, 0.0, 0.0}}" />
-         <Parameter name="domain high coordinate" type="Array(double)" value="{{100.0, 1.0, 10.0}}" />
+         <Parameter name="number of cells" type="Array(int)" value="{100, 1, 100}"/>
+         <Parameter name="domain low coordinate" type="Array(double)" value="{0.0, 0.0, 0.0}" />
+         <Parameter name="domain high coordinate" type="Array(double)" value="{100.0, 1.0, 10.0}" />
        </ParameterList>
      </ParameterList>   
    </ParameterList>   

src/executables/coordinator.cc

@@ -215,7 +215,6 @@ void Coordinator::initialize()
     }
   }
 
-
   // Final checks.
   //S_->CheckNotEvaluatedFieldsInitialized();
   S_->InitializeEvaluators();
@@ -482,6 +481,7 @@ Coordinator::get_dt(bool after_fail)
 
   // ask the step manager if this step is ok
   dt = tsm_->TimeStep(S_->get_time(Amanzi::Tags::NEXT), dt, after_fail);
+
   // note, I believe this can go away (along with the input spec flag) once
   // amanzi/amanzi#685 is closed --etc
   if (subcycled_ts_) dt = std::min(dt, dt_pk);
@@ -583,5 +583,4 @@ void Coordinator::checkpoint(bool force)
   }
 }
 
-
 } // close namespace Amanzi

src/executables/main.cc

@@ -88,6 +88,9 @@ int main(int argc, char *argv[])
   std::string verbosity;
   clp.setOption("verbosity", &verbosity, "Default verbosity level: \"none\", \"low\", \"medium\", \"high\", \"extreme\".");
 
+  std::string writing_rank;
+  clp.setOption("write_on_rank", &writing_rank, "Rank on which to write VerboseObjects");
+
   clp.throwExceptions(false);
   clp.recogniseAllOptions(true);
 
@@ -144,6 +147,24 @@ int main(int argc, char *argv[])
     return 1;
   }
 
+  // parse the writing rank
+  if (writing_rank.empty()) {
+    // pass
+  } else {
+    int writing_rank_j;
+    try {
+      writing_rank_j = std::stoi(writing_rank);
+    } catch (std::invalid_argument& e) {
+      std::cerr << "ERROR: invalid writing rank \"" << writing_rank << "\"" << std::endl;
+      clp.printHelpMessage("ats", std::cerr);
+    }
+    if (writing_rank_j < 0) {
+      std::cerr << "ERROR: invalid writing rank \"" << writing_rank << "\"" << std::endl;
+      clp.printHelpMessage("ats", std::cerr);
+    }
+    Amanzi::VerboseObject::global_writing_rank = writing_rank_j;
+  }
+
   // parse the input file and check validity
   if (input_filename.empty() && !opt_input_filename.empty()) input_filename = opt_input_filename;
   if (input_filename.empty()) {

src/pks/deform/volumetric_deformation.cc

@@ -29,10 +29,10 @@ VolumetricDeformation::VolumetricDeformation(Teuchos::ParameterList& pk_tree,
         const Teuchos::RCP<TreeVector>& solution)
   : PK(pk_tree, glist,  S, solution),
     PK_Physical_Default(pk_tree, glist,  S, solution),
-    surf_mesh_(Teuchos::null)
+    surf_mesh_(Teuchos::null),
+    deformed_this_step_(false)
 {
   dt_max_ = plist_->get<double>("max time step [s]", std::numeric_limits<double>::max());
-  dt_ = dt_max_;
 
   // The deformation mode describes how to calculate new cell volume from a
   // provided function and the old cell volume.

src/pks/deform/volumetric_deformation.hh

@@ -150,9 +150,7 @@ class VolumetricDeformation : public PK_Physical_Default {
     return dt_max_;
   }
 
-  virtual void set_dt(double dt) override {
-    dt_ = dt;
-  }
+  virtual void set_dt(double dt) override {}
 
  private:
 
@@ -188,7 +186,7 @@ class VolumetricDeformation : public PK_Physical_Default {
   double time_scale_, structural_vol_frac_;
 
   // PK timestep control
-  double dt_, dt_max_;
+  double dt_max_;
 
   // meshes
   Key domain_surf_, domain_surf_3d_;