diff --git a/docs/documentation/source/ats_demos b/docs/documentation/source/ats_demos index 8297c37a3..8b8e59a15 160000 --- a/docs/documentation/source/ats_demos +++ b/docs/documentation/source/ats_demos @@ -1 +1 @@ -Subproject commit 8297c37a3e39885ea2931c3791b86b3aa68fe06f +Subproject commit 8b8e59a15aa6558f3064ae6fd4814869a15bb705 diff --git a/src/executables/CMakeLists.txt b/src/executables/CMakeLists.txt index adf7e704f..dc56f496c 100644 --- a/src/executables/CMakeLists.txt +++ b/src/executables/CMakeLists.txt @@ -40,7 +40,7 @@ include_directories(${ATS_SOURCE_DIR}/src/pks) include_directories(${ATS_SOURCE_DIR}/src/pks/mpc) include_directories(${ATS_SOURCE_DIR}/src/pks/energy) include_directories(${ATS_SOURCE_DIR}/src/pks/flow) -include_directories(${ATS_SOURCE_DIR}/src/pks/deformation) +include_directories(${ATS_SOURCE_DIR}/src/pks/deform) include_directories(${ATS_SOURCE_DIR}/src/pks/transport) include_directories(${ATS_SOURCE_DIR}/src/operators/upwinding) include_directories(${ATS_SOURCE_DIR}/src/operators/advection) @@ -182,3 +182,12 @@ add_amanzi_executable(ats LINK_LIBS ats_executable ${fates_link_libs} ${tpl_link_libs} ${ats_link_libs} ${amanzi_link_libs} OUTPUT_NAME ats OUTPUT_DIRECTORY ${ATS_BINARY_DIR}) + + +# set this here for now +# add option to bootstrap before completion +set (BUILD_ELM_ATS_API true) +if (BUILD_ELM_ATS_API) + message("building elm_ats api") + add_subdirectory(elm_ats_api) +endif() diff --git a/src/executables/elm_ats_api/CMakeLists.txt b/src/executables/elm_ats_api/CMakeLists.txt new file mode 100644 index 000000000..37db82422 --- /dev/null +++ b/src/executables/elm_ats_api/CMakeLists.txt @@ -0,0 +1,47 @@ +# -*- mode: cmake -*- +# +# ELM_ATS_API +# library +# +# builds both dynamic and static libs for now +# need to add build logic to make configurable + +# verify the compatibility of the C/Fortran and C++/Fortran compilers +# Amanzi probably does this somewhere upstream, should check +include(FortranCInterface) +FortranCInterface_VERIFY(CXX) + +include_directories(${ATS_SOURCE_DIR}/src/executables) + +set(elm_ats_src_files + ${ATS_SOURCE_DIR}/src/executables/coordinator.cc + ${ATS_SOURCE_DIR}/src/executables/ats_mesh_factory.cc + elm_ats_driver.cc + elm_ats_api.cc + ) + +set(elm_ats_inc_files + ${ATS_SOURCE_DIR}/src/executables/coordinator.hh + ${ATS_SOURCE_DIR}/src/executables/ats_mesh_factory.hh + elm_ats_driver.hh + elm_ats_api.h + ) + +## build shared lib +add_amanzi_library(elm_ats + SOURCE ${elm_ats_src_files} + HEADERS ${elm_ats_inc_files} + LINK_LIBS ${fates_link_libs} ${tpl_link_libs} ${ats_link_libs} ${amanzi_link_libs}) +if (APPLE AND BUILD_SHARED_LIBS) + set_target_properties(elm_ats PROPERTIES LINK_FLAGS "-Wl,-undefined,dynamic_lookup") +endif() + +# ## build static lib +# add_amanzi_library(elm_ats_static +# STATIC +# SOURCE ${elm_ats_src_files} +# HEADERS ${elm_ats_inc_files} +# LINK_LIBS ${fates_link_libs} ${tpl_link_libs} ${ats_link_libs} ${amanzi_link_libs}) + +# re-enable later +##add_subdirectory(test) diff --git a/src/executables/elm_ats_api/Indexer.hh b/src/executables/elm_ats_api/Indexer.hh new file mode 100644 index 000000000..f7be7d43c --- /dev/null +++ b/src/executables/elm_ats_api/Indexer.hh @@ -0,0 +1,108 @@ +/* + 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 +*/ +//! Tools for indexing arbitrary arrays in a variety of ways. + +#pragma once + + +namespace ATS { +namespace Utils { + +enum class Indexer_kind { + SCALAR = 0, + ELM, + ATS +}; + +enum class Domain_kind { + SURF, + SUBSURF +}; + + +template struct Indexer; + +// +// Indexer for scalars. +// +template +struct Indexer { + template + typename std::enable_if::type + get(T const * const& val, const int col, const int cic) const { return *val_; } + + template + typename std::enable_if::type + get(T const * const& val, const int i) const { return *val_; } + +private: + T* val_; +}; + + +// +// Indexer for ELM +// +template +struct Indexer { + + Indexer(int ncells_per_col=-1) + : ncells_per_col_(ncells_per_col) {} + + // for use by subsurface + template + typename std::enable_if::type + get(T * const& val, const int col, const int cic) const { + return val[col * ncells_per_col_ + cic]; + } + + // for use by surface + template + typename std::enable_if::type + get(T * const& val, const int i) const { + return val[i]; + } + + private: + int ncells_per_col_; +}; + + +// +// Indexer for ATS +// +template +struct Indexer { + + Indexer(const Amanzi::AmanziMesh::Mesh& mesh) + : mesh_(mesh) {} + + // for use by subsurface + template + typename std::enable_if::type + get(T * const& val, const int col, const int cic) const { + return val[mesh_.cells_of_column(col)[cic]]; + } + + // for use by surface + template + typename std::enable_if::type + get(T * const& val, const int i) const { + return val[i]; + } + + private: + const Amanzi::AmanziMesh::Mesh& mesh_; +}; + + + + + +} // namespace Utils +} // namespace ATS diff --git a/src/executables/elm_ats_api/elm_ats_api.cc b/src/executables/elm_ats_api/elm_ats_api.cc new file mode 100644 index 000000000..c2bebe3d7 --- /dev/null +++ b/src/executables/elm_ats_api/elm_ats_api.cc @@ -0,0 +1,157 @@ +/* + 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: Joe Beisman + Fenming Yuan + Ethan Coon +*/ +//! Wrapper for driving ATS from ELM. + +#include "elm_ats_driver.hh" +#include "elm_ats_api.h" + +#ifdef __cplusplus +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) +{ + return reinterpret_cast(ATS::createELM_ATSDriver(f_comm, input_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; +} + +// 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); +} + +// 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_setup(ELM_ATSDriver_ptr ats) +{ + reinterpret_cast(ats)->setup(); +} + +// 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_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); +} + + +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); +} + + +void ats_set_soil_hydrologic_properties(ELM_ATSDriver_ptr ats, + double const * const effective_porosity) +{ + reinterpret_cast(ats) + ->set_soil_hydrologic_properties(effective_porosity); +} + + +// 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 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_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); +} + + +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 +} +#endif diff --git a/src/executables/elm_ats_api/elm_ats_api.h b/src/executables/elm_ats_api/elm_ats_api.h new file mode 100644 index 000000000..0b2a26545 --- /dev/null +++ b/src/executables/elm_ats_api/elm_ats_api.h @@ -0,0 +1,133 @@ +/* + 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: Joe Beisman + Fenming Yuan + Ethan Coon + +This file defines a C interface to the ATS library +for use with LSMs. +*/ +//! Wrapper for driving ATS from ELM. + +#ifndef ELM_ATS_API_HH_ +#define ELM_ATS_API_HH_ + +#ifdef __cplusplus +extern "C" { + +// opaque pointer +// 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; + +#else +// calling code should not dereference the pointer to the ATS object +// pointer hidden behind typedef to discourage +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); + +// 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 +// ----------------------------------------------------------------------------- + +// call driver get_mesh_info() +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); + +// 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); + +// +// 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); + +// 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); + + +// +// Called after timestep advance +// ----------------------------------------------------------------------------- +// 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); + +#ifdef __cplusplus +} +#endif + +// include guard +#endif diff --git a/src/executables/elm_ats_api/elm_ats_driver.cc b/src/executables/elm_ats_api/elm_ats_driver.cc new file mode 100644 index 000000000..a8782a48b --- /dev/null +++ b/src/executables/elm_ats_api/elm_ats_driver.cc @@ -0,0 +1,739 @@ +#include +#include +#include +#include "errors.hh" +#include "dbc.hh" + +#include "Teuchos_ParameterList.hpp" +#include "Teuchos_ParameterXMLFileReader.hpp" +#include "Teuchos_XMLParameterListHelpers.hpp" +#include "Teuchos_TimeMonitor.hpp" +#include "Teuchos_StandardParameterEntryValidators.hpp" +#include "Teuchos_VerboseObjectParameterListHelpers.hpp" +#include "VerboseObject.hh" + +// registration files +#include "ats_registration_files.hh" + +// include fenv if it exists +#include "boost/version.hpp" +#if (BOOST_VERSION / 100 % 1000 >= 46) +#include "boost/config.hpp" +#ifndef BOOST_NO_FENV_H +#ifdef _GNU_SOURCE +#define AMANZI_USE_FENV +#include "boost/detail/fenv.hpp" +#endif +#endif +#endif +#include "boost/filesystem.hpp" + +#include "AmanziComm.hh" +#include "CompositeVector.hh" +#include "IO.hh" +#include "UnstructuredObservations.hh" +#include "pk_helpers.hh" +#include "elm_ats_driver.hh" + +namespace ATS { + +ELM_ATSDriver* +createELM_ATSDriver(MPI_Fint *f_comm, const char *infile, int npfts) { + // -- create communicator & get process rank + //auto comm = getDefaultComm(); + auto c_comm = MPI_Comm_f2c(*f_comm); + auto comm = getComm(c_comm); + auto rank = comm->MyPID(); + + // convert input file to std::string for easier handling + // infile must be null-terminated + std::string input_filename(infile); + + // check validity of input file name + if (input_filename.empty()) { + if (rank == 0) + std::cerr << "ERROR: no input file provided" << std::endl; + } else if (!boost::filesystem::exists(input_filename)) { + if (rank == 0) + std::cerr << "ERROR: input file \"" << input_filename << "\" does not exist." << std::endl; + } + + // -- parse input file + Teuchos::RCP plist = Teuchos::getParametersFromXmlFile(input_filename); + return new ELM_ATSDriver(plist, comm, npfts); +} + + +ELM_ATSDriver::ELM_ATSDriver(const Teuchos::RCP& plist, + const Comm_ptr_type& comm, + int npfts) + : Coordinator(plist, comm), + npfts_(npfts), + ncolumns_(-1), + ncells_per_col_(-1) +{ + // -- set default verbosity level to no output + // -- TODO make the verbosity level an input argument + VerboseObject::global_default_level = Teuchos::VERB_NONE; + + // domain names + domain_subsurf_ = Keys::readDomain(*plist_, "domain", "domain"); + domain_surf_ = Keys::readDomainHint(*plist_, domain_subsurf_, "subsurface", "surface"); + + // 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"); + + // 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"); + + // 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"); + + // 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"); + + // 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 + + // -- build columns to allow indexing by column + mesh_subsurf_->build_columns(); + + // -- check that number of surface cells = number of columns + ncolumns_ = mesh_surf_->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED); + AMANZI_ASSERT(ncolumns_ == mesh_subsurf_->num_columns(false)); + + // -- 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 + auto& col_zero = mesh_subsurf_->cells_of_column(0); + ncells_per_col_ = col_zero.size(); + for (int col=0; col!=ncolumns_; ++col) + AMANZI_ASSERT(mesh_subsurf_->cells_of_column(col).size() == ncells_per_col_); +} + + +void +ELM_ATSDriver::setup() +{ + // potential fluxes (ELM -> ATS) + requireAtNext(pot_infilt_key_, Amanzi::Tags::NEXT, *S_, pot_infilt_key_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireAtNext(pot_evap_key_, Amanzi::Tags::NEXT, *S_, pot_evap_key_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireAtNext(pot_trans_key_, Amanzi::Tags::NEXT, *S_, pot_trans_key_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + + // subsurface properties + requireAtNext(base_poro_key_, Amanzi::Tags::NEXT, *S_, base_poro_key_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireAtNext(perm_key_, Amanzi::Tags::NEXT, *S_, perm_key_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + + // dynamic subsurface properties + requireAtNext(root_frac_key_, Amanzi::Tags::NEXT, *S_, root_frac_key_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireAtNext(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) + requireAtNext(ch_b_key_, Amanzi::Tags::NEXT, *S_, ch_b_key_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireAtNext(ch_smpsat_key_, Amanzi::Tags::NEXT, *S_, ch_smpsat_key_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireAtNext(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) + requireAtNext(pd_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireAtNext(wtd_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + + // per cell ATS water state + requireAtNext(pc_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireAtNext(sat_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + + // mesh data + // requireAtNext(lat_key_, Amanzi::Tags::NEXT, *S_) + // .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + // requireAtNext(lon_key_, Amanzi::Tags::NEXT, *S_) + // .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + // requireAtNext(elev_key_, Amanzi::Tags::NEXT, *S_) + // .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireAtNext(surf_cv_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireAtNext(cv_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + + // may not be necessary? any PK that utilizes this should already have density + requireAtNext(surf_mol_dens_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireAtNext(surf_mass_dens_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireAtNext(subsurf_mol_dens_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireAtNext(subsurf_mass_dens_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + + requireAtNext(total_trans_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireAtNext(wc_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + + + + requireAtNext(evap_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireAtNext(trans_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + + + Coordinator::setup(); + // NOTE: These must be called after Coordinator::setup() until PK_Phys_Default modifies the state eval + // list in constructor -- otherwise this primary variable is not available + // yet. See amanzi/ats#167 --ETC + requireAtNext(infilt_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireAtNext(pres_key_, Amanzi::Tags::NEXT, *S_, "flow") + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); +} + + + +// +// 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) +{ + ncols_local = ncolumns_; + ncols_global = mesh_surf_->cell_map(AmanziMesh::Entity_kind::CELL).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_->cells_of_column(0); + const auto& fc = mesh_subsurf_->face_centroid(mesh_subsurf_->faces_of_column(0)[0]); + for (int i=0; i!=ncells_per_col_; ++i) { + depth[i] = fc[2] - mesh_subsurf_->cell_centroid(cells_in_col[i])[2]; + } + + // hard-coded Toledo OH for now... + for (int i=0; i!=ncolumns_; ++i) { + lat[i] = 41.65; + lon[i] = -83.54; + } +} + + +void ELM_ATSDriver::initialize(double t, + double const * const elm_water_content, + double const * const elm_pressure) +{ + // Assign start time to ATS + t0_ = t; + + // initialize ATS data, commit initial conditions + Coordinator::initialize(); + + // initialize pressure field + ELM_ATSDriver::init_pressure_from_wc_(elm_water_content); + + // 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_); + + // 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_->cells_of_column(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); + } + } + } + + // 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) +{ + double dt_subcycle = dt; + double t_end = S_->get_time() + dt_subcycle; + + 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); + + // solve model for a single timestep + 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); + } + + dt_subcycle = Coordinator::get_dt(fail); + } // while + + 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); + } +} // advance() + + +// simulates external timeloop with dt coming from calling model +void ELM_ATSDriver::advance_test() +{ + while (S_->get_time() < t1_) { + // use dt from ATS for testing + double dt = Coordinator::get_dt(false); + // call main method + advance(dt, false, false); + } +} + +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) +{ + // 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; + } + + 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_->faces_of_column(i); + const auto& cells_of_col = mesh_subsurf_->cells_of_column(i); + for (int j=0; j!=ncells_per_col_; ++j) { + const double dz = mesh_subsurf_->face_centroid(faces[j])[2] - mesh_subsurf_->face_centroid(faces[j + 1])[2]; + sat_liq[j * ncolumns_ + i] = satl[0][cells_of_col[j]] * por[0][cells_of_col[j]] * dens[0][cells_of_col[j]] * dz; + } + } + + // 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_of_col = mesh_subsurf_->cells_of_column(i); +// for (int j=0; j!=ncells_per_col_; ++j) { +// matric_potential[j * ncolumns_ + i] = pc[0][cells_of_col[j]] * g_inv; +// soil_water_potential[j * ncolumns_ + i] = 0.101325 - 1.e-6 * pres[0][cells_of_col[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) +{ + // 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_->faces_of_column(i); + const auto& cells = mesh_subsurf_->cells_of_column(i); + for (int j=0; j!=ncells_per_col_; ++j) { + double dz = mesh_subsurf_->face_centroid(faces[j])[2] - mesh_subsurf_->face_centroid(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 +} + + +void ELM_ATSDriver::initZero_(const Key& key) +{ + auto& vec = S_->GetW(key, Amanzi::Tags::NEXT, key); + vec.PutScalar(0.); + S_->GetRecordW(key, Amanzi::Tags::NEXT, key).set_initialized(); +} + + +void ELM_ATSDriver::copyToSurf_(double const * const in, const Key& key, Key owner) +{ + 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_); +} + + +// +// 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) +{ + 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_of_col = mesh_subsurf_->cells_of_column(i); + for (int j=0; j!=ncells_per_col_; ++j) { + vec[0][cells_of_col[j]] = in[j * ncolumns_ + i]; + } + } + + 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]; +} + + +// +// 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 +{ + 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) { + const auto& cells_of_col = mesh_subsurf_->cells_of_column(i); + for (int j=0; j!=ncells_per_col_; ++j) { + out[j * ncolumns_ + i] = vec[0][cells_of_col[j]]; + } + } +} + + +} // namespace ATS diff --git a/src/executables/elm_ats_api/elm_ats_driver.hh b/src/executables/elm_ats_api/elm_ats_driver.hh new file mode 100644 index 000000000..9544b7089 --- /dev/null +++ b/src/executables/elm_ats_api/elm_ats_driver.hh @@ -0,0 +1,169 @@ +/* + 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: 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 + +#include "Teuchos_RCP.hpp" +#include "Epetra_MultiVector.h" + +#include "Mesh.hh" +#include "MeshPartition.hh" +#include "Key.hh" +#include "coordinator.hh" + +namespace ATS { + +using namespace Amanzi; + +class ELM_ATSDriver : public Coordinator { + +public: + + ELM_ATSDriver(const Teuchos::RCP& plist, + const Comm_ptr_type& comm, + 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); + + void setup(); + + 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); + + private: + void init_pressure_from_wc_(double const * const elm_water_content); + + 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); + + + private: + Teuchos::RCP elm_list_; + + int ncolumns_; + int ncells_per_col_; + int npfts_; + + Key domain_subsurf_; + Key domain_surf_; + + 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 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 sat_key_; + //Key sat_gas_key_; + //Key sat_ice_key_; + Key infilt_key_; + Key trans_key_; + Key evap_key_; + Key total_trans_key_; + Key surf_mol_dens_key_; + Key surf_mass_dens_key_; + Key subsurf_mol_dens_key_; + Key subsurf_mass_dens_key_; + Key surf_cv_key_; + Key cv_key_; + +}; + + +// +// 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); + +} // namespace ATS + + + + diff --git a/src/executables/elm_ats_api/test/CMakeLists.txt b/src/executables/elm_ats_api/test/CMakeLists.txt new file mode 100644 index 000000000..02b3b84de --- /dev/null +++ b/src/executables/elm_ats_api/test/CMakeLists.txt @@ -0,0 +1,18 @@ +# -*- mode: cmake -*- +# +# ELM_ATS_API +# example/test drivers +# +# build both tests for now +# need to make configurable + +include_directories(${ATS_SOURCE_DIR}/src/executables) +include_directories(${ATS_SOURCE_DIR}/src/executables/elm_ats_api) + +add_amanzi_executable(elm_ats_cc_test + SOURCE elm_ats_test.cc + LINK_LIBS elm_ats ${fates_link_libs} ${tpl_link_libs} ${ats_link_libs} ${amanzi_link_libs} + OUTPUT_NAME elm_ats_cc_test + OUTPUT_DIRECTORY ${ATS_BINARY_DIR}) + +add_subdirectory(fortran) diff --git a/src/executables/elm_ats_api/test/elm_ats_test.cc b/src/executables/elm_ats_api/test/elm_ats_test.cc new file mode 100644 index 000000000..438fc7523 --- /dev/null +++ b/src/executables/elm_ats_api/test/elm_ats_test.cc @@ -0,0 +1,140 @@ + +#include + +#include +#include +#include "Epetra_SerialComm.h" + +#include "Teuchos_ParameterList.hpp" +#include "Teuchos_ParameterXMLFileReader.hpp" +#include "Teuchos_XMLParameterListHelpers.hpp" +#include "Teuchos_CommandLineProcessor.hpp" +#include "Teuchos_StandardParameterEntryValidators.hpp" +#include "Teuchos_VerboseObjectParameterListHelpers.hpp" + +#include "VerboseObject_objs.hh" + +#include "ats_version.hh" +#include "amanzi_version.hh" +#include "tpl_versions.h" + +#include "dbc.hh" +#include "errors.hh" + +// registration files +#include "ats_registration_files.hh" + + +// include fenv if it exists +#include "boost/version.hpp" +#if (BOOST_VERSION / 100 % 1000 >= 46) +#include "boost/config.hpp" +#ifndef BOOST_NO_FENV_H +#ifdef _GNU_SOURCE +#define AMANZI_USE_FENV +#include "boost/detail/fenv.hpp" +#endif +#endif +#endif + +#include "boost/filesystem.hpp" + +#include "elm_ats_driver.hh" +#include "elm_ats_api.h" + +int main(int argc, char *argv[]) +{ + +#ifdef AMANZI_USE_FENV + // feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); + feraiseexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); +#endif + + Teuchos::GlobalMPISession mpiSession(&argc,&argv,0); + int rank = mpiSession.getRank(); + + std::string input_filename; + if ((argc >= 2) && (argv[argc-1][0] != '-')) { + input_filename = std::string(argv[argc-1]); + argc--; + } + + Teuchos::CommandLineProcessor clp; + clp.setDocString("Run ATS simulations for ecosystem hydrology.\n\nStandard usage: ats input.xml\n"); + + std::string opt_input_filename = ""; + clp.setOption("xml_file", &opt_input_filename, "XML input file"); + + bool version(false); + clp.setOption("version", "no_version", &version, "Print version number and exit."); + + bool print_version(false); + clp.setOption("print_version", "no_print_version", &print_version, "Print full version info and exit."); + + std::string verbosity; + clp.setOption("verbosity", &verbosity, "Default verbosity level: \"none\", \"low\", \"medium\", \"high\", \"extreme\"."); + + clp.throwExceptions(false); + clp.recogniseAllOptions(true); + + auto parseReturn = clp.parse(argc, argv); + if (parseReturn == Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED) { + return 0; + } + if (parseReturn != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) { + return 1; + } + + if (input_filename.empty() && !opt_input_filename.empty()) input_filename = opt_input_filename; + + // dummy data + // 1 col, 100 cells + double time = 0.; + int n = 1; + int m = 100; + int ncols_local, ncols_global, ncells_per_col; + std::vector soil_infil(n, 0.0); + std::vector soil_evap(n, 0.0); + std::vector air_pres(n); + std::vector surf_pres(n); + std::vector elev(n); + std::vector surf_area_m2(n); + std::vector lat(n); + std::vector lon(n); + + std::vector dz(m); + std::vector depth(m); + std::vector root_tran(m); + std::vector soil_pres(m); + std::vector soil_pot(m); + std::vector satl(m); + std::vector sati(m); + std::vector pft_i(m); + + // dummy fortran comm + MPI_Fint comm = 0; + + // test driver directly + auto driver = std::unique_ptr(ATS::createELM_ATSDriver(&comm, input_filename.data())); + driver->setup(); + //driver->get_mesh_info(ncols_local, ncols_global, lat.data(), lon.data(), elev.data(), surf_area_m2.data(), pft_i.data(), ncells_per_col, depth.data()); + driver->initialize(); + //driver->set_potential_sources(soil_infil.data(), soil_evap.data(), root_tran.data()); + driver->advance_test(); + //driver->get_waterstate(surf_pres.data(), soil_pres.data(), soil_pot.data(), satl.data(), sati.data()); + driver->finalize(); + + // test api + //auto driver_api = ats_create_c(&comm, input_filename.data()); + //ats_setup_c(driver_api); + //ats_get_mesh_info_c(driver_api, &ncols_local, &ncols_global, lat.data(), lon.data(), + // elev.data(), surf_area_m2.data(), pft_i.data(), &ncells_per_col, depth.data()); + //ats_initialize_c(driver_api, &time, air_pres.data(), soil_pres.data()); + //ats_set_potential_sources_c(driver_api, soil_infil.data(), soil_evap.data(), root_tran.data()); + //ats_advance_test_c(driver_api); + //ats_get_waterstate_c(driver_api, surf_pres.data(), soil_pres.data(), soil_pot.data(), satl.data(), sati.data()); + //ats_delete_c(driver_api); + std::cout << "DONE WITH ELM-ATS C++ TEST" << std::endl; + + return 0; +} diff --git a/src/executables/elm_ats_api/test/fortran/CMakeLists.txt b/src/executables/elm_ats_api/test/fortran/CMakeLists.txt new file mode 100644 index 000000000..1dedda87a --- /dev/null +++ b/src/executables/elm_ats_api/test/fortran/CMakeLists.txt @@ -0,0 +1,17 @@ +# -*- mode: cmake -*- +# +# ELM_ATS_API +# example/test Fortran driver +# + +set(fort_src_files + ats_mod.f90 + test.f90 + ) + +add_amanzi_executable(fort_driver + SOURCE ${fort_src_files} + LINK_LIBS elm_ats + OUTPUT_NAME elm_ats_f90_test + OUTPUT_DIRECTORY ${ATS_BINARY_DIR}) + diff --git a/src/executables/elm_ats_api/test/fortran/ats.f90 b/src/executables/elm_ats_api/test/fortran/ats.f90 new file mode 100644 index 000000000..85891b5d8 --- /dev/null +++ b/src/executables/elm_ats_api/test/fortran/ats.f90 @@ -0,0 +1,86 @@ + +!! Defines a Fortran interface to ATS library +!! Provides Fortran method declarations and bindings to ATS C API functions +!! The implementation of the methods declared here can be found in ats_mod.f90 +interface + + function ats_create_c(comm, infile) bind(c, name="ats_create") + use iso_c_binding + implicit none + type(c_ptr) :: ats_create_c + integer, intent(in) :: comm + character(len=1, kind=C_CHAR), intent(in) :: infile(*) + end function + + subroutine ats_delete_c(ats) bind(c, name="ats_delete") + use iso_c_binding + implicit none + type(c_ptr) :: ats + end subroutine + + subroutine ats_setup_c(ats) bind(c, name="ats_setup") + use iso_c_binding + implicit none + type(c_ptr), value :: ats + end subroutine + + subroutine ats_initialize_c(ats) bind(c, name="ats_initialize") + use iso_c_binding + implicit none + type(c_ptr), value :: ats + end subroutine + + subroutine ats_advance_test_c(ats) bind(c, name="ats_advance_test") + use iso_c_binding + implicit none + type(c_ptr), value :: ats + end subroutine + + subroutine ats_advance_c(ats, dt) bind(c, name="ats_advance") + use iso_c_binding + implicit none + type(c_ptr), value :: ats + real(c_double), intent(in) :: dt + end subroutine + + subroutine ats_set_sources_c(ats, soil_infil, soil_evap, root_trans, ncols, ncells) & + bind(c, name="ats_set_sources") + use iso_c_binding + implicit none + type(c_ptr), value :: ats + real(c_double), dimension(*), intent(in) :: soil_infil + real(c_double), dimension(*), intent(in) :: soil_evap + real(c_double), dimension(*), intent(in) :: root_trans + integer(c_int), intent(in) :: ncols + integer(c_int), intent(in) :: ncells + end subroutine + + subroutine ats_get_waterstate_c(ats, surf_pres, soil_pres, satur, ncols, ncells) & + bind(c, name="ats_get_waterstate") + use iso_c_binding + implicit none + type(c_ptr), value :: ats + real(c_double), dimension(*), intent(in) :: surf_pres + real(c_double), dimension(*), intent(in) :: soil_pres + real(c_double), dimension(*), intent(in) :: satur + integer(c_int), intent(in) :: ncols + integer(c_int), intent(in) :: ncells + end subroutine + + subroutine ats_get_mesh_info_c(ats, ncols_local, ncols_global, ncells_per_col, dz, & + depth, elev, surf_area_m2, lat, lon) bind(c, name="ats_get_mesh_info") + use iso_c_binding + implicit none + type(c_ptr), value :: ats + real(c_double), dimension(*), intent(in) :: dz + real(c_double), dimension(*), intent(in) :: depth + real(c_double), dimension(*), intent(in) :: elev + real(c_double), dimension(*), intent(in) :: surf_area_m2 + real(c_double), dimension(*), intent(in) :: lat + real(c_double), dimension(*), intent(in) :: lon + integer(c_int), intent(in) :: ncols_local + integer(c_int), intent(in) :: ncols_global + integer(c_int), intent(in) :: ncells_per_col + end subroutine + +end interface diff --git a/src/executables/elm_ats_api/test/fortran/ats_mod.f90 b/src/executables/elm_ats_api/test/fortran/ats_mod.f90 new file mode 100644 index 000000000..01abed891 --- /dev/null +++ b/src/executables/elm_ats_api/test/fortran/ats_mod.f90 @@ -0,0 +1,130 @@ + +!! Defines a Fortran type that can store a pointer to ATS +!! Provides implementations of methods declared in ats.f90 +module libats + use iso_c_binding + private + public :: ats + include "ats.f90" + + type ats + private + type(c_ptr) :: ptr ! pointer to ats class + contains + final :: ats_delete + procedure :: setup => ats_setup + procedure :: initialize => ats_initialize + procedure :: advance => ats_advance + procedure :: advance_test => ats_advance_test + procedure :: set_sources => ats_set_sources + procedure :: get_waterstate => ats_get_waterstate + procedure :: get_mesh_info => ats_get_mesh_info + end type + + ! constructor + interface ats + procedure ats_create + end interface + +!! Definitions of Fortran interface methods +!! These methods act like wrappers around their C bindings. +contains + + !! Allocate and construct a new ATS object + !! Returns a pointer to new ATS object + function ats_create(comm, infile) + implicit none + type(ats) :: ats_create + integer, intent(in) :: comm + character(len=*), intent(in) :: infile + character(len=1, kind=C_CHAR) :: c_str_infile(len_trim(infile) + 1) + integer :: n_char, i + + !! basic parsing + n_char = len_trim(infile) + do i = 1, n_char + c_str_infile(i) = infile(i:i) + end do + c_str_infile(n_char + 1) = C_NULL_CHAR + + ats_create%ptr = ats_create_c(comm, c_str_infile) + end function + + !! Destroy ATS object + subroutine ats_delete(this) + implicit none + type(ats) :: this + call ats_delete_c(this%ptr) + end subroutine + + !! Setup ATS for LSM simulation + subroutine ats_setup(this) + implicit none + class(ats) :: this + call ats_setup_c(this%ptr) + end subroutine + + !! Initialize ATS for this simulation + subroutine ats_initialize(this) + implicit none + class(ats) :: this + call ats_initialize_c(this%ptr) + end subroutine + + !! Advance ATS by dt + subroutine ats_advance(this, dt) + implicit none + class(ats) :: this + double precision, intent(in) :: dt + call ats_advance_c(this%ptr, dt) + end subroutine + + !! Temporary for testing - implements time loop and calls ats_advance() + subroutine ats_advance_test(this) + implicit none + class(ats) :: this + call ats_advance_test_c(this%ptr) + end subroutine + + !! WIP + !! Sets ATS coupling variables infiltration, evaporation, and transpiration + subroutine ats_set_sources(this, soil_infil, soil_evap, root_trans, ncols, ncells) + implicit none + class(ats) :: this + double precision, dimension(*), intent(in) :: soil_infil + double precision, dimension(*), intent(in) :: soil_evap + double precision, dimension(*), intent(in) :: root_trans + integer, intent(in) :: ncols + integer, intent(in) :: ncells + call ats_set_sources_c(this%ptr, soil_infil, soil_evap, root_trans, ncols, ncells) + end subroutine + + !! WIP + subroutine ats_get_waterstate(this, surf_pres, soil_pres, satur, ncols, ncells) + implicit none + class(ats) :: this + double precision, dimension(*), intent(in) :: surf_pres + double precision, dimension(*), intent(in) :: soil_pres + double precision, dimension(*), intent(in) :: satur + integer, intent(in) :: ncols + integer, intent(in) :: ncells + call ats_get_waterstate_c(this%ptr, surf_pres, soil_pres, satur, ncols, ncells) + end subroutine + + !! WIP + subroutine ats_get_mesh_info(this, ncols_local, ncols_global, ncells_per_col, dz, depth, elev, surf_area_m2, lat, lon) + implicit none + class(ats) :: this + double precision, dimension(*), intent(in) :: dz + double precision, dimension(*), intent(in) :: depth + double precision, dimension(*), intent(in) :: elev + double precision, dimension(*), intent(in) :: surf_area_m2 + double precision, dimension(*), intent(in) :: lat + double precision, dimension(*), intent(in) :: lon + integer, intent(in) :: ncols_local + integer, intent(in) :: ncols_global + integer, intent(in) :: ncells_per_col + call ats_get_mesh_info_c(this%ptr, ncols_local, ncols_global, ncells_per_col, dz, depth, elev, surf_area_m2, lat, lon) + end subroutine + +end module diff --git a/src/executables/elm_ats_api/test/fortran/test.f90 b/src/executables/elm_ats_api/test/fortran/test.f90 new file mode 100644 index 000000000..3bfcadff3 --- /dev/null +++ b/src/executables/elm_ats_api/test/fortran/test.f90 @@ -0,0 +1,65 @@ + +program elm_test + use libats + implicit none + include 'mpif.h' + type(ats) :: ats_driver + Character(len = 200) :: infile_name + integer :: ierror, i + + ! dummy data + ! 1 column, 100 cells + integer, parameter :: ncol = 1 + integer, parameter :: ncell = 100 + double precision, dimension(ncol) :: infil + double precision, dimension(ncol) :: evap + double precision, dimension(ncol) :: surf_pres + double precision, dimension(ncol) :: elev + double precision, dimension(ncol) :: surf_area_m2 + double precision, dimension(ncol) :: lat + double precision, dimension(ncol) :: lon + + double precision, dimension(ncell) :: dz + double precision, dimension(ncell) :: depth + double precision, dimension(ncell) :: soil_pres + double precision, dimension(ncell) :: satur + double precision, dimension(ncell) :: tran + + integer :: ncols_local, ncols_global, ncells_per_col + character (len = 40) :: test_prefix + + infil(:) = 10.0 + evap(:) = 3.0 + + do i=1,ncell + tran(i) = (1.0 - (1.0/ncell)*i) + end do + + call get_command_argument(1, infile_name) + + ! remove prefix if sent from test suite + if (infile_name(1:11) == '--xml_file=') then + infile_name = trim(adjustl(infile_name(12:))) + end if + + call MPI_INIT(ierror) + + ! create an ATS driver object + ats_driver = ats(MPI_COMM_WORLD, infile_name) + + ! call ATS methods + call ats_driver%setup() + !!call ats_driver%get_mesh_info(ncols_local, ncols_global, ncells_per_col, dz, depth, elev, surf_area_m2, lat, lon) + call ats_driver%initialize() + !!call ats_driver%set_sources(infil, evap, tran, ncol, ncell) + call ats_driver%advance_test() + !!call ats_driver%get_waterstate(surf_pres, soil_pres, satur, ncol, ncell) + + ! don't need to call ats_driver%delete now that final is used + + print*,"DONE WITH ELM-ATS FORTRAN TEST" + + call MPI_FINALIZE(ierror) + +end program elm_test +