diff --git a/src/constitutive_relations/CMakeLists.txt b/src/constitutive_relations/CMakeLists.txt index 147c7e796..3d6d42220 100644 --- a/src/constitutive_relations/CMakeLists.txt +++ b/src/constitutive_relations/CMakeLists.txt @@ -44,6 +44,11 @@ register_evaluator_with_factory( LISTNAME ATS_RELATIONS_REG ) +register_evaluator_with_factory( + HEADERFILE surface_subsurface_fluxes/flux_divergence_evaluator_reg.hh + LISTNAME ATS_RELATIONS_REG + ) + register_evaluator_with_factory( HEADERFILE generic_evaluators/MultiplicativeEvaluator_reg.hh diff --git a/src/constitutive_relations/surface_subsurface_fluxes/CMakeLists.txt b/src/constitutive_relations/surface_subsurface_fluxes/CMakeLists.txt index 168b568b8..d4c7a1805 100644 --- a/src/constitutive_relations/surface_subsurface_fluxes/CMakeLists.txt +++ b/src/constitutive_relations/surface_subsurface_fluxes/CMakeLists.txt @@ -10,6 +10,7 @@ set(ats_surf_subsurf_src_files surface_top_cells_evaluator.cc top_cells_surface_evaluator.cc volumetric_darcy_flux_evaluator.cc + flux_divergence_evaluator.cc ) file(GLOB ats_surf_subsurf_inc_files "*.hh") diff --git a/src/constitutive_relations/surface_subsurface_fluxes/flux_divergence_evaluator.cc b/src/constitutive_relations/surface_subsurface_fluxes/flux_divergence_evaluator.cc new file mode 100644 index 000000000..16969577e --- /dev/null +++ b/src/constitutive_relations/surface_subsurface_fluxes/flux_divergence_evaluator.cc @@ -0,0 +1,65 @@ +/* + Copyright 2010-202x held jointly by participating institutions. + ATS is released under the three-clause BSD License. + The terms of use and "as is" disclaimer for this license are + provided in the top-level COPYRIGHT file. + + Authors: Ethan Coon (coonet@ornl.gov) +*/ + +#include "flux_divergence_evaluator.hh" + +namespace Amanzi { +namespace Relations { + +FluxDivergenceEvaluator::FluxDivergenceEvaluator(Teuchos::ParameterList& plist) + : EvaluatorSecondaryMonotypeCV(plist) +{ + // determine the domain + Key akey = my_keys_.front().first; + auto tag = my_keys_.front().second; + Key domain = Keys::getDomain(akey); + + flux_key_ = Keys::readKey(plist, domain, "flux", "water_flux"); + dependencies_.insert(KeyTag{ flux_key_, tag }); +} + +Teuchos::RCP +FluxDivergenceEvaluator::Clone() const +{ + return Teuchos::rcp(new FluxDivergenceEvaluator(*this)); +} + +void +FluxDivergenceEvaluator::EnsureCompatibility_ToDeps_(State& S) +{ + const auto& my_fac = S.Require(my_keys_.front().first, my_keys_.front().second); + if (my_fac.Mesh() != Teuchos::null) { + for (auto dep : dependencies_) { + auto& fac = S.Require(dep.first, dep.second); + fac.SetMesh(my_fac.Mesh())->AddComponent("face", AmanziMesh::Entity_kind::FACE, 1); + } + } +} + + +void +FluxDivergenceEvaluator::Evaluate_(const State& S, const std::vector& result) +{ + auto tag = my_keys_.front().second; + const auto& flux_f = *S.Get(flux_key_, tag).ViewComponent("face", true); + const AmanziMesh::Mesh& m = *result[0]->Mesh(); + auto& div_flux = *result[0]->ViewComponent("cell", false); + div_flux.PutScalar(0.); + + for (AmanziMesh::Entity_ID c = 0; c != div_flux.MyLength(); ++c) { + auto [cfaces, cfdirs] = m.getCellFacesAndDirections(c); + for (int i = 0; i != cfaces.size(); ++i) { + div_flux[0][c] += flux_f[0][cfaces[i]] * cfdirs[i]; + } + } +} + + +} // namespace Relations +} // namespace Amanzi diff --git a/src/constitutive_relations/surface_subsurface_fluxes/flux_divergence_evaluator.hh b/src/constitutive_relations/surface_subsurface_fluxes/flux_divergence_evaluator.hh new file mode 100644 index 000000000..e730d9ed1 --- /dev/null +++ b/src/constitutive_relations/surface_subsurface_fluxes/flux_divergence_evaluator.hh @@ -0,0 +1,59 @@ +/* + Copyright 2010-202x held jointly by participating institutions. + ATS is released under the three-clause BSD License. + The terms of use and "as is" disclaimer for this license are + provided in the top-level COPYRIGHT file. + + Authors: Ethan Coon (coonet@ornl.gov) +*/ + +/* + An evaluator for compute the divergence of a flux field. +*/ +#pragma once + +#include "dbc.hh" +#include "Evaluator_Factory.hh" +#include "EvaluatorSecondaryMonotype.hh" + +namespace Amanzi { +namespace Relations { + +class FluxDivergenceEvaluator : public EvaluatorSecondaryMonotypeCV { + public: + explicit FluxDivergenceEvaluator(Teuchos::ParameterList& plist); + FluxDivergenceEvaluator(const FluxDivergenceEvaluator& other) = default; + Teuchos::RCP Clone() const override; + + // turn off all derivatives manually + virtual bool IsDifferentiableWRT(const State& S, + const Key& wrt_key, + const Tag& wrt_tag) const override { + return false; + } + +protected: + // custom ensure compatibility as all data is not just on the same components + virtual 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 { + AMANZI_ASSERT(false); + } + + protected: + Key flux_key_; + + private: + static Utils::RegisteredFactory reg_; +}; + + +} // namespace Relations +} // namespace Amanzi + + diff --git a/src/constitutive_relations/surface_subsurface_fluxes/flux_divergence_evaluator_reg.hh b/src/constitutive_relations/surface_subsurface_fluxes/flux_divergence_evaluator_reg.hh new file mode 100644 index 000000000..70bbbfc47 --- /dev/null +++ b/src/constitutive_relations/surface_subsurface_fluxes/flux_divergence_evaluator_reg.hh @@ -0,0 +1,20 @@ +/* + Copyright 2010-202x held jointly by participating institutions. + ATS is released under the three-clause BSD License. + The terms of use and "as is" disclaimer for this license are + provided in the top-level COPYRIGHT file. + + Authors: Ethan Coon (coonet@ornl.gov) +*/ + +#include "flux_divergence_evaluator.hh" + +namespace Amanzi { +namespace Relations { + +// registry of method +Utils::RegisteredFactory + FluxDivergenceEvaluator::reg_("flux divergence"); + +} // namespace Relations +} // namespace Amanzi diff --git a/src/executables/ats_driver.cc b/src/executables/ats_driver.cc index 498569c18..733b69bc8 100644 --- a/src/executables/ats_driver.cc +++ b/src/executables/ats_driver.cc @@ -69,10 +69,16 @@ ATSDriver::cycle_driver() << "Beginning setup stage..." << std::endl << std::flush; } + parseParameterList(); + if (vo_->os_OK(Teuchos::VERB_LOW)) { + *vo_->os() << " ... completed: "; + reportOneTimer_("2a: parseParameterList"); + } + setup(); if (vo_->os_OK(Teuchos::VERB_LOW)) { *vo_->os() << " ... completed: "; - reportOneTimer_("2: setup"); + reportOneTimer_("2b: setup"); } // @@ -182,7 +188,6 @@ ATSDriver::cycle_driver() reportOneTimer_("4: solve"); } - // finalizing simulation if (vo_->os_OK(Teuchos::VERB_LOW)) { *vo_->os() << "================================================================================" diff --git a/src/executables/coordinator.cc b/src/executables/coordinator.cc index b2023c8e6..60310319c 100644 --- a/src/executables/coordinator.cc +++ b/src/executables/coordinator.cc @@ -78,7 +78,8 @@ Coordinator::Coordinator(const Teuchos::RCP& plist, timers_["0: create mesh"] = Teuchos::TimeMonitor::getNewCounter("0: create mesh"); timers_["1: create run"] = Teuchos::TimeMonitor::getNewCounter("1: create run"); - timers_["2: setup"] = Teuchos::TimeMonitor::getNewCounter("2: setup"); + timers_["2a: parseParameterList"] = Teuchos::TimeMonitor::getNewCounter("2a: parseParameterList"); + timers_["2b: setup"] = Teuchos::TimeMonitor::getNewCounter("2b: setup"); timers_["3: initialize"] = Teuchos::TimeMonitor::getNewCounter("3: initialize"); timers_["4: solve"] = Teuchos::TimeMonitor::getNewCounter("4: solve"); timers_["4a: advance step"] = Teuchos::TimeMonitor::getNewCounter("4a: advance step"); @@ -272,24 +273,31 @@ Coordinator::Coordinator(const Teuchos::RCP& plist, } +void +Coordinator::parseParameterList() +{ + Teuchos::TimeMonitor monitor(*timers_.at("2a: parseParameterList")); + + // require needed times -- these may be targets as dependencies + S_->require_time(Amanzi::Tags::CURRENT); + S_->require_time(Amanzi::Tags::NEXT); + + pk_->set_tags(Amanzi::Tags::CURRENT, Amanzi::Tags::NEXT); + pk_->parseParameterList(); +} + void Coordinator::setup() { - Teuchos::TimeMonitor monitor(*timers_.at("2: setup")); + Teuchos::TimeMonitor monitor(*timers_.at("2b: setup")); // common constants S_->Require("atmospheric_pressure", Amanzi::Tags::DEFAULT, "coordinator"); S_->Require("gravity", Amanzi::Tags::DEFAULT, "coordinator"); - // needed other times - S_->require_time(Amanzi::Tags::CURRENT); - S_->require_time(Amanzi::Tags::NEXT); - // order matters here -- PK::Setup() set the leaves, then observations can // use those if provided, and State::Setup finally deals with all secondaries // and allocates memory - pk_->set_tags(Amanzi::Tags::CURRENT, Amanzi::Tags::NEXT); - pk_->parseParameterList(); pk_->Setup(); for (auto& obs : observations_) obs->Setup(S_.ptr()); S_->Setup(); diff --git a/src/executables/coordinator.hh b/src/executables/coordinator.hh index d7d2c5350..b7f8741fa 100644 --- a/src/executables/coordinator.hh +++ b/src/executables/coordinator.hh @@ -46,6 +46,7 @@ class Coordinator { const Amanzi::Comm_ptr_type& comm); // PK methods + void parseParameterList(); void setup(); void initialize(); void finalize(); diff --git a/src/executables/elm_ats_api/ats_variables.F90 b/src/executables/elm_ats_api/ats_variables.F90 new file mode 100644 index 000000000..992c8df8b --- /dev/null +++ b/src/executables/elm_ats_api/ats_variables.F90 @@ -0,0 +1,34 @@ +module ats_variables + use, intrinsic :: iso_c_binding + implicit none + + type :: ats_var_id_type + integer(c_int), parameter :: NULL = 0 + integer(c_int), parameter :: BASE_POROSITY = 1 + integer(c_int), parameter :: HYDRAULIC_CONDUCTIVITY = 2 + integer(c_int), parameter :: CLAPP_HORNBERGER_B = 3 + integer(c_int), parameter :: CLAPP_HORNBERGER_PSI_SAT = 4 + integer(c_int), parameter :: RESIDUAL_SATURATION = 5 + integer(c_int), parameter :: EFFECTIVE_POROSITY = 6 + integer(c_int), parameter :: ROOT_FRACTION = 7 + integer(c_int), parameter :: SURFACE_WATER_CONTENT = 8 + integer(c_int), parameter :: WATER_CONTENT = 9 + integer(c_int), parameter :: PRESSURE = 10 + integer(c_int), parameter :: GROSS_SURFACE_WATER_SOURCE = 11 + integer(c_int), parameter :: POTENTIAL_EVAPORATION = 12 + integer(c_int), parameter :: POTENTIAL_TRANSPIRATION = 13 + integer(c_int), parameter :: EVAPORATION = 14 + integer(c_int), parameter :: COLUMN_TRANSPIRATION = 15 + integer(c_int), parameter :: RUNOFF = 16 + integer(c_int), parameter :: BASEFLOW = 17 + integer(c_int), parameter :: ELEVATION = 18 + integer(c_int), parameter :: TIME = 19 + integer(c_int), parameter :: INITIAL_WATER_CONTENT = 20 + integer(c_int), parameter :: SATURATION_LIQUID = 21 + integer(c_int), parameter :: PONDED_DEPTH = 22 + end type ats_var_id_type + + type(ats_var_id_type), parameter :: ats_var_id = ats_var_id_type() + +end module ats_variables + diff --git a/src/executables/elm_ats_api/ats_variables.hh b/src/executables/elm_ats_api/ats_variables.hh new file mode 100644 index 000000000..f2c13a466 --- /dev/null +++ b/src/executables/elm_ats_api/ats_variables.hh @@ -0,0 +1,34 @@ +namespace ATS { +namespace ELM { + +enum class VarID : int { + NO_VARIABLE = 0, + BASE_POROSITY = 1, + HYDRAULIC_CONDUCTIVITY = 2, + CLAPP_HORNBERGER_B = 3, + CLAPP_HORNBERGER_PSI_SAT = 4, + RESIDUAL_SATURATION = 5, + EFFECTIVE_POROSITY = 6, + ROOT_FRACTION = 7, + SURFACE_WATER_CONTENT = 8, + WATER_CONTENT = 9, + PRESSURE = 10, + GROSS_SURFACE_WATER_SOURCE = 11, + POTENTIAL_EVAPORATION = 12, + POTENTIAL_TRANSPIRATION = 13, + EVAPORATION = 14, + TRANSPIRATION = 15, + RUNOFF = 16, + BASEFLOW = 17, + ELEVATION = 18, + TIME = 19, + INITIAL_WATER_CONTENT = 20, + SATURATION_LIQUID = 21, + PONDED_DEPTH = 22, + DEPTH_TO_WATER_TABLE = 23, + WATER_CONTENT_OLD = 24, + SURFACE_WATER_CONTENT_OLD = 25 +}; + +} // namespace ELM +} // namespace ATS diff --git a/src/executables/elm_ats_api/elm_ats_api.cc b/src/executables/elm_ats_api/elm_ats_api.cc index b113e24f1..9d2c38139 100644 --- a/src/executables/elm_ats_api/elm_ats_api.cc +++ b/src/executables/elm_ats_api/elm_ats_api.cc @@ -21,149 +21,127 @@ namespace { extern "C" { + #endif - // allocate, call constructor and cast ptr to opaque ELM_ATSDriver_ptr - ELM_ATSDriver_ptr ats_create(MPI_Fint* f_comm, const char* input_filename) - { - // Initialize Kokkos if ELM hasn't already (for near future, it won't) - if (!Kokkos::is_initialized()) { - Kokkos::initialize(); - ats_kokkos_init = true; - } - return reinterpret_cast(ATS::createELM_ATSDriver(f_comm, input_filename)); +// allocate, call constructor and cast ptr to opaque ELM_ATSDriver_ptr +ELM_ATSDriver_ptr ats_create(MPI_Fint* f_comm, + const char* input_filename, + const char* logfile_filename) +{ + // Initialize Kokkos if ELM hasn't already (for near future, it won't) + if (!Kokkos::is_initialized()) { + Kokkos::initialize(); + ats_kokkos_init = true; } + return reinterpret_cast(ATS::createELM_ATSDriver(f_comm, input_filename, logfile_filename)); +} - // reinterpret as elm_ats_driver and delete (calls destructor) - void ats_delete(ELM_ATSDriver_ptr ats) - { - auto ats_ptr = reinterpret_cast(ats); - ats_ptr->finalize(); - delete ats_ptr; - - // If ATS initialized Kokkos, then finalize it here - if (ats_kokkos_init && Kokkos::is_initialized()) { - Kokkos::finalize(); - ats_kokkos_init = false; - } - } - // call driver advance() - void ats_advance(ELM_ATSDriver_ptr ats, - double const* const dt, - bool const* const checkpoint, - bool const* const visualize) - { - reinterpret_cast(ats)->advance(*dt, *checkpoint, *visualize); +// reinterpret as elm_ats_driver and delete (calls destructor) +void ats_delete(ELM_ATSDriver_ptr ats) +{ + auto ats_ptr = reinterpret_cast(ats); + ats_ptr->finalize(); + delete ats_ptr; + + // If ATS initialized Kokkos, then finalize it here + if (ats_kokkos_init && Kokkos::is_initialized()) { + Kokkos::finalize(); + ats_kokkos_init = false; } +} - // call driver advance_test() - void ats_advance_test(ELM_ATSDriver_ptr ats) - { - reinterpret_cast(ats)->advance_test(); - } - void ats_get_mesh_info(ELM_ATSDriver_ptr ats, - int* const ncols_local, - int* const ncols_global, - double* const lat, - double* const lon, - double* const elev, - double* const surf_area, - int* const pft, - int* const nlevgrnd, - double* const depth) - { - reinterpret_cast(ats)->get_mesh_info( - *ncols_local, *ncols_global, lat, lon, elev, surf_area, pft, *nlevgrnd, depth); - } +// call driver setup() +void ats_parse_parameter_list(ELM_ATSDriver_ptr ats) +{ + reinterpret_cast(ats)->parseParameterList(); +} - // call driver setup() - void ats_setup(ELM_ATSDriver_ptr ats) - { - reinterpret_cast(ats)->setup(); - } +void ats_get_mesh_info(ELM_ATSDriver_ptr ats, + int * const ncols_local) +{ + *ncols_local = reinterpret_cast(ats)->ncolumns; +} - // call driver initialize() - void ats_initialize(ELM_ATSDriver_ptr ats, - double const* const t, - double const* const patm, - double const* const soilp) - { - reinterpret_cast(ats)->initialize(*t, patm, soilp); +void ats_get_mesh_info2(ELM_ATSDriver_ptr ats, + int * const ncols_local, + int * const ncols_global, + int * const nlevgrnd, + double * const dzs, + double * const areas, + double * const lat, + double * const lon) +{ + ATS::ELM_ATSDriver::MeshInfo info = reinterpret_cast(ats)->getMeshInfo(); + *ncols_local = info.ncols_local; + *ncols_global = info.ncols_global; + *nlevgrnd = info.nlevgrnd; + for (int i = 0; i != info.nlevgrnd; ++i) { + dzs[i] = info.dzs[i]; } - - - void ats_set_soil_hydrologic_parameters(ELM_ATSDriver_ptr ats, - double const* const base_porosity, - double const* const hydraulic_conductivity, - double const* const clapp_horn_b, - double const* const clapp_horn_smpsat, - double const* const clapp_horn_sr) - { - reinterpret_cast(ats)->set_soil_hydrologic_parameters( - base_porosity, hydraulic_conductivity, clapp_horn_b, clapp_horn_smpsat, clapp_horn_sr); + for (int i = 0; i != info.ncols_local; ++i) { + areas[i] = info.areas[i]; + lat[i] = info.latitudes[i]; + lon[i] = info.longitudes[i]; } +} - void ats_set_veg_parameters(ELM_ATSDriver_ptr ats, - double const* const mafic_potential_full_turgor, - double const* const mafic_potential_wilt_point) - { - reinterpret_cast(ats)->set_veg_parameters(mafic_potential_full_turgor, - mafic_potential_wilt_point); - } +// call driver setup() +void ats_setup(ELM_ATSDriver_ptr ats) +{ + reinterpret_cast(ats)->setup(); +} +// call driver initialize() +void ats_initialize(ELM_ATSDriver_ptr ats) +{ + reinterpret_cast(ats)->initialize(); +} - void ats_set_soil_hydrologic_properties(ELM_ATSDriver_ptr ats, - double const* const effective_porosity) - { - reinterpret_cast(ats)->set_soil_hydrologic_properties(effective_porosity); - } +// call driver advance(dt) +void ats_advance(ELM_ATSDriver_ptr ats, + double dt, + bool checkpoint, + bool visualize) +{ + reinterpret_cast(ats)->advance(dt, checkpoint, visualize); +} - // set veg properties, non-constant in time - void ats_set_veg_properties(ELM_ATSDriver_ptr ats, double const* const rooting_fraction) - { - reinterpret_cast(ats)->set_veg_properties(rooting_fraction); - } +// call driver advance_test() +void ats_advance_test(ELM_ATSDriver_ptr ats) +{ + reinterpret_cast(ats)->advanceTest(); +} +// +// Memory movement/data passing +// ----------------------------------------------------------------------------- +double ats_get_scalar(ELM_ATSDriver_ptr ats, int var_id) +{ + return reinterpret_cast(ats)->getScalar(static_cast(var_id)); +} - // call driver set_sources() - void ats_set_sources(ELM_ATSDriver_ptr ats, - double const* const surface_infiltration, - double const* const surface_evaporation, - double const* const subsurface_transpiration) - { - reinterpret_cast(ats)->set_potential_sources( - surface_infiltration, surface_evaporation, subsurface_transpiration); - } +void ats_set_scalar(ELM_ATSDriver_ptr ats, int var_id, double in) +{ + reinterpret_cast(ats)->setScalar(static_cast(var_id), in); +} +double const * ats_get_field_ptr(ELM_ATSDriver_ptr ats, int var_id) +{ + return reinterpret_cast(ats)->getFieldPtr(static_cast(var_id)); +} - void ats_get_waterstate(ELM_ATSDriver_ptr ats, - double* const surface_ponded_depth, - double* const water_table_depth, - double* const soil_pressure, - double* const soil_psi, - double* const sat_liq, - double* const sat_ice) - { - reinterpret_cast(ats)->get_waterstate( - surface_ponded_depth, water_table_depth, soil_pressure, soil_psi, sat_liq, sat_ice); - } +double * ats_get_field_ptr_w(ELM_ATSDriver_ptr ats, int var_id) +{ + return reinterpret_cast(ats)->getFieldPtrW(static_cast(var_id)); +} - void ats_get_water_fluxes(ELM_ATSDriver_ptr ats, - double* const soil_infiltration, - double* const evaporation, - double* const transpiration, - double* net_subsurface_fluxes, - double* net_runon) - { - reinterpret_cast(ats)->get_water_fluxes( - soil_infiltration, evaporation, transpiration, net_subsurface_fluxes, net_runon); - } #ifdef __cplusplus } diff --git a/src/executables/elm_ats_api/elm_ats_api.h b/src/executables/elm_ats_api/elm_ats_api.h index 0b2a26545..84b2ae359 100644 --- a/src/executables/elm_ats_api/elm_ats_api.h +++ b/src/executables/elm_ats_api/elm_ats_api.h @@ -19,7 +19,7 @@ for use with LSMs. extern "C" { // opaque pointer -// external caller only sees *ELM_ATSDriver_ptr - similar to void*, but better type safety +// external caller only sees *ELM_ATSDriver_ptr - similar to void*, but better type safety // ATS resolves ELM_ATSDriver_ptr as real ELM_ATSDriver during linking class ELM_ATSDriver; typedef ELM_ATSDriver *ELM_ATSDriver_ptr; @@ -32,98 +32,71 @@ typedef struct ELM_ATSDriver_ptr *ELM_ATSDriver_ptr; #endif // allocate, call constructor and cast ptr to opaque ELM_ATSDriver_ptr -ELM_ATSDriver_ptr ats_create(MPI_Fint *f_comm, const char *input_filename); +ELM_ATSDriver_ptr ats_create(MPI_Fint *f_comm, const char *input_filename, const char *logfile_filename); // reinterpret as elm_ats_driver and delete (calls destructor) void ats_delete(ELM_ATSDriver_ptr ats); -// call driver advance(dt) -void ats_advance(ELM_ATSDriver_ptr ats, - double const * const dt, - bool const * const checkpoint, - bool const * const visualize); - -// call driver advance_test() -void ats_advance_test(ELM_ATSDriver_ptr ats); // -// Called prior to run +// Control // ----------------------------------------------------------------------------- -// call driver get_mesh_info() +// +// simulation setup +// +void ats_parse_parameter_list(ELM_ATSDriver_ptr ats); + + +// +// This sets the decomposition +// void ats_get_mesh_info(ELM_ATSDriver_ptr ats, - int * const ncols_local, - int * const ncols_global, - double * const lat, - double * const lon, - double * const elev, - double * const surf_area, - int * const pft, - int * const nlevgrnd, - double * const depth); - -// call driver setup() -void ats_setup(ELM_ATSDriver_ptr ats); + int * const ncols_local); + +// +// These quantities should be compared against ELM to ensure consistent setup +// +void ats_get_mesh_info2(ELM_ATSDriver_ptr ats, + int * const ncols_local, + int * const ncols_global, + int * const nlevgrnd, + double * const dzs, + double * const areas, + double * const lat, + double * const lon); -// call driver initialize() -void ats_initialize(ELM_ATSDriver_ptr ats, - double const * const t, - double const * const patm, - double const * const soilp); - -// set material parameters, which are constant in time -void ats_set_soil_hydrologic_parameters(ELM_ATSDriver_ptr ats, - double const * const base_porosity, - double const * const hydraulic_conductivity, - double const * const clapp_horn_b, - double const * const clapp_horn_smpsat, - double const * const clapp_horn_sr); - -// set veg parameters, which are constant in time -void ats_set_veg_parameters(ELM_ATSDriver_ptr ats, - double const * const mafic_potential_full_turgor, - double const * const mafic_potential_wilt_point); +// +// simulation setup +// +void ats_setup(ELM_ATSDriver_ptr ats); // -// Called prior to timestep advance -// ----------------------------------------------------------------------------- -// set hydrologic properties, non-constant in time -void ats_set_soil_hydrologic_properties(ELM_ATSDriver_ptr ats, - double const * const effective_porosity); +// initial conditions +// +void ats_initialize(ELM_ATSDriver_ptr ats); -// set veg properties, non-constant in time -void ats_set_veg_properties(ELM_ATSDriver_ptr ats, - double const * const rooting_fraction); -// call driver set_sources() -// soil_infiltration & soil_evaporation are 1D arrays of length ncols -// root_transpiration is a 1D array array of length (ncells) -void ats_set_sources(ELM_ATSDriver_ptr ats, - double const * const surface_infiltration, - double const * const surface_evaporation, - double const * const subsurface_transpiration); +// call driver advance(dt) +void ats_advance(ELM_ATSDriver_ptr ats, + double dt, + bool checkpoint, + bool visualize); + +// call driver advance_test() +void ats_advance_test(ELM_ATSDriver_ptr ats); // -// Called after timestep advance +// Memory movement/data passing // ----------------------------------------------------------------------------- -// Get the water solution after the step is complete -// surface_pressure is a 1D array of length ncols -// soil_pressure & saturation are 1D arrays array of length (ncells) -void ats_get_waterstate(ELM_ATSDriver_ptr ats, - double * const surface_ponded_depth, - double * const water_table_depth, - double * const soil_pressure, - double * const soil_psi, - double * const sat_liq, - double * const sat_ice); - -void ats_get_water_fluxes(ELM_ATSDriver_ptr ats, - double * const soil_infiltration, - double * const evaporation, - double * const transpiration, - double * net_subsurface_fluxes, - double * net_runon); +double ats_get_scalar(ELM_ATSDriver_ptr ats, int var_id); +void ats_set_scalar(ELM_ATSDriver_ptr ats, int var_id, double in); + +void ats_get_field(ELM_ATSDriver_ptr ats, int var_id, double * const out); +double const * ats_get_field_ptr(ELM_ATSDriver_ptr ats, int var_id); +double * ats_get_field_ptr_w(ELM_ATSDriver_ptr ats, int var_id); +void ats_set_field(ELM_ATSDriver_ptr ats, int var_id, double const * const in); #ifdef __cplusplus } diff --git a/src/executables/elm_ats_api/elm_ats_driver.cc b/src/executables/elm_ats_api/elm_ats_driver.cc index 00d3b73de..55f89ca6a 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.cc +++ b/src/executables/elm_ats_api/elm_ats_driver.cc @@ -27,8 +27,7 @@ namespace ATS { ELM_ATSDriver* -createELM_ATSDriver(MPI_Fint* f_comm, const char* infile, int npfts) -{ +createELM_ATSDriver(MPI_Fint *f_comm, const char *infile, const char *logfile, int npfts) { // -- create communicator & get process rank //auto comm = getDefaultComm(); auto c_comm = MPI_Comm_f2c(*f_comm); @@ -39,10 +38,12 @@ createELM_ATSDriver(MPI_Fint* f_comm, const char* infile, int npfts) // convert input file to std::string for easier handling // infile must be null-terminated std::string input_filename(infile); + std::string logfile_filename(logfile); // check validity of input file name if (input_filename.empty()) { - if (rank == 0) std::cerr << "ERROR: no input file provided" << std::endl; + if (rank == 0) + std::cerr << "ERROR: no input file provided" << std::endl; } else if (!std::filesystem::exists(input_filename)) { if (rank == 0) std::cerr << "ERROR: input file \"" << input_filename << "\" does not exist." << std::endl; @@ -52,434 +53,355 @@ createELM_ATSDriver(MPI_Fint* f_comm, const char* infile, int npfts) Teuchos::RCP plist = Teuchos::getParametersFromXmlFile(input_filename); auto wallclock_timer = Teuchos::TimeMonitor::getNewCounter("wallclock duration"); - return new ELM_ATSDriver(plist, wallclock_timer, teuchos_comm, comm, npfts); + return new ELM_ATSDriver(plist, wallclock_timer, teuchos_comm, comm, logfile_filename, npfts); } + ELM_ATSDriver::ELM_ATSDriver(const Teuchos::RCP& plist, const Teuchos::RCP& wallclock_timer, const Teuchos::RCP>& teuchos_comm, const Amanzi::Comm_ptr_type& comm, + const std::string& logfile, int npfts) : Coordinator(plist, wallclock_timer, teuchos_comm, comm), npfts_(npfts), - ncolumns_(-1), - ncells_per_col_(-1) + ncolumns(-1), + ncells_per_col_(-1), + elm_plist_(Teuchos::sublist(Teuchos::sublist(plist, "cycle driver"), "ELM driver")) { // -- set default verbosity level to no output // -- TODO make the verbosity level an input argument - VerboseObject::global_default_level = Teuchos::VERB_NONE; + VerboseObject::global_default_level = Teuchos::VERB_HIGH; + VerboseObject::global_logfile = logfile; +} + +void +ELM_ATSDriver::parseParameterList() +{ + if (vo_->os_OK(Teuchos::VERB_LOW)) { + *vo_->os() << "================================================================================" + << std::endl + << "Beginning parseParameterList stage..." << std::endl + << std::flush; + } + // parse my parameter list // domain names - domain_subsurf_ = Keys::readDomain(*plist_, "domain", "domain"); - domain_surf_ = Keys::readDomainHint(*plist_, domain_subsurf_, "subsurface", "surface"); + domain_subsurf_ = Keys::readDomain(*elm_plist_, "domain", "domain"); + domain_surf_ = Keys::readDomainHint(*elm_plist_, domain_subsurf_, "subsurface", "surface"); - // meshes + // meshes mesh_subsurf_ = S_->GetMesh(domain_subsurf_); mesh_surf_ = S_->GetMesh(domain_surf_); - // keys for fields in mesh info - lat_key_ = Keys::readKey(*plist_, domain_surf_, "latitude", "latitude"); - lon_key_ = Keys::readKey(*plist_, domain_surf_, "longitude", "longitude"); - elev_key_ = Keys::readKey(*plist_, domain_surf_, "elevation", "elevation"); - - // soil parameters/properties - base_poro_key_ = Keys::readKey(*plist_, domain_subsurf_, "base porosity", "base_porosity"); - poro_key_ = Keys::readKey(*plist_, domain_subsurf_, "porosity", "porosity"); - perm_key_ = Keys::readKey(*plist_, domain_subsurf_, "permeability", "permeability"); - ch_b_key_ = Keys::readKey(*plist_, domain_subsurf_, "Clapp and Hornberger b", "clapp_horn_b"); - ch_smpsat_key_ = Keys::readKey(*plist_, - domain_subsurf_, - "Clapp and Hornberger soil mafic potential at saturation", - "clapp_horn_smpsat"); - ch_sr_key_ = Keys::readKey( - *plist_, domain_subsurf_, "Clapp and Hornberger residual saturation", "clapp_horn_sr"); + key_map_[ELM::VarID::TIME] = {"time", Tags::NEXT}; + + // parameters + poro_key_ = Keys::readKey(*elm_plist_, domain_subsurf_, "porosity", "porosity"); + key_map_[ELM::VarID::EFFECTIVE_POROSITY] = { poro_key_, Tags::NEXT }; // potential sources - root_frac_key_ = - Keys::readKey(*plist_, domain_subsurf_, "rooting depth fraction", "rooting_depth_fraction"); - pot_infilt_key_ = - Keys::readKey(*plist_, - domain_surf_, - "potential infiltration mps", - "potential_infiltration_mps"); // inputs onto surface (rain, snowmelt) - pot_evap_key_ = - Keys::readKey(*plist_, domain_surf_, "potential evaporation mps", "potential_evaporation_mps"); - pot_trans_key_ = Keys::readKey( - *plist_, domain_surf_, "potential transpiration mps", "potential_transpiration_mps"); + gross_water_source_key_ = Keys::readKey(*elm_plist_, domain_surf_, "gross water source", "gross_water_source"); + key_map_[ELM::VarID::GROSS_SURFACE_WATER_SOURCE] = { gross_water_source_key_, Tags::NEXT }; + pot_evap_key_ = Keys::readKey(*elm_plist_, domain_surf_, "potential evaporation", "potential_evaporation"); + key_map_[ELM::VarID::POTENTIAL_EVAPORATION] = { pot_evap_key_, Tags::NEXT }; + pot_trans_key_ = Keys::readKey(*elm_plist_, domain_surf_, "potential transpiration mps", "potential_transpiration_mps"); + key_map_[ELM::VarID::POTENTIAL_TRANSPIRATION] = { pot_trans_key_, Tags::NEXT }; // water state - pd_key_ = Keys::readKey(*plist_, domain_surf_, "ponded depth", "ponded_depth"); - wtd_key_ = Keys::readKey(*plist_, domain_surf_, "water table depth", "water_table_depth"); - pres_key_ = Keys::readKey(*plist_, domain_subsurf_, "pressure", "pressure"); - wc_key_ = Keys::readKey(*plist_, domain_subsurf_, "conserved", "water_content"); - pc_key_ = Keys::readKey( - *plist_, domain_subsurf_, "capillary_pressure_gas_liq", "capillary_pressure_gas_liq"); - sat_key_ = Keys::readKey(*plist_, domain_subsurf_, "saturation", "saturation_liquid"); - - // water fluxes - infilt_key_ = - Keys::readKey(*plist_, domain_surf_, "surface-subsurface flux", "surface_subsurface_flux"); - evap_key_ = Keys::readKey(*plist_, domain_surf_, "evaporation", "evaporation"); - trans_key_ = Keys::readKey(*plist_, domain_subsurf_, "transpiration", "transpiration"); + pres_key_ = Keys::readKey(*elm_plist_, domain_subsurf_, "pressure", "pressure"); + key_map_[ELM::VarID::PRESSURE] = { pres_key_, Tags::NEXT }; + sat_key_ = Keys::readKey(*elm_plist_, domain_subsurf_, "saturation liquid", "saturation_liquid"); + key_map_[ELM::VarID::SATURATION_LIQUID] = { sat_key_, Tags::NEXT }; + zwt_key_ = Keys::readKey(*elm_plist_, domain_surf_, "water table depth", "water_table_depth"); + key_map_[ELM::VarID::DEPTH_TO_WATER_TABLE] = { zwt_key_, Tags::NEXT }; + pd_key_ = Keys::readKey(*elm_plist_, domain_surf_, "ponded depth", "ponded_depth"); + key_map_[ELM::VarID::PONDED_DEPTH] = { pd_key_, Tags::NEXT }; + + surf_wc_key_ = Keys::readKey(*elm_plist_, domain_surf_, "surface water content", "water_content"); + key_map_[ELM::VarID::SURFACE_WATER_CONTENT_OLD] = { surf_wc_key_, Tags::CURRENT }; + wc_key_ = Keys::readKey(*elm_plist_, domain_subsurf_, "water content", "water_content"); + key_map_[ELM::VarID::WATER_CONTENT_OLD] = { wc_key_, Tags::CURRENT }; + + lat_key_ = Keys::readKey(*elm_plist_, domain_surf_, "latitude", "latitude"); + lon_key_ = Keys::readKey(*elm_plist_, domain_surf_, "longitude", "longitude"); + + // actual water fluxes + evap_key_ = Keys::readKey(*elm_plist_, domain_surf_, "evaporation", "evaporation"); + key_map_[ELM::VarID::EVAPORATION] = { evap_key_, Tags::NEXT }; + col_trans_key_ = Keys::readKey(*elm_plist_, domain_surf_, "total transpiration", "total_transpiration"); + key_map_[ELM::VarID::TRANSPIRATION] = { col_trans_key_, Tags::NEXT }; + col_baseflow_key_ = Keys::readKey(*elm_plist_, domain_surf_, "baseflow generation", "baseflow_mps"); + key_map_[ELM::VarID::BASEFLOW] = { col_baseflow_key_, Tags::NEXT }; + col_runoff_key_ = Keys::readKey(*elm_plist_, domain_surf_, "runoff generation", "runoff_generation_mps"); + key_map_[ELM::VarID::RUNOFF] = { col_runoff_key_, Tags::NEXT }; // keys for fields used to convert ELM units to ATS units - surf_mol_dens_key_ = - Keys::readKey(*plist_, domain_surf_, "surface molar density", "molar_density_liquid"); - surf_mass_dens_key_ = - Keys::readKey(*plist_, domain_surf_, "surface mass density", "mass_density_liquid"); - subsurf_mol_dens_key_ = - Keys::readKey(*plist_, domain_subsurf_, "molar density", "molar_density_liquid"); - subsurf_mass_dens_key_ = - Keys::readKey(*plist_, domain_subsurf_, "mass density", "mass_density_liquid"); - - // need to put into observations or explicitly update if value other than 0.0 is desired - total_trans_key_ = - Keys::readKey(*plist_, domain_surf_, "total transpiration", "total_transpiration"); + surf_mol_dens_key_ = Keys::readKey(*elm_plist_, domain_surf_, "surface molar density", "molar_density_liquid"); + mol_dens_key_ = Keys::readKey(*elm_plist_, domain_subsurf_, "molar density", "molar_density_liquid"); + // surf_mass_dens_key_ = Keys::readKey(*elm_plist_, domain_surf_, "surface mass density", "mass_density_liquid"); + // subsurf_mass_dens_key_ = Keys::readKey(*elm_plist_, domain_subsurf_, "mass density", "mass_density_liquid"); // cell vol keys surf_cv_key_ = Keys::getKey(domain_surf_, "cell_volume"); cv_key_ = Keys::getKey(domain_subsurf_, "cell_volume"); - // currently unused keys - //pres_atm_key_ = Keys::readKey(*plist_, domain_surf_, "atmospheric_pressure", "atmospheric_pressure"); // hardwired as 101325 - //sat_gas_key_ = Keys::readKey(*plist_, domain_subsurf_, "saturation gas", "saturation_gas"); // probably never needed - //sat_ice_key_ = Keys::readKey(*plist_, domain_subsurf_, "saturation ice", "saturation_ice"); // not until energy - - // -- check that number of surface cells = number of columns - ncolumns_ = mesh_surf_->getNumEntities(AmanziMesh::CELL, AmanziMesh::Parallel_kind::OWNED); - AMANZI_ASSERT(ncolumns_ == mesh_subsurf_->columns.num_columns_owned); - - // -- get num cells per column - include consistency check later need to know - // if coupling zone is the entire subsurface mesh (as currently coded) or - // a portion of the total depth specified by # of cells into the - // subsurface + // -- mesh info + ncolumns = mesh_surf_->getNumEntities(AmanziMesh::CELL, AmanziMesh::Parallel_kind::OWNED); auto col_zero = mesh_subsurf_->columns.getCells(0); ncells_per_col_ = col_zero.size(); - for (int col = 0; col != ncolumns_; ++col) - AMANZI_ASSERT(mesh_subsurf_->columns.getCells(col).size() == ncells_per_col_); + + // require my primary variables + // parameters set by ELM + requireEvaluatorAtNext(poro_key_, Amanzi::Tags::NEXT, *S_, poro_key_); + + // potential fluxes (ELM -> ATS) + requireEvaluatorAtNext(gross_water_source_key_, Amanzi::Tags::NEXT, *S_, gross_water_source_key_); + requireEvaluatorAtNext(pot_evap_key_, Amanzi::Tags::NEXT, *S_, pot_evap_key_); + requireEvaluatorAtNext(pot_trans_key_, Amanzi::Tags::NEXT, *S_, pot_trans_key_); + + Coordinator::parseParameterList(); + if (vo_->os_OK(Teuchos::VERB_LOW)) { + *vo_->os() << " ... completed: "; + reportOneTimer_("2a: parseParameterList"); + } } void ELM_ATSDriver::setup() { - // potential fluxes (ELM -> ATS) - requireEvaluatorAtNext(pot_infilt_key_, Amanzi::Tags::NEXT, *S_, pot_infilt_key_) - .SetMesh(mesh_surf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - requireEvaluatorAtNext(pot_evap_key_, Amanzi::Tags::NEXT, *S_, pot_evap_key_) - .SetMesh(mesh_surf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - requireEvaluatorAtNext(pot_trans_key_, Amanzi::Tags::NEXT, *S_, pot_trans_key_) - .SetMesh(mesh_surf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - - // subsurface properties - requireEvaluatorAtNext(base_poro_key_, Amanzi::Tags::NEXT, *S_, base_poro_key_) - .SetMesh(mesh_subsurf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - requireEvaluatorAtNext(perm_key_, Amanzi::Tags::NEXT, *S_, perm_key_) - .SetMesh(mesh_subsurf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - - // dynamic subsurface properties - requireEvaluatorAtNext(root_frac_key_, Amanzi::Tags::NEXT, *S_, root_frac_key_) - .SetMesh(mesh_subsurf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - requireEvaluatorAtNext(poro_key_, - Amanzi::Tags::NEXT, - *S_) // use base_porosity from elm and ATS model for compressibility - .SetMesh(mesh_subsurf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - - // Clapp and Hornberger water retention params (ELM -> ATS) - requireEvaluatorAtNext(ch_b_key_, Amanzi::Tags::NEXT, *S_, ch_b_key_) - .SetMesh(mesh_subsurf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - requireEvaluatorAtNext(ch_smpsat_key_, Amanzi::Tags::NEXT, *S_, ch_smpsat_key_) - .SetMesh(mesh_subsurf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - requireEvaluatorAtNext(ch_sr_key_, Amanzi::Tags::NEXT, *S_, ch_sr_key_) - .SetMesh(mesh_subsurf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - - // per-column ATS water state (ATS -> ELM) - requireEvaluatorAtNext(pd_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_surf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - requireEvaluatorAtNext(wtd_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_surf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - - // per cell ATS water state - requireEvaluatorAtNext(pc_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_subsurf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - requireEvaluatorAtNext(sat_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_subsurf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - - // mesh data - // requireEvaluatorAtNext(lat_key_, Amanzi::Tags::NEXT, *S_) - // .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - // requireEvaluatorAtNext(lon_key_, Amanzi::Tags::NEXT, *S_) - // .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - // requireEvaluatorAtNext(elev_key_, Amanzi::Tags::NEXT, *S_) - // .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + if (vo_->os_OK(Teuchos::VERB_LOW)) { + *vo_->os() << "================================================================================" + << std::endl + << "Beginning setup stage..." << std::endl + << std::flush; + } + + // variables needed for unit changes requireEvaluatorAtNext(surf_cv_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_surf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); requireEvaluatorAtNext(cv_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_subsurf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); - // may not be necessary? any PK that utilizes this should already have density requireEvaluatorAtNext(surf_mol_dens_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_surf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - requireEvaluatorAtNext(surf_mass_dens_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_surf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - requireEvaluatorAtNext(subsurf_mol_dens_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_subsurf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - requireEvaluatorAtNext(subsurf_mass_dens_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_subsurf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - - requireEvaluatorAtNext(total_trans_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_surf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - requireEvaluatorAtNext(wc_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_subsurf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireEvaluatorAtNext(mol_dens_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + + // ELM --> ATS variables + // -- parameters set by ELM + requireEvaluatorAtNext(poro_key_, Amanzi::Tags::NEXT, *S_, poro_key_) + .SetMesh(mesh_subsurf_)->SetComponent("cell", AmanziMesh::CELL, 1); + + // -- potential fluxes + requireEvaluatorAtNext(gross_water_source_key_, Amanzi::Tags::NEXT, *S_, gross_water_source_key_) + .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); + requireEvaluatorAtNext(pot_evap_key_, Amanzi::Tags::NEXT, *S_, pot_evap_key_) + .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); + requireEvaluatorAtNext(pot_trans_key_, Amanzi::Tags::NEXT, *S_, pot_trans_key_) + .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); + + // -- water contents at the OLD time -- note, these are already primary and + // -- owned by overland flow. + requireEvaluatorAtCurrent(wc_key_, Amanzi::Tags::CURRENT, *S_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireEvaluatorAtCurrent(surf_wc_key_, Amanzi::Tags::CURRENT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + + // ATS --> ELM variables + requireEvaluatorAtNext(pd_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireEvaluatorAtNext(zwt_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireEvaluatorAtNext(sat_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireEvaluatorAtNext(pres_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); requireEvaluatorAtNext(evap_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_surf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - requireEvaluatorAtNext(trans_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_subsurf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - - requireEvaluatorAtNext(infilt_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_surf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - requireEvaluatorAtNext(pres_key_, Amanzi::Tags::NEXT, *S_, "flow") - .SetMesh(mesh_subsurf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireEvaluatorAtNext(col_trans_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireEvaluatorAtNext(col_baseflow_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireEvaluatorAtNext(col_runoff_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + + // mesh info -- optional! + if (S_->HasEvaluatorList(lat_key_)) { + requireEvaluatorAtNext(lat_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireEvaluatorAtNext(lon_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + } + // This must be last always -- it allocates memory calling State::setup, so + // all other setup must be done. Coordinator::setup(); + if (vo_->os_OK(Teuchos::VERB_LOW)) { + *vo_->os() << " ... completed: "; + reportOneTimer_("2b: setup"); + } } // // dz & depth are currently ignored -- presumed identical between ATS & ELM // -void -ELM_ATSDriver::get_mesh_info(int& ncols_local, - int& ncols_global, - double* const lat, - double* const lon, - double* const elev, - double* const surface_area, - int* const pft, - int& nlevgrnd, - double* const depth) +ELM_ATSDriver::MeshInfo ELM_ATSDriver::getMeshInfo() { - ncols_local = ncolumns_; - ncols_global = mesh_surf_->getMap(AmanziMesh::Entity_kind::CELL, false).NumGlobalElements(); - - // copyFromSurf_(elev, elev_key_); - // copyFromSurf_(surface_area, surf_cv_key_); - // copyFromSurf_(lat, lat_key_); - // copyFromSurf_(lon, lon_key_); - - // NOTE: figure out how to map from ATS LC types to ELM PFT... --etc - // hard-coded veg type for now... - for (int i = 0; i != ncolumns_; ++i) pft[i] = 1; - - nlevgrnd = ncells_per_col_; - const auto& cells_in_col = mesh_subsurf_->columns.getCells(0); - double ldepth = 0.; - double top_face_z = mesh_subsurf_->getFaceCentroid(mesh_subsurf_->columns.getFaces(0)[0])[2]; + ELM_ATSDriver::MeshInfo info; + info.ncols_local = ncolumns; + info.ncols_global = mesh_surf_->getMap(AmanziMesh::Entity_kind::CELL, false).NumGlobalElements(); + info.nlevgrnd = ncells_per_col_; + + // compute dzs on one column only -- presumed terrain following! + info.dzs.resize(ncells_per_col_); + const auto& col_cells = mesh_subsurf_->columns.getCells(0); + double surface_area = mesh_surf_->getCellVolume(0); for (int i = 0; i != ncells_per_col_; ++i) { - double bottom_face_z = - mesh_subsurf_->getFaceCentroid(mesh_subsurf_->columns.getFaces(0)[i + 1])[2]; - double dz = top_face_z - bottom_face_z; - ldepth += dz / 2.; - depth[i] = ldepth; - top_face_z = bottom_face_z; - ldepth += dz / 2.; + info.dzs[i] = mesh_subsurf_->getCellVolume(i) / surface_area; } - // hard-coded Toledo OH for now... - for (int i = 0; i != ncolumns_; ++i) { - lat[i] = 41.65; - lon[i] = -83.54; + // surface area + info.areas.resize(ncolumns); + for (int i = 0; i != ncolumns; ++i) { + info.areas[i] = mesh_surf_->getCellVolume(i); } + + // lat/lon - optional + info.latitudes.resize(ncolumns, -1.); + info.longitudes.resize(ncolumns, -1.); + + if (S_->HasEvaluator(lat_key_, Tags::NEXT)) { + S_->GetEvaluator(lat_key_, Tags::NEXT).Update(*S_, "ELM ATS driver"); + S_->GetEvaluator(lon_key_, Tags::NEXT).Update(*S_, "ELM ATS driver"); + const Epetra_MultiVector& lat = *S_->Get(lat_key_, Tags::NEXT) + .ViewComponent("cell", false); + const Epetra_MultiVector& lon = *S_->Get(lon_key_, Tags::NEXT) + .ViewComponent("cell", false); + + for (int i = 0; i != ncolumns; ++i) { + info.latitudes[i] = lat[0][i]; + info.longitudes[i] = lon[0][i]; + } + } + return info; } -void -ELM_ATSDriver::initialize(double t, - double const* const elm_water_content, - double const* const elm_pressure) +void ELM_ATSDriver::initialize() { - // Assign start time to ATS - t0_ = t; - - // initialize ATS data, commit initial conditions - Coordinator::initialize(); + if (vo_->os_OK(Teuchos::VERB_LOW)) { + *vo_->os() << "================================================================================" + << std::endl + << "Beginning initialize stage..." << std::endl + << std::flush; + } - // initialize pressure field - ELM_ATSDriver::init_pressure_from_wc_(elm_water_content); + // some parameters have been initialized in ELM + S_->GetRecordW(poro_key_, Tags::NEXT, poro_key_).set_initialized(); // set as zero and tag as initialized - initZero_(root_frac_key_); - initZero_(pot_infilt_key_); - initZero_(pot_trans_key_); - initZero_(pot_evap_key_); - initZero_(infilt_key_); - initZero_(trans_key_); - initZero_(evap_key_); - initZero_(total_trans_key_); + initValue_(gross_water_source_key_); + initValue_(pot_trans_key_); + initValue_(pot_evap_key_); + + // no current evaluators for these, treat as primary + initValue_(col_baseflow_key_); + initValue_(col_runoff_key_); + + // initialize ATS data, commit initial conditions + Coordinator::initialize(); // visualization at IC -- TODO remove this or place behind flag visualize(); checkpoint(); -} -// use incoming water content to initialize pressure field -void -ELM_ATSDriver::init_pressure_from_wc_(double const* const elm_water_content) -{ - // gravity, atmospheric pressure, and liquid water density - // hardwired for now - const double g = 9.80665; - const double p_atm = 101325.0; - const double rho = 1000.0; - - // evaluators - S_->GetEvaluator(subsurf_mass_dens_key_, Amanzi::Tags::NEXT).Update(*S_, subsurf_mass_dens_key_); - const auto& mass_d = *S_->Get(subsurf_mass_dens_key_, Amanzi::Tags::NEXT) - .ViewComponent("cell", false); - S_->GetEvaluator(poro_key_, Amanzi::Tags::NEXT).Update(*S_, poro_key_); - const auto& por = - *S_->Get(poro_key_, Amanzi::Tags::NEXT).ViewComponent("cell", false); - S_->GetEvaluator(cv_key_, Amanzi::Tags::NEXT).Update(*S_, cv_key_); - const auto& volume = - *S_->Get(cv_key_, Amanzi::Tags::NEXT).ViewComponent("cell", false); - S_->GetEvaluator(surf_cv_key_, Amanzi::Tags::NEXT).Update(*S_, surf_cv_key_); - const auto& area = - *S_->Get(surf_cv_key_, Amanzi::Tags::NEXT).ViewComponent("cell", false); - - // writable to pressure - auto& pres = - *S_->GetW(pres_key_, Amanzi::Tags::NEXT, "flow").ViewComponent("cell", false); - - // WRM model - auto& wrm_eval = S_->GetEvaluator(sat_key_, Amanzi::Tags::NEXT); - auto wrm_ptr = dynamic_cast(&wrm_eval); - AMANZI_ASSERT(wrm_ptr != nullptr); - auto wrms_ = wrm_ptr->get_WRMs(); - AMANZI_ASSERT(wrms_->second.size() == 1); // only supports one WRM for now - Teuchos::RCP wrm_ = wrms_->second[0]; - - // initialize pressure field from ELM water content - // per-column hydrostatic pressure in areas of continuous total saturation - // unsaturated areas are considered to be in contact with atmosphere - for (int i = 0; i != ncolumns_; ++i) { - const auto& cells_of_col = mesh_subsurf_->columns.getCells(i); - int top_sat_idx = -1; - double sat_depth = 0.0; - for (int j = 0; j != ncells_per_col_; ++j) { - // convert ELM water content (kg/m2] to saturation of pore space (0 to 1) [-] - // VWC = elm_wc * 1/dz * 1/porosity * 1/mass density - // [-] = [kg/m2] * [m^-1] * [-] * [m3/kg] - const double dz = volume[0][cells_of_col[j]] / area[0][i]; - const double factor = 1 / (dz * por[0][cells_of_col[j]] * mass_d[0][cells_of_col[j]]); - const double satl = elm_water_content[j * ncolumns_ + i] * factor; - if (satl < 1.0) { - pres[0][cells_of_col[j]] = p_atm - wrm_->capillaryPressure(satl); - top_sat_idx = -1; - } else { - if (top_sat_idx == -1) { - top_sat_idx = j; - sat_depth = 0.0; - } - sat_depth += dz; - pres[0][cells_of_col[j]] = p_atm + rho * g * (sat_depth - dz / 2); - } - } + if (vo_->os_OK(Teuchos::VERB_LOW)) { + *vo_->os() << " ... completed: "; + reportOneTimer_("3: initialize"); } - - // mark pressure as changed and update face values - changedEvaluatorPrimary(pres_key_, Amanzi::Tags::NEXT, *S_); - DeriveFaceValuesFromCellValues(S_->GetW(pres_key_, Amanzi::Tags::NEXT, "flow")); - S_->GetRecordW(pres_key_, Amanzi::Tags::NEXT, "flow").set_initialized(); - - // update saturation and water content - S_->GetEvaluator(sat_key_, Amanzi::Tags::NEXT).Update(*S_, sat_key_); - S_->GetEvaluator(wc_key_, Amanzi::Tags::NEXT).Update(*S_, wc_key_); } -void -ELM_ATSDriver::advance(double dt, bool do_chkp, bool do_vis) + +void ELM_ATSDriver::advance(double dt, bool force_chkp, bool force_vis) { - double dt_subcycle = dt; - double t_end = S_->get_time() + dt_subcycle; + { + Teuchos::TimeMonitor timer(*timers_.at("4: solve")); + + AMANZI_ASSERT(std::abs(S_->get_time(Amanzi::Tags::NEXT) - S_->get_time(Amanzi::Tags::CURRENT)) < + 1.e-4); + + if (vo_->os_OK(Teuchos::VERB_LOW)) { + *vo_->os() + << "================================================================================" + << std::endl + << std::endl + << "Cycle = " << S_->get_cycle() << ", Time [days] = " << std::setprecision(16) + << S_->get_time() / (60 * 60 * 24) << ", dt [days] = " << std::setprecision(16) + << dt / (60 * 60 * 24) << std::endl + << "--------------------------------------------------------------------------------" + << std::endl; + } + S_->Assign("dt", Amanzi::Tags::DEFAULT, "dt", dt); + S_->advance_time(Amanzi::Tags::NEXT, dt); + + // old water_contents are set by ELM, using ELM units. Change them to ATS units + // + // NOTE: the change in units is done here, and not in evaluators like + // inputs (e.g. transpiration_mps), because WC evaluators are expected to be + // primary by flow PKs, and then get overwritten here. + // + // NOTE: these were just marked as changed in a call to getFieldPtrW so no + // need to mark them as changed again. + // -- subsurface wc in volumetric water content --> mols + { + auto& wc = *S_->GetW(wc_key_, Tags::CURRENT, + S_->GetRecord(wc_key_, Tags::CURRENT).owner()).ViewComponent("cell", false); + const auto& cv = *S_->Get(cv_key_, Tags::NEXT).ViewComponent("cell", false); + const auto& n_liq = *S_->Get(mol_dens_key_, Tags::NEXT).ViewComponent("cell", false); + for (int c = 0; c != wc.MyLength(); ++c) { + wc[0][c] = wc[0][c] * cv[0][c] * n_liq[0][c]; + } + } - bool fail{ false }; - while (S_->get_time() < t_end && dt_subcycle > 0.0) { - S_->Assign("dt", Amanzi::Tags::DEFAULT, "dt", dt_subcycle); - S_->advance_time(Amanzi::Tags::NEXT, dt_subcycle); + // -- surface wc in m ponded depth --> mols + { + auto& surf_wc = *S_->GetW(surf_wc_key_, Tags::CURRENT, + S_->GetRecord(surf_wc_key_, Tags::CURRENT).owner()).ViewComponent("cell", false); + const auto& surf_cv = *S_->Get(surf_cv_key_, Tags::NEXT).ViewComponent("cell", false); + const auto& surf_n_liq = *S_->Get(surf_mol_dens_key_, Tags::NEXT).ViewComponent("cell", false); + for (int c = 0; c != surf_wc.MyLength(); ++c) { + surf_wc[0][c] = surf_wc[0][c] * surf_cv[0][c] * surf_n_liq[0][c]; + } + } // solve model for a single timestep - fail = Coordinator::advance(); - + bool fail = Coordinator::advance(); if (fail) { - // reset t_new - S_->set_time(Amanzi::Tags::NEXT, S_->get_time(Amanzi::Tags::CURRENT)); - } else { - S_->set_time(Amanzi::Tags::CURRENT, S_->get_time(Amanzi::Tags::NEXT)); - S_->advance_cycle(); - - // make observations, vis, and checkpoints - for (const auto& obs : observations_) obs->MakeObservations(S_.ptr()); - - // vis/checkpoint if EITHER ATS or ELM request it - //if (do_vis && !visualize()) visualize(true); - //if (do_chkp && !checkpoint()) checkpoint(true); - - visualize(do_vis); - checkpoint(do_chkp); + Errors::Message msg("ELM_ATSDriver: advance(dt) failed. Make ATS subcycle for proper ELM use."); + Exceptions::amanzi_throw(msg); } - dt_subcycle = Coordinator::get_dt(fail); - } // while + S_->set_time(Amanzi::Tags::CURRENT, S_->get_time(Amanzi::Tags::NEXT)); + S_->advance_cycle(); - if (fail) { - // write one more vis for help debugging - S_->advance_cycle(Amanzi::Tags::NEXT); - visualize(true); // force vis - - // flush observations to make sure they are saved - for (const auto& obs : observations_) obs->Flush(); - - // dump a post_mortem checkpoint file for debugging - checkpoint_->set_filebasename("post_mortem"); - checkpoint_->Write(*S_, Checkpoint::WriteType::POST_MORTEM); - - Errors::Message msg("ELM_ATSDriver: advance(dt) failed."); - Exceptions::amanzi_throw(msg); + visualize(force_vis); + checkpoint(force_chkp); + observe(); + } + if (vo_->os_OK(Teuchos::VERB_LOW)) { + *vo_->os() << " ... completed: "; + reportOneTimer_("4: solve"); } -} // advance() +} // simulates external timeloop with dt coming from calling model -void -ELM_ATSDriver::advance_test() +void ELM_ATSDriver::advanceTest() { while (S_->get_time() < t1_) { // use dt from ATS for testing @@ -489,301 +411,62 @@ ELM_ATSDriver::advance_test() } } -void -ELM_ATSDriver::finalize() -{ - WriteStateStatistics(*S_, *vo_); - report_memory(); - Teuchos::TimeMonitor::summarize(*vo_->os()); - Coordinator::finalize(); -} - - -void -ELM_ATSDriver::set_soil_hydrologic_parameters(double const* const base_porosity, - double const* const hydraulic_conductivity, - double const* const clapp_horn_b, - double const* const clapp_horn_smpsat, - double const* const clapp_horn_sr) -{ - // convert Ksat to perm via rho * g / visc using default rho and visc values. - copyToSub_(hydraulic_conductivity, perm_key_); - double factor = 8.9e-4 / (1000 * 9.80665); - S_->GetW(perm_key_, Amanzi::Tags::NEXT, perm_key_).Scale(factor); - S_->GetRecordW(perm_key_, Amanzi::Tags::NEXT, perm_key_).set_initialized(); - - copyToSub_(base_porosity, base_poro_key_); - copyToSub_(clapp_horn_b, ch_b_key_); - copyToSub_(clapp_horn_smpsat, ch_smpsat_key_); - copyToSub_(clapp_horn_sr, ch_sr_key_); - - S_->GetRecordW(base_poro_key_, Amanzi::Tags::NEXT, base_poro_key_).set_initialized(); - S_->GetRecordW(ch_b_key_, Amanzi::Tags::NEXT, ch_b_key_).set_initialized(); - S_->GetRecordW(ch_smpsat_key_, Amanzi::Tags::NEXT, ch_smpsat_key_).set_initialized(); - S_->GetRecordW(ch_sr_key_, Amanzi::Tags::NEXT, ch_sr_key_).set_initialized(); -} - - -void -ELM_ATSDriver::set_veg_parameters(double const* const mafic_potential_full_turgor, - double const* const mafic_potential_wilt_point) -{ - // pass for now! FIXME --etc -} - - -void -ELM_ATSDriver::set_soil_hydrologic_properties(double const* const effective_porosity) -{ - // this isn't well defined -- pass for now --etc -} - - -void -ELM_ATSDriver::set_veg_properties(double const* const rooting_fraction) -{ - copyToSub_(rooting_fraction, root_frac_key_); -} - - -void -ELM_ATSDriver::set_potential_sources(double const* const elm_surface_input, - double const* const elm_evaporation, - double const* const elm_transpiration) +void ELM_ATSDriver::finalize() { - // ELM's infiltration, evaporation, and transpiration are in units of mm s^-1 - // ATS units are in mol m^-2 s^-1 - - // kg / m2 / s * m3/kg * mol/m3 = mol m^-2 s^-1 - // or - // mm / s * m/mm * mol/m3 = mol m^-2 s^-1 - auto& in = *S_->GetW(pot_infilt_key_, Amanzi::Tags::NEXT, pot_infilt_key_) - .ViewComponent("cell", false); - auto& ev = *S_->GetW(pot_evap_key_, Amanzi::Tags::NEXT, pot_evap_key_) - .ViewComponent("cell", false); - auto& tr = *S_->GetW(pot_trans_key_, Amanzi::Tags::NEXT, pot_trans_key_) - .ViewComponent("cell", false); - AMANZI_ASSERT(in.MyLength() == ncolumns_); - AMANZI_ASSERT(ev.MyLength() == ncolumns_); - AMANZI_ASSERT(tr.MyLength() == ncolumns_); - - S_->GetEvaluator(surf_mol_dens_key_, Amanzi::Tags::NEXT).Update(*S_, surf_mol_dens_key_); - const auto& surf_d = - *S_->Get(surf_mol_dens_key_, Amanzi::Tags::NEXT).ViewComponent("cell", false); - - for (int i = 0; i != ncolumns_; ++i) { - const double factor = 0.001 * surf_d[0][i]; - in[0][i] = elm_surface_input[i] * factor; - ev[0][i] = elm_evaporation[i] * factor; - tr[0][i] = elm_transpiration[i] * factor; + if (vo_->os_OK(Teuchos::VERB_LOW)) { + *vo_->os() << "================================================================================" + << std::endl + << "Beginning finalize stage..." << std::endl + << std::flush; } - - changedEvaluatorPrimary(pot_infilt_key_, Amanzi::Tags::NEXT, *S_); - changedEvaluatorPrimary(pot_evap_key_, Amanzi::Tags::NEXT, *S_); - changedEvaluatorPrimary(pot_trans_key_, Amanzi::Tags::NEXT, *S_); -} - - -void -ELM_ATSDriver::get_waterstate(double* const ponded_depth, - double* const wt_depth, - double* const soil_water_potential, - double* const matric_potential, - double* const sat_liq, - double* const sat_ice) // remove? -{ - // convert saturation into [kg/m2] from [-] - S_->GetEvaluator(sat_key_, Amanzi::Tags::NEXT).Update(*S_, sat_key_); - const auto& satl = - *S_->Get(sat_key_, Amanzi::Tags::NEXT).ViewComponent("cell", false); - - S_->GetEvaluator(poro_key_, Amanzi::Tags::NEXT).Update(*S_, poro_key_); - const auto& por = - *S_->Get(poro_key_, Amanzi::Tags::NEXT).ViewComponent("cell", false); - - S_->GetEvaluator(subsurf_mass_dens_key_, Amanzi::Tags::NEXT).Update(*S_, subsurf_mass_dens_key_); - const auto& dens = *S_->Get(subsurf_mass_dens_key_, Amanzi::Tags::NEXT) - .ViewComponent("cell", false); - - // TODO look into ELM effective porosity, ATS ice density, ice saturation - for (int i = 0; i != ncolumns_; ++i) { - const auto faces = mesh_subsurf_->columns.getFaces(i); - const auto cells = mesh_subsurf_->columns.getCells(i); - for (int j = 0; j != ncells_per_col_; ++j) { - const double dz = mesh_subsurf_->getFaceCentroid(faces[j])[2] - - mesh_subsurf_->getFaceCentroid(faces[j + 1])[2]; - sat_liq[j * ncolumns_ + i] = satl[0][cells[j]] * por[0][cells[j]] * dens[0][cells[j]] * dz; - } + Coordinator::finalize(); + if (vo_->os_OK(Teuchos::VERB_LOW)) { + *vo_->os() << " ... completed: "; + reportOneTimer_("5: finalize"); } - - // Ponded depth - // convert ATS m to mm - S_->GetEvaluator(pd_key_, Amanzi::Tags::NEXT).Update(*S_, "ELM"); - const auto& pd = - *S_->Get(pd_key_, Amanzi::Tags::NEXT).ViewComponent("cell", false); - for (int i = 0; i != ncolumns_; ++i) ponded_depth[i] = pd[0][i] * 1000.0; - - // water table depth - S_->GetEvaluator(wtd_key_, Amanzi::Tags::NEXT).Update(*S_, "ELM"); - copyFromSurf_(wt_depth, wtd_key_); - - // Soil matric potential - // convert ATS Pa to -mmH2O - // int z_index = mesh_subsurf_->space_dimension() - 1; - // const auto& gravity = S_->Get("gravity", Amanzi::Tags::DEFAULT); - // const double g_inv = 1.0 / gravity[z_index]; // should be -9.80665 m s^-2 - // - // S_->GetEvaluator(pres_key_, Amanzi::Tags::NEXT).Update(*S_, pres_key_); - // const auto& pres = *S_->Get(pres_key_, Amanzi::Tags::NEXT) - // .ViewComponent("cell", false); - // - // S_->GetEvaluator(pc_key_, Amanzi::Tags::NEXT).Update(*S_, "ELM"); - // const auto& pc = *S_->Get(pc_key_, Amanzi::Tags::NEXT) - // .ViewComponent("cell", false); - // - // for (int i=0; i!=ncolumns_; ++i) { - // const auto& cells = mesh_subsurf_->columns.getCells(i); - // for (int j=0; j!=ncells_per_col_; ++j) { - // matric_potential[j * ncolumns_ + i] = pc[0][cells[j]] * g_inv; - // soil_water_potential[j * ncolumns_ + i] = 0.101325 - 1.e-6 * pres[0][cells[j]]; - // } - // } } -void -ELM_ATSDriver::get_water_fluxes(double* const surf_subsurf_flx, - double* const evaporation, - double* const transpiration, - double* const net_subsurface_fluxes, - double* const net_runon) +void ELM_ATSDriver::initValue_(const Key& key, double value) { - // Convert and send ATS fluxes to ELM - - // Flux units: ATS ELM - // surface-evaporation [mol / m2 / s] [mmH2o/s] - // surface-infiltration [mol / m2 / s] [mmH2o/s] - // transpiration [mol / m3 / s] [mmH2o/s] - - // Surface fluxes - S_->GetEvaluator(infilt_key_, Amanzi::Tags::NEXT).Update(*S_, infilt_key_); - const auto& infilt = - *S_->Get(infilt_key_, Amanzi::Tags::NEXT).ViewComponent("cell", false); - - S_->GetEvaluator(evap_key_, Amanzi::Tags::NEXT).Update(*S_, evap_key_); - const auto& evap = - *S_->Get(evap_key_, Amanzi::Tags::NEXT).ViewComponent("cell", false); - - S_->GetEvaluator(surf_mol_dens_key_, Amanzi::Tags::NEXT).Update(*S_, surf_mol_dens_key_); - const auto& surfdens = - *S_->Get(surf_mol_dens_key_, Amanzi::Tags::NEXT).ViewComponent("cell", false); - - // convert mol/m2/s to mmH2O/s for ELM - // mol/m2/s * m3/mol = m/s * mm/m = mm/s - for (int i = 0; i != ncolumns_; ++i) { - const double mm_per_mol = 1000.0 / surfdens[0][i]; - surf_subsurf_flx[i] = infilt[0][i] * mm_per_mol; - evaporation[i] = evap[0][i] * mm_per_mol; - } - - // Subsurface flux - S_->GetEvaluator(trans_key_, Amanzi::Tags::NEXT).Update(*S_, "ELM"); - const auto& trans = - *S_->Get(trans_key_, Amanzi::Tags::NEXT).ViewComponent("cell", false); - - S_->GetEvaluator(subsurf_mol_dens_key_, Amanzi::Tags::NEXT).Update(*S_, subsurf_mol_dens_key_); - const auto& subsurfdens = *S_->Get(subsurf_mol_dens_key_, Amanzi::Tags::NEXT) - .ViewComponent("cell", false); - - // convert mol/m3/s to mmH2O/s by integrating over dz - NO? - // treat the same as surface fluxes? - for (int i = 0; i != ncolumns_; ++i) { - const auto& faces = mesh_subsurf_->columns.getFaces(i); - const auto& cells = mesh_subsurf_->columns.getCells(i); - for (int j = 0; j != ncells_per_col_; ++j) { - double dz = mesh_subsurf_->getFaceCentroid(faces[j])[2] - - mesh_subsurf_->getFaceCentroid(faces[j + 1])[2]; - AMANZI_ASSERT(dz > 0.); - const double factor = dz * 1000.0 / subsurfdens[0][cells[j]]; - transpiration[j * ncolumns_ + i] = trans[0][cells[j]] * factor; - } - } - // unclear how to implement net_subsurface_fluxes and net_runon... --ETC + auto& vec = S_->GetW(key, Tags::NEXT, key); + vec.PutScalar(value); + S_->GetRecordW(key, Tags::NEXT, key).set_initialized(); } -void -ELM_ATSDriver::initZero_(const Key& key) +void ELM_ATSDriver::setScalar(const ELM::VarID& scalar_id, double in) { - auto& vec = S_->GetW(key, Amanzi::Tags::NEXT, key); - vec.PutScalar(0.); - S_->GetRecordW(key, Amanzi::Tags::NEXT, key).set_initialized(); + KeyTag scalar_key = key_map_.at(scalar_id); + S_->GetW(scalar_key.first, scalar_key.second, scalar_key.first) = in; } - -void -ELM_ATSDriver::copyToSurf_(double const* const in, const Key& key, Key owner) +double ELM_ATSDriver::getScalar(const ELM::VarID& scalar_id) { - if (owner.empty() ) owner = key; - // surf maps directly into columns - auto& vec = - *S_->GetW(key, Amanzi::Tags::NEXT, owner).ViewComponent("cell", false); - AMANZI_ASSERT(vec.MyLength() == ncolumns_); - - for (int i = 0; i != ncolumns_; ++i) vec[0][i] = in[i]; - - changedEvaluatorPrimary(key, Amanzi::Tags::NEXT, *S_); + KeyTag scalar_key = key_map_.at(scalar_id); + return S_->Get(scalar_key.first, scalar_key.second); } -// -// ELM data is defined as var(col,lev), meaning that, due to Fortran ordering, -// the column is fasted varying, not the grid cell. -// -void -ELM_ATSDriver::copyToSub_(double const* const in, const Key& key, Key owner) +double const * +ELM_ATSDriver::getFieldPtr(const ELM::VarID& var_id) { - if (owner.empty() ) owner = key; - auto& vec = - *S_->GetW(key, Amanzi::Tags::NEXT, owner).ViewComponent("cell", false); - - for (int i = 0; i != ncolumns_; ++i) { - const auto& cells = mesh_subsurf_->columns.getCells(i); - for (int j = 0; j != ncells_per_col_; ++j) { - vec[0][cells[j]] = in[j * ncolumns_ + i]; - } + Amanzi::KeyTag var_key = key_map_.at(var_id); + if (S_->HasEvaluator(var_key.first, var_key.second)) { + S_->GetEvaluator(var_key.first, var_key.second).Update(*S_, std::string("elm_ats_driver_on_"+domain_subsurf_)); } - changedEvaluatorPrimary(key, Amanzi::Tags::NEXT, *S_); -} - - -void -ELM_ATSDriver::copyFromSurf_(double* const out, const Key& key) const -{ - S_->GetEvaluator(key, Amanzi::Tags::NEXT).Update(*S_, "ELM"); - const auto& vec = *S_->Get(key, Amanzi::Tags::NEXT).ViewComponent("cell", false); - - for (int i = 0; i != ncolumns_; ++i) out[i] = vec[0][i]; + return &(*S_->Get(var_key.first, var_key.second).ViewComponent("cell", false))[0][0]; } - -// -// ELM data is defined as var(col,lev), meaning that, due to Fortran ordering, -// the column is fasted varying, not the grid cell. -// -void -ELM_ATSDriver::copyFromSub_(double* const out, const Key& key) const +double * +ELM_ATSDriver::getFieldPtrW(const ELM::VarID& var_id) { - S_->GetEvaluator(key, Amanzi::Tags::NEXT).Update(*S_, "ELM"); - const auto& vec = *S_->Get(key, Amanzi::Tags::NEXT).ViewComponent("cell", false); + Amanzi::KeyTag var_key = key_map_.at(var_id); + Amanzi::changedEvaluatorPrimary(var_key.first, var_key.second, *S_); - for (int i = 0; i != ncolumns_; ++i) { - const auto cells = mesh_subsurf_->columns.getCells(i); - for (int j = 0; j != ncells_per_col_; ++j) { - out[j * ncolumns_ + i] = vec[0][cells[j]]; - } - } + Amanzi::Key owner = S_->GetRecord(var_key.first, var_key.second).owner(); + return &(*S_->GetW(var_key.first, var_key.second, owner).ViewComponent("cell", false))[0][0]; } diff --git a/src/executables/elm_ats_api/elm_ats_driver.hh b/src/executables/elm_ats_api/elm_ats_driver.hh index 910f45ef2..72bb8ce21 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.hh +++ b/src/executables/elm_ats_api/elm_ats_driver.hh @@ -5,30 +5,6 @@ Authors: Joe Beisman */ -//! Simulation control from ELM. - -/*! - -The expected order of evaluation is: - -..code:: - - ELM_ATSDriver ats; - ats.setup(); - ats.get_mesh_info(); - ats.set_soil_parameters(); - ats.set_veg_parameters(); - ats.initialize(); - - for (step) { - ats.set_soil_properties(); - ats.set_veg_properties(); - ats.advance(); - ats.get_water_state(); - ats.get_water_fluxes(); - } - -*/ #pragma once @@ -40,76 +16,55 @@ The expected order of evaluation is: #include "Key.hh" #include "coordinator.hh" +#include "ats_variables.hh" + namespace ATS { using namespace Amanzi; class ELM_ATSDriver : public Coordinator { + public: + + struct MeshInfo { + int ncols_local; + int ncols_global; + int nlevgrnd; + std::vector dzs; + std::vector areas; + std::vector latitudes; + std::vector longitudes; + }; + ELM_ATSDriver(const Teuchos::RCP& plist, const Teuchos::RCP& wallclock_timer, const Teuchos::RCP>& teuchos_comm, const Amanzi::Comm_ptr_type& comm, + const std::string& logfile, int npfts = 17); - void finalize(); - void advance(double dt, bool checkpoint, bool vis); - void advance_test(); - - void get_mesh_info(int& ncols_local, - int& ncols_global, - double* const lat, - double* const lon, - double* const elev, - double* const surf_area, - int* const pft, - int& nlevgrnd, - double* const depth); - + MeshInfo getMeshInfo(); + void parseParameterList(); void setup(); + void initialize(); + void advance(double dt, bool checkpoint, bool vis); + void advanceTest(); + void finalize(); - void initialize(double t, double const* const p_atm, double const* const pressure); - - void set_soil_hydrologic_parameters(double const* const base_porosity, - double const* const hydraulic_conductivity, - double const* const clapp_horn_b, - double const* const clapp_horn_smpsat, - double const* const clapp_horn_sr); - void set_veg_parameters(double const* const mafic_potential_full_turgor, - double const* const mafic_potential_wilt_point); - void set_soil_hydrologic_properties(double const* const effective_porosity); - void set_veg_properties(double const* const rooting_fraction); - void set_potential_sources(double const* const elm_surface_input, - double const* const elm_evaporation, - double const* const elm_transpiration); - - void get_waterstate(double* const surface_ponded_depth, - double* const water_table_depth, - double* const soil_pressure, - double* const soil_psi, - double* const sat_liq, - double* const sat_ice); - - void get_water_fluxes(double* const soil_infiltration, - double* const evaporation, - double* const transpiration, - double* const net_subsurface_fluxes, - double* const net_runon); + void setScalar(const ELM::VarID& scalar_id, double in); + double getScalar(const ELM::VarID& scalar_id); - private: - void init_pressure_from_wc_(double const* const elm_water_content); + double const * getFieldPtr(const ELM::VarID& var_id); + double * getFieldPtrW(const ELM::VarID& var_id); - void copyToSurf_(double const* const in, const Key& key, Key owner = ""); - void copyToSub_(double const* const in, const Key& key, Key owner = ""); - void copyFromSurf_(double* const out, const Key& key) const; - void copyFromSub_(double* const out, const Key& key) const; - void initZero_(const Key& key); + int ncolumns; + private: + void initValue_(const Key& key, double value = 0.); private: - Teuchos::RCP elm_list_; + Teuchos::RCP elm_plist_; - int ncolumns_; int ncells_per_col_; int npfts_; @@ -119,43 +74,48 @@ class ELM_ATSDriver : public Coordinator { Teuchos::RCP mesh_subsurf_; Teuchos::RCP mesh_surf_; - Key lat_key_; - Key lon_key_; - Key elev_key_; - Key base_poro_key_; - Key perm_key_; - Key ch_b_key_; - Key ch_smpsat_key_; - Key ch_sr_key_; Key poro_key_; - Key root_frac_key_; + + Key gross_water_source_key_; Key pot_evap_key_; Key pot_trans_key_; - Key pot_infilt_key_; - Key pd_key_; - Key wtd_key_; - Key pres_key_; + Key wc_key_; - Key pc_key_; + Key surf_wc_key_; + + Key pres_key_; Key sat_key_; - //Key sat_gas_key_; - //Key sat_ice_key_; - Key infilt_key_; - Key trans_key_; + Key pd_key_; + Key zwt_key_; + Key evap_key_; - Key total_trans_key_; + Key col_trans_key_; + Key col_baseflow_key_; + Key col_runoff_key_; + Key surf_mol_dens_key_; - Key surf_mass_dens_key_; - Key subsurf_mol_dens_key_; - Key subsurf_mass_dens_key_; + Key mol_dens_key_; + // Key surf_mass_dens_key_; + // Key subsurf_mass_dens_key_; + Key surf_cv_key_; Key cv_key_; + Key lat_key_; + Key lon_key_; + + std::map key_map_; + }; // // Nonmember constructor/factory reads file, converts comm to the right type. // -ELM_ATSDriver* createELM_ATSDriver(MPI_Fint* f_comm, const char* infile, int npfts = 17); +ELM_ATSDriver* +createELM_ATSDriver(MPI_Fint *f_comm, const char *infile, const char *logfile, int npfts=17); } // namespace ATS + + + + diff --git a/src/pks/flow/overland_pressure_pk.cc b/src/pks/flow/overland_pressure_pk.cc index 4b69ac95e..0937479bd 100644 --- a/src/pks/flow/overland_pressure_pk.cc +++ b/src/pks/flow/overland_pressure_pk.cc @@ -124,6 +124,7 @@ OverlandPressureFlow::parseParameterList() // require a few primary variable keys now to set the leaf node in the dep graph requireEvaluatorAtCurrent(pd_key_, tag_current_, *S_, name_); + requireEvaluatorAtNext(flux_key_, tag_next_, *S_, name_); } diff --git a/src/pks/mpc/mpc_coupled_water.cc b/src/pks/mpc/mpc_coupled_water.cc index 6eb2644ab..6ef1bb08e 100644 --- a/src/pks/mpc/mpc_coupled_water.cc +++ b/src/pks/mpc/mpc_coupled_water.cc @@ -45,6 +45,9 @@ MPCCoupledWater::parseParameterList() exfilt_key_ = Keys::readKey(*plist_, domain_surf_, "exfiltration flux", "surface_subsurface_flux"); + // require the primary variable coupling field, claim ownership + requireEvaluatorAtNext(exfilt_key_, tag_next_, *S_, name_); + StrongMPC::parseParameterList(); } @@ -60,14 +63,14 @@ MPCCoupledWater::Setup() domain_flow_pk_ = sub_pks_[0]; surf_flow_pk_ = sub_pks_[1]; - // call the MPC's setup, which calls the sub-pk's setups - StrongMPC::Setup(); - - // require the coupling fields, claim ownership + // set structure of the primary variable requireEvaluatorAtNext(exfilt_key_, tag_next_, *S_, name_) .SetMesh(surf_mesh_) ->SetComponent("cell", AmanziMesh::Entity_kind::CELL, 1); + // call the MPC's setup, which calls the sub-pk's setups + StrongMPC::Setup(); + // Create the preconditioner. // -- collect the preconditioners precon_ = domain_flow_pk_->preconditioner(); @@ -163,7 +166,7 @@ MPCCoupledWater::FunctionalResidual(double t_old, t_old, t_new, u_old->SubVector(1), u_new->SubVector(1), g->SubVector(1)); // The residual of the surface flow equation provides the water flux from - // subsurface to surface. + // subsurface to surface, in [mols s^-1] Epetra_MultiVector& source = *S_->GetW(exfilt_key_, tag_next_, name_).ViewComponent("cell", false); source = *g->SubVector(1)->Data()->ViewComponent("cell", false); diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.cc b/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.cc index 57920c337..943de30cb 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.cc +++ b/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.cc @@ -20,7 +20,7 @@ namespace SurfaceBalance { double readPositiveLandCoverParameter(Teuchos::ParameterList& plist, const std::string& name) { - double res = plist.get(name, NAN); + double res = plist.get(name); if (res < 0) { Errors::Message msg; msg << "Invalid land cover parameter \"" << name << "\" in land cover type \"" << plist.name() @@ -33,7 +33,7 @@ readPositiveLandCoverParameter(Teuchos::ParameterList& plist, const std::string& double readNegativeLandCoverParameter(Teuchos::ParameterList& plist, const std::string& name) { - double res = plist.get(name, NAN); + double res = plist.get(name); if (res > 0) { Errors::Message msg; msg << "Invalid land cover parameter \"" << name << "\" in land cover type \"" << plist.name() @@ -58,35 +58,69 @@ readZeroOneLandCoverParameter(Teuchos::ParameterList& plist, const std::string& LandCover::LandCover(Teuchos::ParameterList& plist) - : rooting_depth_max(readPositiveLandCoverParameter(plist, "rooting depth max [m]")), - rooting_profile_alpha(readPositiveLandCoverParameter(plist, "rooting profile alpha [-]")), - rooting_profile_beta(plist.get("rooting profile beta [-]", NAN)), - stomata_closed_capillary_pressure( - readPositiveLandCoverParameter(plist, "capillary pressure at fully closed stomata [Pa]")), - stomata_open_capillary_pressure( - readPositiveLandCoverParameter(plist, "capillary pressure at fully open stomata [Pa]")), - maximum_xylem_capillary_pressure( - readPositiveLandCoverParameter(plist, "maximum xylem capillary pressure [Pa]")), - leaf_on_doy(plist.get("leaf on time [doy]", NAN)), - leaf_off_doy(plist.get("leaf off time [doy]", NAN)), - pt_alpha_snow(readPositiveLandCoverParameter(plist, "Priestley-Taylor alpha of snow [-]")), - pt_alpha_canopy(readPositiveLandCoverParameter(plist, "Priestley-Taylor alpha of canopy [-]")), - pt_alpha_ground( - readPositiveLandCoverParameter(plist, "Priestley-Taylor alpha of bare ground [-]")), - pt_alpha_transpiration( - readPositiveLandCoverParameter(plist, "Priestley-Taylor alpha of transpiration [-]")), - albedo_ground(readZeroOneLandCoverParameter(plist, "albedo of bare ground [-]")), - emissivity_ground(readZeroOneLandCoverParameter(plist, "emissivity of bare ground [-]")), - albedo_canopy(readZeroOneLandCoverParameter(plist, "albedo of canopy [-]")), - beers_k_sw( - readPositiveLandCoverParameter(plist, "Beer's law extinction coefficient, shortwave [-]")), - beers_k_lw( - readPositiveLandCoverParameter(plist, "Beer's law extinction coefficient, longwave [-]")), - snow_transition_depth(readPositiveLandCoverParameter(plist, "snow transition depth [m]")), - water_transition_depth(readPositiveLandCoverParameter(plist, "water transition depth [m]")), - roughness_ground(readPositiveLandCoverParameter(plist, "roughness length of bare ground [m]")), - roughness_snow(readPositiveLandCoverParameter(plist, "roughness length of snow [m]")), - mannings_n(readPositiveLandCoverParameter(plist, "Manning's n [?]")) + : rooting_depth_max(plist.isParameter("rooting depth max [m]") ? + readPositiveLandCoverParameter(plist, "rooting depth max [m]") : + Relations::NaN), + rooting_profile_alpha(plist.isParameter("rooting profile alpha [-]") ? + readPositiveLandCoverParameter(plist, "rooting profile alpha [-]") : + Relations::NaN), + rooting_profile_beta(plist.isParameter("rooting profile beta [-]") ? + readPositiveLandCoverParameter(plist, "rooting profile beta [-]") : + Relations::NaN), + stomata_closed_capillary_pressure(plist.isParameter("capillary pressure at fully closed stomata [Pa]") ? + readPositiveLandCoverParameter(plist, "capillary pressure at fully closed stomata [Pa]") : + Relations::NaN), + stomata_open_capillary_pressure(plist.isParameter("capillary pressure at fully open stomata [Pa]") ? + readPositiveLandCoverParameter(plist, "capillary pressure at fully open stomata [Pa]") : + Relations::NaN), + maximum_xylem_capillary_pressure(plist.isParameter("maximum xylem capillary pressure [Pa]") ? + readPositiveLandCoverParameter(plist, "maximum xylem capillary pressure [Pa]") : + Relations::NaN), + leaf_on_doy(plist.isParameter("leaf on time [doy]") ? + plist.get("leaf on time [doy]") : + Relations::NaN), + leaf_off_doy(plist.isParameter("leaf off time [doy]") ? + plist.get("leaf off time [doy]") : + Relations::NaN), + pt_alpha_snow(plist.isParameter("Priestley-Taylor alpha of snow [-]") ? + readPositiveLandCoverParameter(plist, "Priestley-Taylor alpha of snow [-]") : + Relations::NaN), + pt_alpha_canopy(plist.isParameter("Priestley-Taylor alpha of canopy [-]") ? + readPositiveLandCoverParameter(plist, "Priestley-Taylor alpha of canopy [-]") : + Relations::NaN), + pt_alpha_ground(plist.isParameter("Priestley-Taylor alpha of bare ground [-]") ? + readPositiveLandCoverParameter(plist, "Priestley-Taylor alpha of bare ground [-]") : + Relations::NaN), + pt_alpha_transpiration(plist.isParameter("Priestley-Taylor alpha of transpiration [-]") ? + readPositiveLandCoverParameter(plist, "Priestley-Taylor alpha of transpiration [-]") : + Relations::NaN), + albedo_ground(plist.isParameter("albedo of bare ground [-]") ? + readZeroOneLandCoverParameter(plist, "albedo of bare ground [-]") : + Relations::NaN), + emissivity_ground(plist.isParameter("emissivity of bare ground [-]") ? + readZeroOneLandCoverParameter(plist, "emissivity of bare ground [-]") : + Relations::NaN), + albedo_canopy(plist.isParameter("albedo of canopy [-]") ? + readZeroOneLandCoverParameter(plist, "albedo of canopy [-]") : + Relations::NaN), + beers_k_sw(plist.isParameter("Beer's law extinction coefficient, shortwave [-]") ? + readPositiveLandCoverParameter(plist, "Beer's law extinction coefficient, shortwave [-]") : + Relations::NaN), + beers_k_lw(plist.isParameter("Beer's law extinction coefficient, longwave [-]") ? + readPositiveLandCoverParameter(plist, "Beer's law extinction coefficient, longwave [-]") : + Relations::NaN), + snow_transition_depth(plist.isParameter("snow transition depth [m]") ? + readPositiveLandCoverParameter(plist, "snow transition depth [m]") : + Relations::NaN), + water_transition_depth(plist.isParameter("water transition depth [m]") ? + readPositiveLandCoverParameter(plist, "water transition depth [m]") : + Relations::NaN), + roughness_ground(plist.isParameter("roughness length of bare ground [m]") ? + readPositiveLandCoverParameter(plist, "roughness length of bare ground [m]") : + Relations::NaN), + roughness_snow(plist.isParameter("roughness length of snow [m]") ? + readPositiveLandCoverParameter(plist, "roughness length of snow [m]") : + Relations::NaN) {}