diff --git a/.github/workflows/ats-ci.yml b/.github/workflows/ats-ci.yml index 3213a2fbb..9cdd77bfb 100644 --- a/.github/workflows/ats-ci.yml +++ b/.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 diff --git a/docs/documentation/source/ats_demos b/docs/documentation/source/ats_demos index bc3584da2..b39f15c8a 160000 --- a/docs/documentation/source/ats_demos +++ b/docs/documentation/source/ats_demos @@ -1 +1 @@ -Subproject commit bc3584da2acd9e5a13d651958225328c4e24fccf +Subproject commit b39f15c8a052ff142725cef0e80fd3883f45e096 diff --git a/src/constitutive_relations/column_integrators/EvaluatorColumnIntegrator.hh b/src/constitutive_relations/column_integrators/EvaluatorColumnIntegrator.hh index 17d73438b..60173e2e9 100644 --- a/src/constitutive_relations/column_integrators/EvaluatorColumnIntegrator.hh +++ b/src/constitutive_relations/column_integrators/EvaluatorColumnIntegrator.hh @@ -138,10 +138,12 @@ EvaluatorColumnIntegrator::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]; } + } diff --git a/src/constitutive_relations/column_integrators/thaw_depth_evaluator.cc b/src/constitutive_relations/column_integrators/thaw_depth_evaluator.cc index 1a62c87c7..e01e94e3a 100644 --- a/src/constitutive_relations/column_integrators/thaw_depth_evaluator.cc +++ b/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}); } diff --git a/src/constitutive_relations/eos/eos_constant.cc b/src/constitutive_relations/eos/eos_constant.cc index 15aaee143..4de114a21 100644 --- a/src/constitutive_relations/eos/eos_constant.cc +++ b/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("molar mass [kg/mol]"); + if (eos_plist_.isParameter("molar mass [kg mol^-1]")) { + M_ = eos_plist_.get("molar mass [kg mol^-1]"); } else { - M_ = eos_plist_.get("molar mass [g/mol]", 18.0153) * 1.e-3; + M_ = eos_plist_.get("molar mass [g mol^-1]", 18.0153) * 1.e-3; } - if (eos_plist_.isParameter("density [mol/m^3]")) { - rho_ = eos_plist_.get("density [mol/m^3]") * M_; + if (eos_plist_.isParameter("density [mol m^-3]")) { + rho_ = eos_plist_.get("density [mol m^-3]") * M_; } else { - rho_ = eos_plist_.get("density [kg/m^3]"); + rho_ = eos_plist_.get("density [kg m^-3]"); } }; diff --git a/src/constitutive_relations/generic_evaluators/AdditiveEvaluator.cc b/src/constitutive_relations/generic_evaluators/AdditiveEvaluator.cc index 0483690e4..145a1cd2e 100644 --- a/src/constitutive_relations/generic_evaluators/AdditiveEvaluator.cc +++ b/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(varname+" coefficient"); - } else if (plist.isParameter(dep.first+" coefficient")) { - coefs_[varname] = plist.get(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(variable_tag+" coefficient"); + } else if (plist.isParameter(variable+" coefficient")) { + coefs_[variable_tag] = plist.get(variable+" coefficient"); + } else if (plist.isParameter(varname+" coefficient")) { + coefs_[variable_tag] = plist.get(varname+" coefficient"); } else { - coefs_[varname] = 1.0; + coefs_[variable_tag] = 1.0; } } shift_ = plist.get("constant shift", 0.); + positive_ = plist.get("enforce positivity", false); } @@ -51,6 +56,13 @@ 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 @@ -58,6 +70,19 @@ AdditiveEvaluator::EvaluatePartialDerivative_(const State& S, const Key& wrt_key, const Tag& wrt_tag, const std::vector& result) { result[0]->PutScalar(coefs_[Keys::getKey(wrt_key, wrt_tag)]); + + if (positive_) { + const auto& value = S.Get(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.; + } + } + } + } } diff --git a/src/constitutive_relations/generic_evaluators/AdditiveEvaluator.hh b/src/constitutive_relations/generic_evaluators/AdditiveEvaluator.hh index f2e820979..6be3b0de2 100644 --- a/src/constitutive_relations/generic_evaluators/AdditiveEvaluator.hh +++ b/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 coefs_; double shift_; + bool positive_; private: static Utils::RegisteredFactory factory_; diff --git a/src/constitutive_relations/generic_evaluators/MultiplicativeEvaluator.cc b/src/constitutive_relations/generic_evaluators/MultiplicativeEvaluator.cc index ecb722166..462a97e3c 100644 --- a/src/constitutive_relations/generic_evaluators/MultiplicativeEvaluator.cc +++ b/src/constitutive_relations/generic_evaluators/MultiplicativeEvaluator.cc @@ -96,8 +96,12 @@ MultiplicativeEvaluator::EvaluatePartialDerivative_(const State& S, } } if (positive_) { + const auto& value_c = *S.Get(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.; + } } } } diff --git a/src/constitutive_relations/generic_evaluators/SubgridAggregateEvaluator.cc b/src/constitutive_relations/generic_evaluators/SubgridAggregateEvaluator.cc index 846efd5c5..faf9b52d0 100644 --- a/src/constitutive_relations/generic_evaluators/SubgridAggregateEvaluator.cc +++ b/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(dependencies_.front().first, dependencies_.front().second); + if (dep_fac.HasComponent("cell")) { + S.Require(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(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,}); + } } diff --git a/src/constitutive_relations/generic_evaluators/SubgridAggregateEvaluator.hh b/src/constitutive_relations/generic_evaluators/SubgridAggregateEvaluator.hh index 05a231f1e..65478d2cc 100644 --- a/src/constitutive_relations/generic_evaluators/SubgridAggregateEvaluator.hh +++ b/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; diff --git a/src/constitutive_relations/generic_evaluators/SubgridDisaggregateEvaluator.cc b/src/constitutive_relations/generic_evaluators/SubgridDisaggregateEvaluator.cc index 88a72b56d..92e5aebfb 100644 --- a/src/constitutive_relations/generic_evaluators/SubgridDisaggregateEvaluator.cc +++ b/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(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); + } } diff --git a/src/executables/ats_mesh_factory.cc b/src/executables/ats_mesh_factory.cc index 75cdf40a7..ac8a46e50 100644 --- a/src/executables/ats_mesh_factory.cc +++ b/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("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(); } diff --git a/src/executables/ats_mesh_factory.hh b/src/executables/ats_mesh_factory.hh index bfe3a24d8..8fa92077c 100644 --- a/src/executables/ats_mesh_factory.hh +++ b/src/executables/ats_mesh_factory.hh @@ -85,9 +85,9 @@ Example: - - - + + + diff --git a/src/executables/coordinator.cc b/src/executables/coordinator.cc index b8f4b8261..ee39ca449 100644 --- a/src/executables/coordinator.cc +++ b/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 diff --git a/src/executables/main.cc b/src/executables/main.cc index 5bf9ff418..8a9de36e1 100644 --- a/src/executables/main.cc +++ b/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()) { diff --git a/src/pks/deform/volumetric_deformation.cc b/src/pks/deform/volumetric_deformation.cc index 231778cff..1b6c3262f 100644 --- a/src/pks/deform/volumetric_deformation.cc +++ b/src/pks/deform/volumetric_deformation.cc @@ -29,10 +29,10 @@ VolumetricDeformation::VolumetricDeformation(Teuchos::ParameterList& pk_tree, const Teuchos::RCP& 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("max time step [s]", std::numeric_limits::max()); - dt_ = dt_max_; // The deformation mode describes how to calculate new cell volume from a // provided function and the old cell volume. diff --git a/src/pks/deform/volumetric_deformation.hh b/src/pks/deform/volumetric_deformation.hh index 2cb0ce53a..fbf0ac77f 100644 --- a/src/pks/deform/volumetric_deformation.hh +++ b/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_; diff --git a/src/pks/energy/energy_bc_factory.hh b/src/pks/energy/energy_bc_factory.hh index 42336e41f..56c3c97a7 100644 --- a/src/pks/energy/energy_bc_factory.hh +++ b/src/pks/energy/energy_bc_factory.hh @@ -1,13 +1,58 @@ -/* -*- mode: c++; indent-tabs-mode: nil -*- */ +/* + 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. -/* ------------------------------------------------------------------------- -ATS + Authors: Ethan Coon (ecoon@lanl.gov) +*/ -License: see $ATS_DIR/COPYRIGHT -Author: Ethan Coon +/*! +Energy boundary conditions must follow the general format shown in +BoundaryConditions_. Energy-specific conditions implemented include: + +Dirichlet (temperature) boundary conditions +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +Provide temperature data in units of [K]. + +Example: + +.. code-block:: xml + + + + + + + + + + + + + +Neumann (diffusive energy flux) boundary conditions +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Note that for an energy equation, there are two mechanisms by which energy can +be fluxed into the domain: through diffusion and through advection. This +boundary condition sets the diffusive flux of energy, and allows the advective +flux to be whatever it may be. Frequently this is used in combination with +boundaries where water is expected to be advected out of the domain, and we +wish to allow the energy of that water to be advected away with it, but wish to +independently specify diffusive fluxes. This can also be used in cases where +the mass flux is prescribed to be zero (e.g. bottom boundaries, where this +might be the geothermal gradient). + +Units are in [W m^-2]. + +Example: + +.. code-block:: xml + + @@ -20,8 +65,41 @@ Author: Ethan Coon +Note that another commonly implemented boundary condition is one where the +diffusive flux is prescribed, and also the temperature of incoming water is +prescribed. This is not currently implemented, but would be straightforward to +do so if requested. + + +Neumann (total energy flux) boundary conditions +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +This boundary condition sets the total flux of energy, from both advection and +diffusion. This is not used all that often in real applications, but is common +for benchmarks or other testing. + +Units are in [W m^-2]. + +Example: + +.. code-block:: xml + + + + + + + + + + + + + + -------------------------------------------------------------------------- */ + +*/ #ifndef AMANZI_ENERGY_BC_FACTORY_HH_ #define AMANZI_ENERGY_BC_FACTORY_HH_ diff --git a/src/pks/energy/energy_surface_ice.cc b/src/pks/energy/energy_surface_ice.cc index 566e8cdae..213182240 100644 --- a/src/pks/energy/energy_surface_ice.cc +++ b/src/pks/energy/energy_surface_ice.cc @@ -225,7 +225,7 @@ void EnergySurfaceIce::AddSources_(const Tag& tag, const Teuchos::Ptros_OK(Teuchos::VERB_EXTREME)) { *vo_->os() << "Adding advection to subsurface" << std::endl; } - db_->WriteVector("res (src post surf-subsurf adv)", g, false); + db_->WriteVector("res (s-s adv src)", g, false); } // -- conduction source @@ -239,9 +239,8 @@ void EnergySurfaceIce::AddSources_(const Tag& tag, const Teuchos::PtrWriteVector("res (src post surf-subsurf diff)", g, false); + db_->WriteVector("res (s-s adv src)", g, false); } - } diff --git a/src/pks/flow/flow_bc_factory.hh b/src/pks/flow/flow_bc_factory.hh index 126c93db4..ed5f2a7d1 100644 --- a/src/pks/flow/flow_bc_factory.hh +++ b/src/pks/flow/flow_bc_factory.hh @@ -1,6 +1,3 @@ -/* -*- mode: c++; indent-tabs-mode: nil -*- */ -// Boundary conditions packaged for use by flow PKs. - /* ATS is released under the three-clause BSD License. The terms of use and "as is" disclaimer for this license are diff --git a/src/pks/flow/overland_pressure_pk.cc b/src/pks/flow/overland_pressure_pk.cc index 0774ef804..ef579bbc1 100644 --- a/src/pks/flow/overland_pressure_pk.cc +++ b/src/pks/flow/overland_pressure_pk.cc @@ -460,6 +460,7 @@ void OverlandPressureFlow::Initialize() // Initialize BC values ComputeBoundaryConditions_(tag_next_); + ChangedSolution(tag_next_); // Set extra fields as initialized -- these don't currently have evaluators. S_->GetW(uw_cond_key_, tag_next_, name_).PutScalar(0.0); @@ -472,7 +473,9 @@ void OverlandPressureFlow::Initialize() S_->GetW(flux_dir_key_, tag_next_, name_).PutScalar(0.); S_->GetRecordW(flux_dir_key_, tag_next_, name_).set_initialized(); + S_->GetW(velocity_key_, Tags::NEXT, name_).PutScalar(0.); + changedEvaluatorPrimary(velocity_key_, Tags::NEXT, *S_); S_->GetRecordW(velocity_key_, Tags::NEXT, name_).set_initialized(); }; diff --git a/src/pks/flow/richards_pk.cc b/src/pks/flow/richards_pk.cc index 704bd8eaa..f37bdbd46 100644 --- a/src/pks/flow/richards_pk.cc +++ b/src/pks/flow/richards_pk.cc @@ -677,13 +677,12 @@ Richards::ValidStep() ->ViewComponent("cell",false); Epetra_MultiVector dsl(sl_new); dsl.Update(-1., sl_old, 1.); - double change = 0.; - dsl.NormInf(&change); + auto change = maxValLoc(*dsl(0)); - if (change > sat_change_limit_) { + if (change.value > sat_change_limit_) { if (vo_->os_OK(Teuchos::VERB_LOW)) *vo_->os() << "Invalid time step, max sl change=" - << change << " > limit=" << sat_change_limit_ << std::endl; + << change.value << " > limit=" << sat_change_limit_ << " at cell GID " << change.gid << std::endl; return false; } } @@ -695,13 +694,12 @@ Richards::ValidStep() ->ViewComponent("cell",false); Epetra_MultiVector dsi(si_new); dsi.Update(-1., si_old, 1.); - double change = 0.; - dsi.NormInf(&change); + auto change = maxValLoc(*dsi(0)); - if (change > sat_ice_change_limit_) { + if (change.value > sat_ice_change_limit_) { if (vo_->os_OK(Teuchos::VERB_LOW)) *vo_->os() << "Invalid time step, max si change=" - << change << " > limit=" << sat_ice_change_limit_ << std::endl; + << change.value << " > limit=" << sat_ice_change_limit_ << " at cell GID " << change.gid << std::endl; return false; } } diff --git a/src/pks/mpc/mpc_permafrost_split_flux.hh b/src/pks/mpc/mpc_permafrost_split_flux.hh index 509f8bb2f..1923eff04 100644 --- a/src/pks/mpc/mpc_permafrost_split_flux.hh +++ b/src/pks/mpc/mpc_permafrost_split_flux.hh @@ -101,14 +101,22 @@ class MPCPermafrostSplitFlux : public MPCSubcycled { void CopyStarToPrimary_Standard_Hybrid_(); Tag get_ds_tag_next_(const std::string& subdomain) { - if (subcycling_[1]) - return Tag{Keys::cleanName(tags_[1].second.get() + "_" + Keys::getDomainSetIndex(subdomain))}; - else return Tag{tags_[1].second}; + if (subcycling_[1]) { + AMANZI_ASSERT(Keys::starts_with(subdomain, "surface_")); + return Tag(Keys::getKey(subdomain.substr(std::string("surface_").size(), + std::string::npos), tags_[1].second.get())); + } else { + return Tag{tags_[1].second}; + } } Tag get_ds_tag_current_(const std::string& subdomain) { - if (subcycling_[1]) - return Tag{Keys::cleanName(tags_[1].first.get() + "_" + Keys::getDomainSetIndex(subdomain))}; - else return Tag{tags_[1].first}; + if (subcycling_[1]) { + AMANZI_ASSERT(Keys::starts_with(subdomain, "surface_")); + return Tag(Keys::getKey(subdomain.substr(std::string("surface_").size(), + std::string::npos), tags_[1].first.get())); + } else { + return Tag{tags_[1].first}; + } } diff --git a/src/pks/mpc/mpc_subcycled.cc b/src/pks/mpc/mpc_subcycled.cc index 6c03acd88..aa5b1b165 100644 --- a/src/pks/mpc/mpc_subcycled.cc +++ b/src/pks/mpc/mpc_subcycled.cc @@ -41,8 +41,7 @@ MPCSubcycled::MPCSubcycled(Teuchos::ParameterList& pk_tree, dts_.resize(sub_pks_.size(), -1); // min dt allowed in subcycling - max_dt_ = plist_->get("subcycling target time step [s]", -1); - min_dt_ = plist_->get("minimum subcycled time step [s]", 1.e-4); + target_dt_ = plist_->get("subcycling target time step [s]", -1); } @@ -101,7 +100,7 @@ MPCSubcycled::Initialize() double MPCSubcycled::get_dt() { double dt = std::numeric_limits::max(); - if (max_dt_ > 0) dt = max_dt_; + if (target_dt_ > 0) dt = target_dt_; int i = 0; for (auto& pk : sub_pks_) { @@ -184,12 +183,6 @@ MPCSubcycled::AdvanceStep_i_(std::size_t i, double t_old, double t_new, bool rei if (vo_->os_OK(Teuchos::VERB_EXTREME)) *vo_->os() << " success, new timestep is " << dt_inner << std::endl; } - - if (dt_inner < min_dt_) { - Errors::Message msg; - msg << "SubPK " << sub_pks_[i]->name() << " crashing timestep in subcycling: dt = " << dt_inner; - Exceptions::amanzi_throw(msg); - } } } else { diff --git a/src/pks/mpc/mpc_subcycled.hh b/src/pks/mpc/mpc_subcycled.hh index cb409db49..a719072a4 100644 --- a/src/pks/mpc/mpc_subcycled.hh +++ b/src/pks/mpc/mpc_subcycled.hh @@ -55,7 +55,7 @@ public: bool AdvanceStep_i_(std::size_t i, double t_old, double t_new, bool reinit); Teuchos::Array subcycling_; - double dt_, min_dt_, max_dt_; + double dt_, target_dt_; std::vector dts_; std::vector> tags_; diff --git a/src/pks/mpc/mpc_weak_subdomain.cc b/src/pks/mpc/mpc_weak_subdomain.cc index a62e23b5d..8926b4de3 100644 --- a/src/pks/mpc/mpc_weak_subdomain.cc +++ b/src/pks/mpc/mpc_weak_subdomain.cc @@ -44,7 +44,6 @@ MPCWeakSubdomain::MPCWeakSubdomain(Teuchos::ParameterList& FElist, subcycled_ = plist_->template get("subcycle", false); if (subcycled_) { subcycled_target_dt_ = plist_->template get("subcycling target time step [s]"); - subcycled_min_dt_ = plist_->template get("minimum subcycled time step [s]", 1.e-4); } }; @@ -174,65 +173,79 @@ MPCWeakSubdomain::AdvanceStep_Subcycled_(double t_old, double t_new, bool reinit int my_pid = solution_->Comm()->MyPID(); int i = 0; + int n_throw = 0; + std::string throw_msg; for (const auto& subdomain : ds) { - if (vo_->os_OK(Teuchos::VERB_EXTREME)) - *vo_->os() << "Beginning subcyling on pk \"" << sub_pks_[i]->name() << "\"" << std::endl; - - bool done = false; + double dt_inner = -1; double t_inner = t_old; - - Tag tag_subcycle_current = get_ds_tag_current_(subdomain); - Tag tag_subcycle_next = get_ds_tag_next_(subdomain); - - S_->set_time(tag_subcycle_current, t_old); - while (!done) { - double dt_inner = std::min(sub_pks_[i]->get_dt(), t_new - t_inner); - S_->Assign("dt", tag_subcycle_next, name(), dt_inner); - S_->set_time(tag_subcycle_next, t_inner + dt_inner); - bool fail_inner = sub_pks_[i]->AdvanceStep(t_inner, t_inner+dt_inner, false); + try { // must catch non-collective throws for TimeStepCrash if (vo_->os_OK(Teuchos::VERB_EXTREME)) - *vo_->os() << " step failed? " << fail_inner << std::endl; - bool valid_inner = sub_pks_[i]->ValidStep(); - if (vo_->os_OK(Teuchos::VERB_EXTREME)) { - *vo_->os() << " step valid? " << valid_inner << std::endl; - } + *vo_->os() << "Beginning subcyling on pk \"" << sub_pks_[i]->name() << "\"" << std::endl; - if (fail_inner || !valid_inner) { - sub_pks_[i]->FailStep(t_old, t_new, tag_subcycle_next); - - dt_inner = sub_pks_[i]->get_dt(); - S_->set_time(tag_subcycle_next, S_->get_time(tag_subcycle_current)); + bool done = false; + Tag tag_subcycle_current = get_ds_tag_current_(subdomain); + Tag tag_subcycle_next = get_ds_tag_next_(subdomain); + S_->set_time(tag_subcycle_current, t_old); + while (!done) { + dt_inner = std::min(sub_pks_[i]->get_dt(), t_new - t_inner); + S_->Assign("dt", tag_subcycle_next, name(), dt_inner); + S_->set_time(tag_subcycle_next, t_inner + dt_inner); + bool fail_inner = sub_pks_[i]->AdvanceStep(t_inner, t_inner+dt_inner, false); if (vo_->os_OK(Teuchos::VERB_EXTREME)) - *vo_->os() << " failed, new timestep is " << dt_inner << std::endl; + *vo_->os() << " step failed? " << fail_inner << std::endl; + bool valid_inner = sub_pks_[i]->ValidStep(); + if (vo_->os_OK(Teuchos::VERB_EXTREME)) { + *vo_->os() << " step valid? " << valid_inner << std::endl; + } - } else { - sub_pks_[i]->CommitStep(t_inner, t_inner + dt_inner, tag_subcycle_next); - t_inner += dt_inner; - if (std::abs(t_new - t_inner) < 1.e-10) done = true; + if (fail_inner || !valid_inner) { + sub_pks_[i]->FailStep(t_old, t_new, tag_subcycle_next); - S_->set_time(tag_subcycle_current, S_->get_time(tag_subcycle_next)); - S_->advance_cycle(tag_subcycle_next); + dt_inner = sub_pks_[i]->get_dt(); + S_->set_time(tag_subcycle_next, S_->get_time(tag_subcycle_current)); - dt_inner = sub_pks_[i]->get_dt(); - if (vo_->os_OK(Teuchos::VERB_EXTREME)) - *vo_->os() << " success, new timestep is " << dt_inner << std::endl; - } + if (vo_->os_OK(Teuchos::VERB_EXTREME)) + *vo_->os() << " failed, new timestep is " << dt_inner << std::endl; + + } else { + sub_pks_[i]->CommitStep(t_inner, t_inner + dt_inner, tag_subcycle_next); + t_inner += dt_inner; + if (std::abs(t_new - t_inner) < 1.e-10) done = true; - if (dt_inner < subcycled_min_dt_) { - Errors::Message msg; - msg << "Subdomain " << subdomain << " on PID " << my_pid << " crashing timestep in subcycling: dt = " << dt_inner; - Exceptions::amanzi_throw(msg); + S_->set_time(tag_subcycle_current, S_->get_time(tag_subcycle_next)); + S_->advance_cycle(tag_subcycle_next); + + dt_inner = sub_pks_[i]->get_dt(); + if (vo_->os_OK(Teuchos::VERB_EXTREME)) + *vo_->os() << " success, new timestep is " << dt_inner << std::endl; + } } + i++; + } catch(Errors::TimeStepCrash& e) { + n_throw++; + throw_msg = e.what(); + break; } - i++; } + // check for any other ranks throwing and, if so, throw ourselves so that all procs throw + int n_throw_g = 0; + comm_->SumAll(&n_throw, &n_throw_g, 1); + if (n_throw > 0) { + // inject more information into the crash message + Errors::TimeStepCrash msg; + msg << throw_msg << " on rank " << comm_->MyPID() << " of " << comm_->NumProc(); + Exceptions::amanzi_throw(msg); + } else if (n_throw_g > 0) { + Errors::TimeStepCrash msg; + msg << "TimeStepCrash on another rank: nprocs failed = " << n_throw_g; + Exceptions::amanzi_throw(msg); + } return false; } - void MPCWeakSubdomain::CommitStep(double t_old, double t_new, const Tag& tag_next) { @@ -291,7 +304,9 @@ MPCWeakSubdomain::init_() PKFactory pk_factory; for (const auto& subdomain : *ds) { auto mesh = S_->GetMesh(subdomain); - // create the solution vector + // create the solution vector, noting that these are created on the + // communicator associated with the mesh of the subdomain, which may differ + // from the coupler's comm. Teuchos::RCP pk_soln = Teuchos::rcp(new TreeVector(mesh->get_comm())); solution_->PushBack(pk_soln); diff --git a/src/pks/mpc/mpc_weak_subdomain.hh b/src/pks/mpc/mpc_weak_subdomain.hh index d9c52979a..bad576c1a 100644 --- a/src/pks/mpc/mpc_weak_subdomain.hh +++ b/src/pks/mpc/mpc_weak_subdomain.hh @@ -49,19 +49,18 @@ class MPCWeakSubdomain : public MPC { Tag get_ds_tag_next_(const std::string& subdomain) { if (subcycled_) - return Tag(Keys::cleanName(tag_next_.get() + "_" + Keys::getDomainSetIndex(subdomain))); + return Tag(Keys::getKey(subdomain, tag_next_.get())); else return tag_next_; } Tag get_ds_tag_current_(const std::string& subdomain) { if (subcycled_) - return Tag(Keys::cleanName(tag_current_.get() + "_" + Keys::getDomainSetIndex(subdomain))); + return Tag(Keys::getKey(subdomain, tag_current_.get())); else return tag_current_; } Comm_ptr_type comm_; bool subcycled_; double subcycled_target_dt_; - double subcycled_min_dt_; double cycle_dt_; Key ds_name_; diff --git a/src/pks/pk_bdf_default.cc b/src/pks/pk_bdf_default.cc index 69c5bdcce..60c0a1c62 100644 --- a/src/pks/pk_bdf_default.cc +++ b/src/pks/pk_bdf_default.cc @@ -38,7 +38,7 @@ void PK_BDF_Default::Setup() } // require data for checkpointing timestep size - S_->Require("dt", Tag(name_), name_); + S_->Require("dt_internal", Tag(name_), name_); } }; @@ -54,14 +54,14 @@ void PK_BDF_Default::Initialize() // Note, this is done here and not in setup because solution is not ready in setup Teuchos::ParameterList& bdf_plist = plist_->sublist("time integrator"); bdf_plist.sublist("verbose object").setParametersNotAlreadySet(plist_->sublist("verbose object")); - bdf_plist.sublist("verbose object").set("name", name()); + bdf_plist.sublist("verbose object").set("name", name()+"_TI"); time_stepper_ = Teuchos::rcp(new BDF1_TI(*this, bdf_plist, solution_, S_)); double dt_init = time_stepper_->initial_timestep(); - S_->Assign("dt", Tag(name_), name_, dt_init); - S_->GetRecordW("dt", Tag(name_), name_).set_initialized(); + S_->Assign("dt_internal", Tag(name_), name_, dt_init); + S_->GetRecordW("dt_internal", Tag(name_), name_).set_initialized(); // -- initialize continuation parameter if needed. if (S_->HasRecord("continuation_parameter", Tag(name_))) { @@ -78,30 +78,20 @@ void PK_BDF_Default::Initialize() } }; -void PK_BDF_Default::ResetTimeStepper(double time) -{ - // -- initialize time derivative - auto solution_dot = Teuchos::rcp(new TreeVector(*solution_)); - solution_dot->PutScalar(0.0); - - // -- set initial state - time_stepper_->SetInitialState(time, solution_, solution_dot); - return; -} // ----------------------------------------------------------------------------- // Initialization of timestepper. // ----------------------------------------------------------------------------- double PK_BDF_Default::get_dt() { if (!strongly_coupled_) - return S_->Get("dt", Tag(name_)); + return S_->Get("dt_internal", Tag(name_)); else return -1.; } void PK_BDF_Default::set_dt(double dt) { if (!strongly_coupled_) - S_->Assign("dt", Tag(name_), name_, dt); + S_->Assign("dt_internal", Tag(name_), name_, dt); } // -- Commit any secondary (dependent) variables. @@ -109,12 +99,8 @@ void PK_BDF_Default::CommitStep(double t_old, double t_new, const Tag& tag) { if (tag == tag_next_) { double dt = t_new - t_old; - if (time_stepper_ != Teuchos::null) { - if (dt <= 0) { - ResetTimeStepper(t_old); - } else { - time_stepper_->CommitSolution(dt, solution_, true); - } + if (time_stepper_ != Teuchos::null && dt > 0) { + time_stepper_->CommitSolution(dt, solution_, true); } } } @@ -141,7 +127,7 @@ bool PK_BDF_Default::AdvanceStep(double t_old, double t_new, bool reinit) // -- 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("dt", Tag(name_)); + double dt_internal = S_->Get("dt_internal", Tag(name_)); // Note, the fact that this is triggering an assertion on old old runs // indicates that there may be a long-standing bug in TimeStepController. @@ -150,37 +136,50 @@ bool PK_BDF_Default::AdvanceStep(double t_old, double t_new, bool reinit) //AMANZI_ASSERT(dt <= dt_internal + 1.e-8); // roundoff double dt_solver = -1; - bool fail = time_stepper_->TimeStep(dt, dt_solver, solution_); - if (!fail) { - // check step validity - bool valid = ValidStep(); - if (valid) { - if (vo_->os_OK(Teuchos::VERB_LOW)) - *vo_->os() << "successful advance" << std::endl; - // update the timestep size - if (dt_solver < dt_internal && dt_solver >= dt) { - // We took a smaller step than we recommended, and it worked fine (not - // suprisingly). Likely this was due to constraints from other PKs or - // vis. Do not reduce our recommendation. + bool fail = false; + try { + fail = time_stepper_->TimeStep(dt, dt_solver, solution_); + + if (!fail) { + // check step validity + bool valid = ValidStep(); + if (valid) { + if (vo_->os_OK(Teuchos::VERB_LOW)) + *vo_->os() << "successful advance" << std::endl; + // update the timestep size + if (dt_solver < dt_internal && dt_solver >= dt) { + // We took a smaller step than we recommended, and it worked fine (not + // suprisingly). Likely this was due to constraints from other PKs or + // vis. Do not reduce our recommendation. + } else { + dt_internal = dt_solver; + } } else { - dt_internal = dt_solver; + if (vo_->os_OK(Teuchos::VERB_LOW)) + *vo_->os() << "successful advance, but not valid" << std::endl; + time_stepper_->CommitSolution(dt_internal, solution_, valid); + dt_internal = 0.5 * dt_internal; + // when including Valid here, make fail = true refs #110 } } else { if (vo_->os_OK(Teuchos::VERB_LOW)) - *vo_->os() << "successful advance, but not valid" << std::endl; - time_stepper_->CommitSolution(dt_internal, solution_, valid); - dt_internal = 0.5 * dt_internal; - // when including Valid here, make fail = true refs #110 + *vo_->os() << "unsuccessful advance" << std::endl; + // take the decreased timestep size + dt_internal = dt_solver; } - } else { - if (vo_->os_OK(Teuchos::VERB_LOW)) - *vo_->os() << "unsuccessful advance" << std::endl; - // take the decreased timestep size - dt_internal = dt_solver; - } - S_->Assign("dt", Tag(name_), name_, dt_internal); + S_->Assign("dt_internal", Tag(name_), name_, dt_internal); + } catch(Errors::TimeStepCrash& e) { + // inject more information into the crash message + std::stringstream msg_str; + msg_str << "TimeStepCrash in PK: \"" << name() << "\"" << std::endl + << " at t = " << t_old << " with dt = " << dt << std::endl + << " error message: " << std::endl << std::endl + << e.what() << std::endl << std::endl; + Errors::TimeStepCrash msg(msg_str.str()); + Exceptions::amanzi_throw(msg); + } return fail; }; diff --git a/src/pks/pk_bdf_default.hh b/src/pks/pk_bdf_default.hh index cd23a9e73..33cfdbdb2 100644 --- a/src/pks/pk_bdf_default.hh +++ b/src/pks/pk_bdf_default.hh @@ -102,11 +102,8 @@ class PK_BDF_Default : public PK_BDF { return AmanziSolvers::FnBaseDefs::CORRECTION_NOT_MODIFIED; } - virtual void ResetTimeStepper(double time); - - // experimental approach -- calling this indicates that the time - // integration scheme is changing the value of the solution in - // state. + // Calling this indicates that the time integration scheme is changing the + // value of the solution in state. virtual void ChangedSolution() override = 0; virtual void ChangedSolution(const Tag& tag) = 0; diff --git a/src/pks/pk_helpers.cc b/src/pks/pk_helpers.cc index 0e5142114..27a794836 100644 --- a/src/pks/pk_helpers.cc +++ b/src/pks/pk_helpers.cc @@ -326,5 +326,24 @@ copyVectorToMeshCoordinates(const CompositeVector& vec, mesh.deform(node_ids, new_positions); } +int commMaxValLoc(const Comm_type& comm, const ValLoc& local, ValLoc& global) { + MpiComm_type const* mpi_comm = dynamic_cast(&comm); + const MPI_Comm& mpi_comm_raw = mpi_comm->Comm(); + return MPI_Allreduce(&local, &global, 1, MPI_DOUBLE_INT, MPI_MAXLOC, mpi_comm_raw); +} + +ValLoc maxValLoc(const Epetra_Vector& vec) { + ValLoc local{0.,0}; + for (int i=0; i!=vec.MyLength(); ++i) { + if (vec[i] > local.value) { + local.value = vec[i]; + local.gid = vec.Map().GID(i); + } + } + ValLoc global{0.,0}; + int ierr = commMaxValLoc(vec.Comm(), local, global); + AMANZI_ASSERT(!ierr); + return global; +} } // namespace Amanzi diff --git a/src/pks/pk_helpers.hh b/src/pks/pk_helpers.hh index 8ceb5fb69..367745eaa 100644 --- a/src/pks/pk_helpers.hh +++ b/src/pks/pk_helpers.hh @@ -128,4 +128,14 @@ copyVectorToMeshCoordinates(const CompositeVector& vec, AmanziMesh::Mesh& mesh); + +// Compute pairs of value + location +typedef struct ValLoc { + double value; + AmanziMesh::Entity_ID gid; +} ENorm_t; + +int commMaxValLoc(const Comm_type& comm, const ValLoc& local, ValLoc& global); +ValLoc maxValLoc(const Epetra_Vector& vec); + } // namespace Amanzi diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_threecomponent_evaluator.cc b/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_threecomponent_evaluator.cc index 5e70fa294..0129a8db6 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_threecomponent_evaluator.cc +++ b/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_threecomponent_evaluator.cc @@ -136,6 +136,15 @@ AreaFractionsThreeComponentEvaluator::Evaluate_(const State& S, } +// custom EC used to set subfield names +void +AreaFractionsThreeComponentEvaluator::EnsureCompatibility_Structure_(State& S) +{ + S.GetRecordSetW(my_keys_.front().first).set_subfieldnames( + {"bare", "water", "snow"}); +} + + void AreaFractionsThreeComponentEvaluator::EnsureCompatibility_ToDeps_(State& S) { diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_threecomponent_evaluator.hh b/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_threecomponent_evaluator.hh index 688319f13..78bfa1fdd 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_threecomponent_evaluator.hh +++ b/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_threecomponent_evaluator.hh @@ -62,14 +62,17 @@ class AreaFractionsThreeComponentEvaluator : public EvaluatorSecondaryMonotypeCV } protected: - virtual void EnsureCompatibility_ToDeps_(State& S) override; + // custom EC used to set subfield names + void EnsureCompatibility_Structure_(State& S) override; + + // custom EC used because deps have 1 component not 3 + void EnsureCompatibility_ToDeps_(State& S) override; // Required methods from EvaluatorSecondaryMonotypeCV virtual void Evaluate_(const State& S, const std::vector& result) override; virtual void EvaluatePartialDerivative_(const State& S, const Key& wrt_key, const Tag& wrt_tag, const std::vector& result) override { - result[0]->PutScalar(0.); //Errors::Message msg("NotImplemented: AreaFractionsThreeComponentEvaluator currently does not provide derivatives."); //Exceptions::amanzi_throw(msg); diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_threecomponent_microtopography_evaluator.cc b/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_threecomponent_microtopography_evaluator.cc index 51e1c9326..401c0c4af 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_threecomponent_microtopography_evaluator.cc +++ b/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_threecomponent_microtopography_evaluator.cc @@ -124,6 +124,16 @@ AreaFractionsThreeComponentMicrotopographyEvaluator::Evaluate_(const State& S, } } + +// custom EC used to set subfield names +void +AreaFractionsThreeComponentMicrotopographyEvaluator::EnsureCompatibility_Structure_(State& S) +{ + S.GetRecordSetW(my_keys_.front().first).set_subfieldnames( + {"bare", "water", "snow"}); +} + + void AreaFractionsThreeComponentMicrotopographyEvaluator::EnsureCompatibility_ToDeps_(State& S) { diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_threecomponent_microtopography_evaluator.hh b/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_threecomponent_microtopography_evaluator.hh index 9aa87c316..43e5480f5 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_threecomponent_microtopography_evaluator.hh +++ b/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_threecomponent_microtopography_evaluator.hh @@ -62,7 +62,11 @@ class AreaFractionsThreeComponentMicrotopographyEvaluator : public EvaluatorSeco } protected: - virtual void EnsureCompatibility_ToDeps_(State& S) override; + // custom EC used to set subfield names + void EnsureCompatibility_Structure_(State& S) override; + + // custom EC used because deps have 1 component not 3 + void EnsureCompatibility_ToDeps_(State& S) override; // Required methods from EvaluatorSecondaryMonotypeCV virtual void Evaluate_(const State& S, diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_twocomponent_evaluator.cc b/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_twocomponent_evaluator.cc index e24764efa..d90ad94d1 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_twocomponent_evaluator.cc +++ b/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_twocomponent_evaluator.cc @@ -96,6 +96,15 @@ AreaFractionsTwoComponentEvaluator::EvaluatePartialDerivative_(const State& S, } +// custom EC used to set subfield names +void +AreaFractionsTwoComponentEvaluator::EnsureCompatibility_Structure_(State& S) +{ + S.GetRecordSetW(my_keys_.front().first).set_subfieldnames( + {"bare/water", "snow"}); +} + + void AreaFractionsTwoComponentEvaluator::EnsureCompatibility_ToDeps_(State& S) { diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_twocomponent_evaluator.hh b/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_twocomponent_evaluator.hh index 7b63b6ed0..253ee1dc6 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_twocomponent_evaluator.hh +++ b/src/pks/surface_balance/constitutive_relations/land_cover/area_fractions_twocomponent_evaluator.hh @@ -60,7 +60,12 @@ class AreaFractionsTwoComponentEvaluator : public EvaluatorSecondaryMonotypeCV { } protected: - virtual void EnsureCompatibility_ToDeps_(State& S) override; + + // custom EC used to set subfield names + void EnsureCompatibility_Structure_(State& S) override; + + // custom EC used because deps have 1 component not 2 + void EnsureCompatibility_ToDeps_(State& S) override; // Required methods from EvaluatorSecondaryMonotypeCV virtual void Evaluate_(const State& S, diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_defs.hh b/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_defs.hh index 14b928eb0..cc3d4901f 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_defs.hh +++ b/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_defs.hh @@ -227,7 +227,6 @@ struct FluxBalance { struct SurfaceParams { double a_tundra, a_water, a_ice; // albedos double e_snow, e_tundra, e_water, e_ice; // emissivities - double Zsmooth, Zrough; // roughness coefs SurfaceParams() : a_tundra(0.135), // [-] Grenfell and Perovich, (2004) @@ -237,9 +236,7 @@ struct SurfaceParams { e_tundra(0.92), // [-] emissivity for tundra, From P. ReVelle // (Thesis), Ling & Zhang, 2004 e_water(0.995), // [-] emissivity of water, EngineeringToolbox.com - e_ice(0.98), // [-] emissivity of ice, EngineeringToolbox.com - Zsmooth(0.005), // [m]? roughness coef of smooth - Zrough(0.03) {} // [m]? roughness coef of rough + e_ice(0.98) {} // [-] emissivity of ice, EngineeringToolbox.com }; diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/seb_threecomponent_evaluator.cc b/src/pks/surface_balance/constitutive_relations/land_cover/seb_threecomponent_evaluator.cc index 1cd3245ca..259501217 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/seb_threecomponent_evaluator.cc +++ b/src/pks/surface_balance/constitutive_relations/land_cover/seb_threecomponent_evaluator.cc @@ -614,8 +614,11 @@ SEBThreeComponentEvaluator::EnsureCompatibility_ToDeps_(State& S) if (land_cover_.size() == 0) land_cover_ = getLandCover(S.ICList().sublist("land cover types"), - {"roughness_snow", "roughness_ground", - "water_transition_depth", "dessicated_zone_thickness"}); + {"roughness_snow", + "roughness_ground", + "water_transition_depth", + "dessicated_zone_thickness", + "clapp_horn_b"}); // use domain name to set the mesh type CompositeVectorSpace domain_fac; diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/seb_twocomponent_evaluator.cc b/src/pks/surface_balance/constitutive_relations/land_cover/seb_twocomponent_evaluator.cc index 339656108..d161d792b 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/seb_twocomponent_evaluator.cc +++ b/src/pks/surface_balance/constitutive_relations/land_cover/seb_twocomponent_evaluator.cc @@ -281,6 +281,7 @@ SEBTwoComponentEvaluator::Evaluate_(const State& S, if (model_1p1_) surf.density_w = 1000.; else surf.density_w = mass_dens[0][c]; surf.dz = lc.second.dessicated_zone_thickness; + surf.clapp_horn_b = lc.second.clapp_horn_b; surf.albedo = sg_albedo[0][c]; surf.emissivity = emissivity[0][c]; @@ -350,6 +351,7 @@ SEBTwoComponentEvaluator::Evaluate_(const State& S, if (model_1p1_) surf.density_w = 1000; else surf.density_w = mass_dens[0][c]; surf.dz = lc.second.dessicated_zone_thickness; + surf.clapp_horn_b = lc.second.clapp_horn_b; surf.albedo = sg_albedo[1][c]; surf.emissivity = emissivity[1][c]; surf.water_transition_depth = lc.second.water_transition_depth; @@ -536,9 +538,12 @@ SEBTwoComponentEvaluator::EnsureCompatibility_ToDeps_(State& S) if (land_cover_.size() == 0) land_cover_ = getLandCover(S.ICList().sublist("land cover types"), - {"roughness_snow", "roughness_ground", - "water_transition_depth", "snow_transition_depth", - "dessicated_zone_thickness"}); + {"roughness_snow", + "roughness_ground", + "water_transition_depth", + "snow_transition_depth", + "dessicated_zone_thickness", + "clapp_horn_b"}); CompositeVectorSpace domain_fac; domain_fac.SetMesh(S.GetMesh(domain_)) diff --git a/testing/ats-regression-tests b/testing/ats-regression-tests index c1c484863..a4593ba66 160000 --- a/testing/ats-regression-tests +++ b/testing/ats-regression-tests @@ -1 +1 @@ -Subproject commit c1c484863b9e651edb42352b3ebab2af5edb0f5b +Subproject commit a4593ba6656e5b6de81583875279bae26b61a16c diff --git a/tools/input_converters/xml-main-new_state.py b/tools/input_converters/xml-main-new_state.py index d91c6ebf9..c882c276e 100644 --- a/tools/input_converters/xml-main-new_state.py +++ b/tools/input_converters/xml-main-new_state.py @@ -124,13 +124,24 @@ def top_cell_evals(xml): if ev_type.getValue() == 'surface top cell evaluator': ev_type.setValue('surface from top cell evaluator') +def density_units(xml): + """[kg/m^3] --> [kg m^-3] to be consistent with all other units in ATS""" + for d in asearch.findall_path(xml, ["molar mass [kg/mol]"]): + d.setName("molar mass [kg mol^-1]") + for d in asearch.findall_path(xml, ["molar mass [g/mol]"]): + d.setName("molar mass [g mol^-1]") + for d in asearch.findall_path(xml, ["density [kg/m^3]"]): + d.setName("density [kg m^-3]") + for d in asearch.findall_path(xml, ["density [mol/m^3]"]): + d.setName("density [mol m^-3]") + def update(xml, has_orig=False): """generic update calls all needed things""" new_state(xml) pk_initial_timestep(xml, has_orig) top_cell_evals(xml) - + density_units(xml) if __name__ == "__main__": import argparse diff --git a/tools/testing/test_manager.py b/tools/testing/test_manager.py index ce4e4e5ff..c69dc4d18 100644 --- a/tools/testing/test_manager.py +++ b/tools/testing/test_manager.py @@ -36,8 +36,10 @@ -aliases = {'surface-water_flux':['surface-mass_flux','surface-mass_flux_next'], - 'water_flux':['mass_flux','mass_flux_next']} +_aliases = {'surface-water_flux':['surface-mass_flux','surface-mass_flux_next'], + 'water_flux':['mass_flux','mass_flux_next'], + 'saturation_liquid':['prev_saturation_liquid'], + } def get_git_hash(directory): @@ -134,6 +136,7 @@ def __init__(self, executable='ats', mpiexec=None, version=None, suffix=None): self._check_performance = False self._num_failed = 0 self._test_name = None + self._new_checkpoint_files = False # docker goodies self._docker_executable = "docker" @@ -241,10 +244,11 @@ def filenames(self, dirname): fname_list2_tmp = sorted(glob.glob(os.path.join(dirname, self._file_prefix+"*", "*"+self._file_suffix))) fname_list2 = [f for f in fname_list2_tmp if 'checkpoint_final' not in f] - all_fnames = fname_list + fname_list2 - #print(f'Files in {dirname}: {len(all_fnames)}') - #print(all_fnames) - return all_fnames + assert(len(fname_list) == 0 or len(fname_list2) == 0) + if len(fname_list2) > 0: + self._new_checkpoint_files = True + + return fname_list + fname_list2 def run(self, dry_run, status, testlog): """ @@ -648,54 +652,116 @@ def _check_gold(self, status, testlog): def _compare(self, h5_current, h5_gold, status, testlog): """Check that output hdf5 file has not changed from the baseline. """ - for key, tolerance in self._criteria.items(): - if key == self._TIME and 'time' in h5_gold.attrs: - self._check_tolerance(h5_current.attrs['time'], h5_gold.attrs['time'], - self._TIME, tolerance, status, testlog) + def split_name(crit_name): + """Parse crit_name for name,tag,entity,dof""" + crit_split = crit_name.split('.') + name_tag = crit_split[0] + if len(crit_split) > 1: + entity = crit_split[1] + else: + entity = None + if len(crit_split) > 2: + dof = crit_split[2] else: - # find all matches of this key - reg_matches = [] - if len(key.split('.')) == 1: - for k in h5_current.keys(): - if k.split('.')[0] == key: - reg_matches.append(k) - elif len(key.split('.')) == 2: - for k in h5_current.keys(): - if '.'.join(k.split('.')[0:2]) == key: - reg_matches.append(k) - else: - for k in h5_current.keys(): - if k == key: - reg_matches.append(k) + dof = None + + if '@' in name_tag: + name, tag = name_tag.split('@') + else: + name = name_tag + tag = None + return name, tag, entity, dof + + def join_name(name, tag=None, entity=None, dof=None): + """Join split name back into single string""" + res = name + if tag is not None: + res = '@'.join((res, tag)) + if entity is not None: + res = '.'.join((res, entity)) + if dof is not None: + res = '.'.join((res, dof)) + return res + + def find_match(target, container): + """Custom search for split_name target in an iterable list of potential matches""" + matches = [] + for k in container: + ks = split_name(k) + if target[0] == ks[0] and \ + (target[1] == None or target[1] == ks[1]) and \ + (target[2] == None or target[2] == ks[2]) and \ + (target[3] == None or target[3] == ks[3]): + matches.append(k) + + # filter out matches -- if a target without a tag is + # matched by a key without a tag, just check those and do + # NOT check keys that DO have a tag. + if crit[1] is None and any('@' not in k for k in matches): + matches = [k for k in matches if '@' not in k] + + # if that failed to find any, try to ignore the tag + if len(matches) == 0 and target[1] is not None: + matches = find_match((target[0], None, target[2], target[3]), container) + + # if that failed to find any, look for aliased versions + if len(matches) == 0 and target[0] in _aliases: + for alias in _aliases[target[0]]: + matches = find_match((alias, target[1], target[2], target[3]), container) + if len(matches) > 0: + break + + # lastly, look for old/dead suffixes + if len(matches) == 0 and target[3] is not None and not target[3].endswith(" conc"): + matches = find_match((target[0], target[1], target[2], target[3]+" conc"), container) + + return matches + + + for crit_string, tolerance in self._criteria.items(): + # if this is new-style checkpoint files, we have to check domain names + if self._new_checkpoint_files: + try: + a_gold_key = list(h5_gold.keys())[0] + except IndexError: + continue + + gold_domain_name = split_name(a_gold_key)[0].split('-')[0] + crit_domain_name = split_name(crit_string)[0].split('-')[0] + if gold_domain_name != crit_domain_name: + continue + + if crit_string == self._TIME: + if 'time' in h5_current.attrs and 'time' in h5_gold.attrs: + self._check_tolerance(h5_current.attrs['time'], h5_gold.attrs['time'], + self._TIME, tolerance, status, testlog) + + else: + # find all matches of this crit_string + crit = split_name(crit_string) + reg_matches = find_match(crit, h5_current.keys()) # find the corresponding matches from gold gold_matches = [] - for key in reg_matches: - if key in h5_gold.keys(): - gold_matches.append(key) + for match in reg_matches: + if match in h5_gold.keys(): + gold_matches.append(match) else: - found_match = False - key_split = key.split('.') - try: - my_aliases = aliases[key_split[0]] - except KeyError: - pass - else: - for alias in my_aliases: - potential_alias = '.'.join([alias,]+key_split[1:]) - if potential_alias in h5_gold.keys(): - gold_matches.append(potential_alias) - found_match = True - break - if not found_match: + this_matches = find_match(split_name(match), h5_gold.keys()) + if len(this_matches) == 0 or len(this_matches) > 1: status.fail = 1 - print(" FAIL: Cannot find {0} or aliased version in the regression.".format(key), + print(" FAIL: Cannot find {0} or aliased version in the regression.".format(match), file=testlog) return + gold_matches.append(this_matches[0]) - # we should not get here if this is not true assert(len(gold_matches) == len(reg_matches)) + if len(gold_matches) == 0: + status.error = 1 + print("ERROR: Cannot find {0} or aliased version in the gold file -- " + "bad criteria.".format(crit_string), file=testlog) + return for (g, r) in zip(gold_matches, reg_matches): self._check_tolerance(h5_current[r][:], h5_gold[g][:], r, tolerance, status, testlog) @@ -1343,20 +1409,12 @@ def search_for_config_files(base_dir, config_file_list): if filename.endswith(".cfg") and os.path.isfile(os.path.join(root, filename)): config_file_list.append(os.path.join(root, filename)) + def check_options(options): """ Run some sanity checks on the commandline options. """ - # prevent the user from updating regression output during a - # recursive search for config files - if hasattr(options, 'update'): - if options.update and config_list_includes_search(options): - raise RuntimeError("ERROR: cannot update gold regression files " - "during a recursive search for config files.") - - if options.update and options.new_tests: - raise RuntimeError("ERROR: cannot create new tests and update gold " - "regression files at the same time.") + pass def check_for_docker_image(docker_image): diff --git a/tools/utils/plot_timestep_history.py b/tools/utils/plot_timestep_history.py index 8a06f460d..2dafbc583 100755 --- a/tools/utils/plot_timestep_history.py +++ b/tools/utils/plot_timestep_history.py @@ -34,6 +34,41 @@ def parse_logfile(fid, wallclock=False): print(f" bad timesteps = {nbad}") return np.array(data), np.array(faildata) +def parse_logfile_subcycle_pk(fid, pk_name, wallclock=False): + """Reads a file, and returns a list of good and bad timesteps. + + Each are a list of 3-tuples: (step number, time of step, dt) + """ + + data = [] + faildata = [] + + # find number of "cycles" + total_cycles = 0 + cycle = 0 + t0_prev = None + for line in fid: + if line.startswith(pk_name) and "Advancing:" in line: + postline = line.split('|')[1] + sline = postline.split() + t0 = float(sline[3]) + dt = float(sline[9]) + if t0_prev is None: + t0_prev = t0 + elif t0 > t0_prev: + t0_prev = t0 + cycle += 1 + else: + faildata.append(data.pop()) + data.append([cycle,t0,dt]) + + ngood = len(data) + nbad = len(faildata) + print(f" total timesteps = {ngood+nbad}") + print(f" bad timesteps = {nbad}") + return np.array(data), np.array(faildata) + + def get_axs(): """Gets a figure and list of axes for plotting.""" @@ -66,9 +101,12 @@ def plot(data, axs, color, label, symbol='x', units='s'): if data[1].shape != (0,): axs[1].semilogy(data[1][:,0], data[1][:,2], symbol, color=color) -def write_to_file(data, fnamebase): +def write_to_file(data, fnamebase, subcycle=None): """Writes the data to a file for future reading""" - fname = fnamebase+".npz" + if subcycle is not None: + fname = fnamebase+subcycle+".npz" + else: + fname = fnamebase+".npz" np.savez(fname, good_timesteps=data[0], bad_timesteps=data[1]) def read_from_file(fname): @@ -98,6 +136,8 @@ def read_from_file(fname): help="Do not use any existing .npz file -- instead reload from the logfile and overwrite the .npz file.") parser.add_argument("--units", type=str, default='d', choices=['y', 'd', 's'], help="Set the units of the x axis -- one of 'y', 'd', 's'") + parser.add_argument("--subcycled", type=str, action="append", + help="View data from a subcycled PK by this name.") args = parser.parse_args() if args.colors is not None: @@ -107,17 +147,29 @@ def read_from_file(fname): args.colors = colors.enumerated_colors(len(args.LOG_FILES)) fig, axs = get_axs() + + if args.subcycled is None or len(args.subcycled) == 0: + args.subcycled = [None,] * len(args.LOG_FILES) + elif len(args.subcycled) == 1: + args.subcycled = args.subcycled * len(args.LOG_FILES) + elif len(args.subcycled) != len(args.LOG_FILES): + raise ValueError('If subcycled option is provided, it must be provided for all files, or exactly once to be used for all files.') - for fname, color in zip(args.LOG_FILES, args.colors): + for fname, color, subcycle in zip(args.LOG_FILES, args.colors, args.subcycled): if fname.endswith(".npz"): data = read_from_file(fname) - elif os.path.isfile(fname+".npz") and not args.overwrite: + elif os.path.isfile(fname+".npz") and not args.overwrite and args.subcycled is None: data = read_from_file(fname+".npz") else: print(f"File: {fname}") + with open(fname,'r') as fid: - data = parse_logfile(fid) - write_to_file(data, fname) + if subcycle is None: + data = parse_logfile(fid) + else: + data = parse_logfile_subcycle_pk(fid, subcycle) + + write_to_file(data, fname, subcycle) plot(data, axs, color, fname, units=args.units) diff --git a/tools/visit_ats/visit_ats/visit_ats.py b/tools/visit_ats/visit_ats/visit_ats.py index 1cc6e4638..813ab72f0 100644 --- a/tools/visit_ats/visit_ats/visit_ats.py +++ b/tools/visit_ats/visit_ats/visit_ats.py @@ -117,7 +117,7 @@ def _exaggerateVertical(self): if "exaggerate_vertical" == op.oname: done = True if not done: - print "transforming plot %d..."%i + print("transforming plot %d..."%i) tr = v.TransformAttributes() tr.doScale = 1 tr.scaleY = self.exaggeration @@ -163,10 +163,10 @@ def createPseudocolor(self, varname, display_name=None, cmap=None, if "temperature" in display_name: display_name = display_name.replace("[K]", "[C]") - print "defining alias: %s = %s"%(display_name, varname) + print("defining alias: %s = %s"%(display_name, varname)) v.DefineScalarExpression(display_name, "<%s> - 273.15"%varname) elif display_name != varname: - print "defining alias: %s = %s"%(display_name, varname) + print("defining alias: %s = %s"%(display_name, varname)) v.DefineScalarExpression(display_name, '<'+varname+'>') v.AddPlot('Pseudocolor', display_name) @@ -248,7 +248,7 @@ def createContour(self, varname, value, color=None, linewidth=None): return plot def draw(self): - print "drawing window %d of dimension %d"%(self.i,self.dim) + print("drawing window %d of dimension %d"%(self.i,self.dim)) v.SetActiveWindow(self.i) v.SetAnnotationAttributes(self.annot) @@ -262,7 +262,7 @@ def draw(self): sliced = True if not sliced: - print "slicing plot %d..."%i + print("slicing plot %d..."%i) v.SetActivePlots(i) v.AddOperator("Slice") sa = self._slice.toAttributes() @@ -271,23 +271,23 @@ def draw(self): if self.exaggeration is not None: - print "exaggerating..." + print("exaggerating...") self._exaggerateVertical() # set the plot options for i, plot in enumerate(self.plots): - print "setting plot options for plot %i..."%i + print("setting plot options for plot %i..."%i) v.SetActivePlots(i) v.SetPlotOptions(plot.patts) # set the view - print "setting the view..." + print("setting the view...") if self.dim == 2: v.SetView2D(self.view) else: v.SetView3D(self.view) - print "drawing..." + print("drawing...") v.DrawPlots() diff --git a/tools/visit_ats/visit_ats/visit_colormaps.py b/tools/visit_ats/visit_ats/visit_colormaps.py index 2867ae601..c4715ad99 100644 --- a/tools/visit_ats/visit_ats/visit_colormaps.py +++ b/tools/visit_ats/visit_ats/visit_colormaps.py @@ -13,7 +13,7 @@ def createTemperaturesColorMap(Tmin, Tmax, zero=0., name="temperature"): """Creates a ColorMap with warm colors above zero, and cold colors below.""" # color control point format: (position, (r,g,b,a)) for (r,g,b,a) in [0,255] if name in _added_temps: - print "Already added a temperature color map of name %s"%name + print("Already added a temperature color map of name %s"%name) return if Tmin >= Tmax: @@ -54,7 +54,7 @@ def _interp(mymin, mymax, colors): cp.position = val[0] cp.colors = val[1] clist.AddControlPoints(cp) - print "Colormap Temp", name, [val[0] for val in cmap] + print("Colormap Temp", name, [val[0] for val in cmap]) v.AddColorTable(name, clist) _added_temps.append(name) diff --git a/tools/visit_ats/visit_ats/visit_time.py b/tools/visit_ats/visit_ats/visit_time.py index 8a5b7ded5..d5db97631 100644 --- a/tools/visit_ats/visit_ats/visit_time.py +++ b/tools/visit_ats/visit_ats/visit_time.py @@ -48,7 +48,7 @@ def roundTime(time, round=None): tround = tprev else: tround = tnext - print "Rounding: %s to %s"%(time.strftime("%c"), tround.strftime("%c")) + print("Rounding: %s to %s"%(time.strftime("%c"), tround.strftime("%c"))) return tround