From d86fd2de91a57c2a86233c0d9208b6859290ccc8 Mon Sep 17 00:00:00 2001 From: Ethan Coon Date: Fri, 10 Oct 2025 14:41:39 -0600 Subject: [PATCH 01/16] working on updates to elm api --- src/executables/elm_ats_api/elm_ats_api.cc | 214 +++---- src/executables/elm_ats_api/elm_ats_api.h | 104 +--- src/executables/elm_ats_api/elm_ats_driver.cc | 552 ++++-------------- src/executables/elm_ats_api/elm_ats_driver.hh | 109 ++-- 4 files changed, 289 insertions(+), 690 deletions(-) diff --git a/src/executables/elm_ats_api/elm_ats_api.cc b/src/executables/elm_ats_api/elm_ats_api.cc index b113e24f1..61bc56c30 100644 --- a/src/executables/elm_ats_api/elm_ats_api.cc +++ b/src/executables/elm_ats_api/elm_ats_api.cc @@ -21,149 +21,123 @@ 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)); - } - // 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; - } +// 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)); +} - // 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(); +// 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; } +} - 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 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 setup() - void ats_setup(ELM_ATSDriver_ptr ats) - { - reinterpret_cast(ats)->setup(); - } +// call driver advance_test() +void ats_advance_test(ELM_ATSDriver_ptr ats) +{ + reinterpret_cast(ats)->advance_test(); +} - // 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_info(ELM_ATSDriver_ptr ats, + int * const ncols_local, + int * const ncols_global, + int * const nlevgrnd) +{ + ATS::ELM_ATSDriver::MeshInfo info = reinterpret_cast(ats)->getMeshInfo(); + *ncols_local = info.ncols_local; + *ncols_global = info.ncols_global; + *nlevgrnd = info.nlevgrnd; +} - 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); - } +// 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_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 advance(dt) +void ats_advance(ELM_ATSDriver_ptr ats, + double dt, + bool checkpoint, + bool visualize) +{ + reinterpret_cast(ats)->advance(dt, checkpoint, visualize); +} +// call driver advance_test() +void ats_advance_test(ELM_ATSDriver_ptr ats) +{ + reinterpret_cast(ats)->advance(dt, checkpoint, visualize); +} - void ats_set_soil_hydrologic_properties(ELM_ATSDriver_ptr ats, - double const* const effective_porosity) - { - reinterpret_cast(ats)->set_soil_hydrologic_properties(effective_porosity); - } +// +// Memory movement/data passing +// ----------------------------------------------------------------------------- +double ats_get_scalar(ELM_ATSDriver_ptr ats, int scalar_id) +{ + return reinterpret_cast(ats)->getScalar(reinterpret_cast(scalar_id)); +} - // 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); - } +void ats_set_scalar(ELM_ATSDriver_ptr ats, int scalar_id, double in) +{ + reinterpret_cast(ats)->setScalar(reinterpret_cast(scalar_id), in); +} +void ats_get_field(ELM_ATSDriver_ptr ats, int var_id, double * const in) +{ + reinterpret_cast(ats)->getField(reinterpret_cast(var_id), in); +} - // 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); - } +double const * ats_get_field_ptr(ELM_ATSDriver_ptr ats, int var_id) +{ + return reinterpret_cast(ats)->getFieldPtr(reinterpret_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(reinterpret_cast(var_id)); +} +void ats_set_field(ELM_ATSDriver_ptr ats, int var_id, double const * const in) +{ + reinterpret_cast(ats)->SetField(reinterpret_cast(var_id), in); +} - 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..2bb2a85fb 100644 --- a/src/executables/elm_ats_api/elm_ats_api.h +++ b/src/executables/elm_ats_api/elm_ats_api.h @@ -37,93 +37,51 @@ 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 +// Control // ----------------------------------------------------------------------------- -// call driver get_mesh_info() +// +// These quantities should be compared against ELM to ensure consistent setup +// 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, + int * const ncols_global, + int * const nlevgrnd, + double * const depth); -// 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 scalar_id); +void ats_set_scalar(ELM_ATSDriver_ptr ats, int scalar_id, double in); + +void ats_get_field(ELM_ATSDriver_ptr ats, int var_id, double * const in); +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..3ec5ba067 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.cc +++ b/src/executables/elm_ats_api/elm_ats_driver.cc @@ -21,14 +21,13 @@ #include "CompositeVector.hh" #include "IO.hh" #include "UnstructuredObservations.hh" -#include "PK_Helpers.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) -{ +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); @@ -42,7 +41,8 @@ createELM_ATSDriver(MPI_Fint* f_comm, const char* infile, int npfts) // 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; @@ -55,6 +55,7 @@ createELM_ATSDriver(MPI_Fint* f_comm, const char* infile, int npfts) return new ELM_ATSDriver(plist, wallclock_timer, teuchos_comm, comm, npfts); } + ELM_ATSDriver::ELM_ATSDriver(const Teuchos::RCP& plist, const Teuchos::RCP& wallclock_timer, const Teuchos::RCP>& teuchos_comm, @@ -73,7 +74,7 @@ ELM_ATSDriver::ELM_ATSDriver(const Teuchos::RCP& plist, domain_subsurf_ = Keys::readDomain(*plist_, "domain", "domain"); domain_surf_ = Keys::readDomainHint(*plist_, domain_subsurf_, "subsurface", "surface"); - // meshes + // meshes mesh_subsurf_ = S_->GetMesh(domain_subsurf_); mesh_surf_ = S_->GetMesh(domain_surf_); @@ -87,54 +88,36 @@ ELM_ATSDriver::ELM_ATSDriver(const Teuchos::RCP& plist, 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"); + 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"); + 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"); + 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"); + 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"); + 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"); + total_trans_key_ = Keys::readKey(*plist_, domain_surf_, "total transpiration", "total_transpiration"); // cell vol keys surf_cv_key_ = Keys::getKey(domain_surf_, "cell_volume"); @@ -155,7 +138,7 @@ ELM_ATSDriver::ELM_ATSDriver(const Teuchos::RCP& plist, // subsurface auto col_zero = mesh_subsurf_->columns.getCells(0); ncells_per_col_ = col_zero.size(); - for (int col = 0; col != ncolumns_; ++col) + for (int col=0; col!=ncolumns_; ++col) AMANZI_ASSERT(mesh_subsurf_->columns.getCells(col).size() == ncells_per_col_); } @@ -164,166 +147,101 @@ 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); + 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 - 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); + 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 - 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); + 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) - 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); + 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) - 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); + 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 - 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); + 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 - // requireEvaluatorAtNext(lat_key_, Amanzi::Tags::NEXT, *S_) + // requireAtNext(lat_key_, Amanzi::Tags::NEXT, *S_) // .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - // requireEvaluatorAtNext(lon_key_, Amanzi::Tags::NEXT, *S_) + // requireAtNext(lon_key_, Amanzi::Tags::NEXT, *S_) // .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - // requireEvaluatorAtNext(elev_key_, Amanzi::Tags::NEXT, *S_) + // requireAtNext(elev_key_, Amanzi::Tags::NEXT, *S_) // .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - requireEvaluatorAtNext(surf_cv_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_surf_) - ->AddComponent("cell", AmanziMesh::CELL, 1); - requireEvaluatorAtNext(cv_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_subsurf_) - ->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 - 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); - - 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); + 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); + + 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); Coordinator::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) +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]; - 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.; - } - - // hard-coded Toledo OH for now... - for (int i = 0; i != ncolumns_; ++i) { - lat[i] = 41.65; - lon[i] = -83.54; - } + MeshInfo info; + info.ncols_local = ncolumns_; + info.ncols_global = mesh_surf_->getMap(AmanziMesh::Entity_kind::CELL, false).NumGlobalElements(); + info.nlevgrnd = ncells_per_col_; + return info; } -void -ELM_ATSDriver::initialize(double t, - double const* const elm_water_content, - double const* const elm_pressure) +void ELM_ATSDriver::initialize(double t, double const * const elm_water_content) { // Assign start time to ATS t0_ = t; @@ -350,8 +268,7 @@ ELM_ATSDriver::initialize(double t, } // use incoming water content to initialize pressure field -void -ELM_ATSDriver::init_pressure_from_wc_(double const* const elm_water_content) +void ELM_ATSDriver::init_pressure_from_wc_(double const * const elm_water_content) { // gravity, atmospheric pressure, and liquid water density // hardwired for now @@ -362,20 +279,20 @@ ELM_ATSDriver::init_pressure_from_wc_(double const* const elm_water_content) // 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); + .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); + 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); + 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); + 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); + 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); @@ -388,11 +305,11 @@ ELM_ATSDriver::init_pressure_from_wc_(double const* const elm_water_content) // 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) { + 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) { + 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] @@ -408,7 +325,7 @@ ELM_ATSDriver::init_pressure_from_wc_(double const* const elm_water_content) sat_depth = 0.0; } sat_depth += dz; - pres[0][cells_of_col[j]] = p_atm + rho * g * (sat_depth - dz / 2); + pres[0][cells_of_col[j]] = p_atm + rho * g * (sat_depth - dz/2); } } } @@ -424,13 +341,12 @@ ELM_ATSDriver::init_pressure_from_wc_(double const* const elm_water_content) } -void -ELM_ATSDriver::advance(double dt, bool do_chkp, bool do_vis) +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 }; + 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); @@ -478,8 +394,7 @@ ELM_ATSDriver::advance(double dt, bool do_chkp, bool do_vis) // simulates external timeloop with dt coming from calling model -void -ELM_ATSDriver::advance_test() +void ELM_ATSDriver::advance_test() { while (S_->get_time() < t1_) { // use dt from ATS for testing @@ -489,8 +404,7 @@ ELM_ATSDriver::advance_test() } } -void -ELM_ATSDriver::finalize() +void ELM_ATSDriver::finalize() { WriteStateStatistics(*S_, *vo_); report_memory(); @@ -499,221 +413,7 @@ ELM_ATSDriver::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_->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; - } - } - - // 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) -{ - // 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 -} - - -void -ELM_ATSDriver::initZero_(const Key& key) +void ELM_ATSDriver::initZero_(const Key& key) { auto& vec = S_->GetW(key, Amanzi::Tags::NEXT, key); vec.PutScalar(0.); @@ -721,16 +421,17 @@ ELM_ATSDriver::initZero_(const Key& key) } -void -ELM_ATSDriver::copyToSurf_(double const* const in, const Key& key, Key owner) + + +void ELM_ATSDriver::copyToSurf_(double const * const in, const Key& key, Key owner) { - if (owner.empty() ) owner = key; + if (owner.empty()) owner = key; // surf maps directly into columns - auto& vec = - *S_->GetW(key, Amanzi::Tags::NEXT, owner).ViewComponent("cell", false); + 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]; + for (int i=0; i!=ncolumns_; ++i) vec[0][i] = in[i]; changedEvaluatorPrimary(key, Amanzi::Tags::NEXT, *S_); } @@ -740,16 +441,15 @@ ELM_ATSDriver::copyToSurf_(double const* const in, const Key& key, Key owner) // 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) +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); + if (owner.empty()) owner = key; + auto& vec = *S_->GetW(key, Amanzi::Tags::NEXT, owner) + .ViewComponent("cell", false); - for (int i = 0; i != ncolumns_; ++i) { + for (int i=0; i!=ncolumns_; ++i) { const auto& cells = mesh_subsurf_->columns.getCells(i); - for (int j = 0; j != ncells_per_col_; ++j) { + for (int j=0; j!=ncells_per_col_; ++j) { vec[0][cells[j]] = in[j * ncolumns_ + i]; } } @@ -758,13 +458,13 @@ ELM_ATSDriver::copyToSub_(double const* const in, const Key& key, Key owner) } -void -ELM_ATSDriver::copyFromSurf_(double* const out, const Key& key) const +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); + const auto& vec = *S_->Get(key, Amanzi::Tags::NEXT) + .ViewComponent("cell", false); - for (int i = 0; i != ncolumns_; ++i) out[i] = vec[0][i]; + for (int i=0; i!=ncolumns_; ++i) out[i] = vec[0][i]; } @@ -772,15 +472,15 @@ ELM_ATSDriver::copyFromSurf_(double* const out, const Key& key) const // 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 +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); + const auto& vec = *S_->Get(key, Amanzi::Tags::NEXT) + .ViewComponent("cell", false); - for (int i = 0; i != ncolumns_; ++i) { + for (int i=0; i!=ncolumns_; ++i) { const auto cells = mesh_subsurf_->columns.getCells(i); - for (int j = 0; j != ncells_per_col_; ++j) { + for (int j=0; j!=ncells_per_col_; ++j) { out[j * ncolumns_ + i] = vec[0][cells[j]]; } } diff --git a/src/executables/elm_ats_api/elm_ats_driver.hh b/src/executables/elm_ats_api/elm_ats_driver.hh index 910f45ef2..eb9627bbb 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 @@ -45,66 +21,51 @@ namespace ATS { using namespace Amanzi; class ELM_ATSDriver : public Coordinator { + public: + + struct MeshInfo { + int ncols_local; + int ncols_global; + int nlevgrnd; + }; + ELM_ATSDriver(const Teuchos::RCP& plist, const Teuchos::RCP& wallclock_timer, const Teuchos::RCP>& teuchos_comm, const Amanzi::Comm_ptr_type& comm, int npfts = 17); - void finalize(); + MeshInfo getMeshInfo(); + void setup(); + void initialize(); 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 advanceTest(); + void finalize(); - void setup(); + void setScalar(const ScalarID& scalar_id, double in); + double getScalar(const ScalarID& scalar_id); - 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 setField(const VarID& var_id, double const * const in); + void getField(const VarID& var_id, double * const out); + double const * getFieldPtr(const VarID& var_id); + double * getFieldPtrW(const VarID& var_id); private: - void init_pressure_from_wc_(double const* const elm_water_content); + template + void setField_(double const * const in); - 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); + template + void getField_(double * const out); + + void copyToSurf_(double const * const in, const Key& key); + void copyToSub_(double const * const in, const Key& key); + void copyFromSurf_(double * const out, const Key& key) const; + void copyFromSub_(double * const out, const Key& key) const; + + void initPressureFromWC_(double const * const elm_water_content); + void initZero_(const Key& key); private: Teuchos::RCP elm_list_; @@ -150,12 +111,18 @@ class ELM_ATSDriver : public Coordinator { 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); +ELM_ATSDriver* +createELM_ATSDriver(MPI_Fint *f_comm, const char *infile, int npfts=17); } // namespace ATS + + + + From 2c8b5316591f49b3ee8d0b8e36aedacd3412f953 Mon Sep 17 00:00:00 2001 From: Ethan Coon Date: Fri, 10 Oct 2025 16:45:05 -0600 Subject: [PATCH 02/16] rethink of elm api --- src/executables/elm_ats_api/ats_variables.hh | 27 ++ src/executables/elm_ats_api/elm_ats_api.cc | 18 +- src/executables/elm_ats_api/elm_ats_api.h | 6 +- src/executables/elm_ats_api/elm_ats_driver.cc | 330 ++++++------------ src/executables/elm_ats_api/elm_ats_driver.hh | 31 +- 5 files changed, 157 insertions(+), 255 deletions(-) create mode 100644 src/executables/elm_ats_api/ats_variables.hh 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..9c2d298d0 --- /dev/null +++ b/src/executables/elm_ats_api/ats_variables.hh @@ -0,0 +1,27 @@ +namespace ATS { +namespace ELM { + +enum class VarID : int { + BASE_POROSITY, + HYDRAULIC_CONDUCTIVITY, + CLAPP_HORBERGER_B, + CLAPP_HORBERGER_PSI_SAT, + RESIDUAL_SATURATION, + EFFECTIVE_POROSITY, + ROOT_FRACTION, + SURFACE_WATER_CONTENT, + WATER_CONTENT, + PRESSURE, + SURFACE_WATER_SOURCE, + POTENTIAL_EVAPORATION, + POTENTIAL_TRANSPIRATION, + EVAPORATION, + TRANSPIRATION, + RUNOFF, + BASEFLOW, + ELEVATION, +}; + + +} // 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 61bc56c30..4eba893df 100644 --- a/src/executables/elm_ats_api/elm_ats_api.cc +++ b/src/executables/elm_ats_api/elm_ats_api.cc @@ -107,15 +107,15 @@ void ats_advance_test(ELM_ATSDriver_ptr ats) // // Memory movement/data passing // ----------------------------------------------------------------------------- -double ats_get_scalar(ELM_ATSDriver_ptr ats, int scalar_id) -{ - return reinterpret_cast(ats)->getScalar(reinterpret_cast(scalar_id)); -} - -void ats_set_scalar(ELM_ATSDriver_ptr ats, int scalar_id, double in) -{ - reinterpret_cast(ats)->setScalar(reinterpret_cast(scalar_id), in); -} +// double ats_get_scalar(ELM_ATSDriver_ptr ats, int scalar_id) +// { +// return reinterpret_cast(ats)->getScalar(reinterpret_cast(scalar_id)); +// } + +// void ats_set_scalar(ELM_ATSDriver_ptr ats, int scalar_id, double in) +// { +// reinterpret_cast(ats)->setScalar(reinterpret_cast(scalar_id), in); +// } void ats_get_field(ELM_ATSDriver_ptr ats, int var_id, double * const in) { diff --git a/src/executables/elm_ats_api/elm_ats_api.h b/src/executables/elm_ats_api/elm_ats_api.h index 2bb2a85fb..4de8a6749 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; @@ -75,8 +75,8 @@ void ats_advance_test(ELM_ATSDriver_ptr ats); // // Memory movement/data passing // ----------------------------------------------------------------------------- -double ats_get_scalar(ELM_ATSDriver_ptr ats, int scalar_id); -void ats_set_scalar(ELM_ATSDriver_ptr ats, int scalar_id, double in); +// double ats_get_scalar(ELM_ATSDriver_ptr ats, int scalar_id); +// void ats_set_scalar(ELM_ATSDriver_ptr ats, int scalar_id, double in); void ats_get_field(ELM_ATSDriver_ptr ats, int var_id, double * const in); double const * ats_get_field_ptr(ELM_ATSDriver_ptr ats, int var_id); diff --git a/src/executables/elm_ats_api/elm_ats_driver.cc b/src/executables/elm_ats_api/elm_ats_driver.cc index 3ec5ba067..d20db03d3 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.cc +++ b/src/executables/elm_ats_api/elm_ats_driver.cc @@ -78,37 +78,29 @@ ELM_ATSDriver::ELM_ATSDriver(const Teuchos::RCP& plist, 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) + gross_water_source_key_ = Keys::readKey(*plist_, domain_surf_, "gross water source", "gross_water_source"); + key_map_[ELM::SURFACE_WATER_SOURCE] = gross_water_source_key_; pot_evap_key_ = Keys::readKey(*plist_, domain_surf_, "potential evaporation mps", "potential_evaporation_mps"); + key_map_[ELM::POTENTIAL_EVAPORATION] = pot_evap_key_; pot_trans_key_ = Keys::readKey(*plist_, domain_surf_, "potential transpiration mps", "potential_transpiration_mps"); + key_map_[ELM::POTENTIAL_TRANSPIRATION] = pot_trans_key_; // 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"); + surf_wc_key_ = Keys::readKey(*plist_, domain_surf_, "water content", "water_content"); + key_map_[ELM::SURFACE_WATER_CONTENT] = surface_wc_key_; + wc_key_ = Keys::readKey(*plist_, domain_subsurf_, "water content", "water_content"); + key_map_[ELM::WATER_CONTENT] = wc_key_; + + // actual water fluxes evap_key_ = Keys::readKey(*plist_, domain_surf_, "evaporation", "evaporation"); - trans_key_ = Keys::readKey(*plist_, domain_subsurf_, "transpiration", "transpiration"); + key_map_[ELM::EVAPORATION] = evap_key_; + col_trans_key_ = Keys::readKey(*plist_, domain_surf_, "total transpiration", "total_transpiration"); + key_map_[ELM::TRANSPIRATION] = col_trans_key_; + col_baseflow_key_ = Keys::readKey(*plist_, domain_surf_, "column total baseflow", "baseflow"); + key_map_[ELM::BASEFLOW] = col_baseflow_key_; + col_runoff_key_ = Keys::readKey(*plist_, domain_surf_, "runoff generation", "runoff_generation"); + key_map_[ELM::RUNOFF] = col_runoff_key_; // 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"); @@ -116,30 +108,14 @@ ELM_ATSDriver::ELM_ATSDriver(const Teuchos::RCP& plist, 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 - - // -- check that number of surface cells = number of columns + // -- mesh info 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 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_); } @@ -147,58 +123,18 @@ void ELM_ATSDriver::setup() { // potential fluxes (ELM -> ATS) - requireAtNext(pot_infilt_key_, Amanzi::Tags::NEXT, *S_, pot_infilt_key_) + requireAtNext(gross_water_source_key_, Amanzi::Tags::NEXT, *S_, gross_water_source_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_) @@ -208,24 +144,22 @@ ELM_ATSDriver::setup() 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); - - requireAtNext(infilt_key_, Amanzi::Tags::NEXT, *S_) + requireAtNext(col_trans_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); + requireAtNext(col_baseflow_key_, Amanzi::Tags::NEXT, *S_, col_baseflow_key_) + .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); + requireAtNext(col_runoff_key_, Amanzi::Tags::NEXT, *S_, col_runoff_key_) + .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); Coordinator::setup(); -} + requireAtNext(pres_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireAtNext(surf_pres_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); +} // @@ -241,7 +175,9 @@ MeshInfo ELM_ATSDriver::getMeshInfo() } -void ELM_ATSDriver::initialize(double t, double const * const elm_water_content) +void ELM_ATSDriver::initialize(double t, + double const * const elm_water_content, + double const * const elm_surf_water_content) { // Assign start time to ATS t0_ = t; @@ -249,148 +185,41 @@ void ELM_ATSDriver::initialize(double t, double const * const elm_water_content) // 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_(gross_water_source_key_); initZero_(pot_trans_key_); initZero_(pot_evap_key_); - initZero_(infilt_key_); - initZero_(trans_key_); - initZero_(evap_key_); - initZero_(total_trans_key_); + + // no current evaluators for these, treat as primary + initZero_(baseflow_key_); + initZero_(runoff_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_->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); - } - } - } - - // 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 + S_->Assign("dt", Amanzi::Tags::DEFAULT, "dt", dt); + S_->advance_time(Amanzi::Tags::NEXT, dt); + // solve model for a single timestep + bool fail = Coordinator::advance(); 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."); + Errors::Message msg("ELM_ATSDriver: advance(dt) failed. Make ATS subcycle for proper ELM use."); Exceptions::amanzi_throw(msg); } -} // advance() + + // 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); +} // simulates external timeloop with dt coming from calling model @@ -487,4 +316,63 @@ void ELM_ATSDriver::copyFromSub_(double * const out, const Key& key) const } +void ELM_ATSDriver::getField(const VarID& var_id, double * const out) +{ + switch(var_id) { + case ELM::VarID::SURFACE_WATER_CONTENT: getField_(out); break; + case ELM::VarID::WATER_CONTENT: getField_(out); break; + case ELM::VarID::EVAPORATION: getField_(out); break; + case ELM::VarID::TRANSPIRATION: getField_(out); break; + case ELM::VarID::RUNOFF: getField_(out); break; + case ELM::VarID::BASEFLOW: getField_(out); break; + default: throw std::runtime_error("Ungettable variable"); +} + + + +void ELM_ATSDriver::setField(const VarID& var_id, double const * const out) +{ + switch(var_id) { + case ELM::VarID::BASE_POROSITY: setField_(out); break; + case ELM::VarID::HYDRAULIC_CONDUCTIVITY: setField_(out); break; + case ELM::VarID::CLAPP_HORNBERGER_B: setField_(out); break; + case ELM::VarID::CLAPP_HORNBERGER_PSI_SAT: setField_(out); break; + case ELM::VarID::RESIDUAL_SATURATION: setField_(out); break; + case ELM::VarID::EFFECTIVE_POROSITY: setField_(out); break; + case ELM::VarID::ROOT_FRACTION: setField_(out); break; + case ELM::VarID::SURFACE_WATER_CONTENT: setField_(out); break; + case ELM::VarID::WATER_CONTENT: setField_(out); break; + case ELM::VarID::SURFACE_WATER_SOURCE: setField_(out); break; + case ELM::VarID::POTENTIAL_EVAPORATION: setField_(out); break; + case ELM::VarID::POTENTIAL_TRANSPIRATION: setField_(out); break; + default: throw std::runtime_error("Unsettable variable"); +} + + +template +void ELM_ATSDriver::getField_(double * const out) +{ + Key field_key = keys_[var_id]; + if (Keys::getDomain(field_key) == domain_subsurf_) { + copyFromSub_(out, field_key); + } else { + copyFromSurf_(out, field_key); + } +} + + +template +void ELM_ATSDriver::setField_(double const * const out) +{ + Key field_key = keys_[var_id]; + if (Keys::getDomain(field_key) == domain_subsurf_) { + copyToSub_(out, field_key); + } else { + copyToSurf_(out, field_key); + } +} + + + + } // namespace ATS diff --git a/src/executables/elm_ats_api/elm_ats_driver.hh b/src/executables/elm_ats_api/elm_ats_driver.hh index eb9627bbb..f6d14bb0a 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.hh +++ b/src/executables/elm_ats_api/elm_ats_driver.hh @@ -64,7 +64,6 @@ class ELM_ATSDriver : public Coordinator { void copyFromSurf_(double * const out, const Key& key) const; void copyFromSub_(double * const out, const Key& key) const; - void initPressureFromWC_(double const * const elm_water_content); void initZero_(const Key& key); private: @@ -80,35 +79,23 @@ 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 sat_key_; - //Key sat_gas_key_; - //Key sat_ice_key_; - Key infilt_key_; - Key trans_key_; + Key surf_wc_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 surf_cv_key_; Key cv_key_; From 7a299451cf8c26c106ef216efba9a9f3ba7b852b Mon Sep 17 00:00:00 2001 From: Ethan Coon Date: Mon, 13 Oct 2025 13:24:57 -0600 Subject: [PATCH 03/16] wip --- src/executables/ats_variables.inc | 31 ++++++++++++++++ src/executables/elm_ats_api/ats_variables.hh | 39 ++++++++++---------- src/executables/elm_ats_api/elm_ats_api.h | 4 +- 3 files changed, 53 insertions(+), 21 deletions(-) create mode 100644 src/executables/ats_variables.inc diff --git a/src/executables/ats_variables.inc b/src/executables/ats_variables.inc new file mode 100644 index 000000000..ac99d294e --- /dev/null +++ b/src/executables/ats_variables.inc @@ -0,0 +1,31 @@ +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 + 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 index 9c2d298d0..c2daed790 100644 --- a/src/executables/elm_ats_api/ats_variables.hh +++ b/src/executables/elm_ats_api/ats_variables.hh @@ -2,26 +2,27 @@ namespace ATS { namespace ELM { enum class VarID : int { - BASE_POROSITY, - HYDRAULIC_CONDUCTIVITY, - CLAPP_HORBERGER_B, - CLAPP_HORBERGER_PSI_SAT, - RESIDUAL_SATURATION, - EFFECTIVE_POROSITY, - ROOT_FRACTION, - SURFACE_WATER_CONTENT, - WATER_CONTENT, - PRESSURE, - SURFACE_WATER_SOURCE, - POTENTIAL_EVAPORATION, - POTENTIAL_TRANSPIRATION, - EVAPORATION, - TRANSPIRATION, - RUNOFF, - BASEFLOW, - ELEVATION, + NULL = 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, + SURFACE_WATER_SOURCE = 11, + POTENTIAL_EVAPORATION = 12, + POTENTIAL_TRANSPIRATION = 13, + EVAPORATION = 14, + TRANSPIRATION = 15, + RUNOFF = 16, + BASEFLOW = 17, + ELEVATION = 18, + TIME = 19 }; - } // namespace ELM } // namespace ATS diff --git a/src/executables/elm_ats_api/elm_ats_api.h b/src/executables/elm_ats_api/elm_ats_api.h index 4de8a6749..03db2e838 100644 --- a/src/executables/elm_ats_api/elm_ats_api.h +++ b/src/executables/elm_ats_api/elm_ats_api.h @@ -49,7 +49,7 @@ void ats_get_mesh_info(ELM_ATSDriver_ptr ats, int * const ncols_local, int * const ncols_global, int * const nlevgrnd, - double * const depth); + double * const dzs); // // simulation setup @@ -78,7 +78,7 @@ void ats_advance_test(ELM_ATSDriver_ptr ats); // double ats_get_scalar(ELM_ATSDriver_ptr ats, int scalar_id); // void ats_set_scalar(ELM_ATSDriver_ptr ats, int scalar_id, double in); -void ats_get_field(ELM_ATSDriver_ptr ats, int var_id, double * const 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); From 0a336178f3b2abfefc2958df822bd8321bcf22e9 Mon Sep 17 00:00:00 2001 From: Ethan Coon Date: Mon, 13 Oct 2025 21:26:14 -0600 Subject: [PATCH 04/16] updating more generic api --- .../ats_variables.F90} | 0 src/executables/elm_ats_api/ats_variables.hh | 4 +- src/executables/elm_ats_api/elm_ats_api.cc | 48 ++++------ src/executables/elm_ats_api/elm_ats_driver.cc | 87 +++++++++---------- src/executables/elm_ats_api/elm_ats_driver.hh | 21 +++-- 5 files changed, 76 insertions(+), 84 deletions(-) rename src/executables/{ats_variables.inc => elm_ats_api/ats_variables.F90} (100%) diff --git a/src/executables/ats_variables.inc b/src/executables/elm_ats_api/ats_variables.F90 similarity index 100% rename from src/executables/ats_variables.inc rename to src/executables/elm_ats_api/ats_variables.F90 diff --git a/src/executables/elm_ats_api/ats_variables.hh b/src/executables/elm_ats_api/ats_variables.hh index c2daed790..d8a4c1155 100644 --- a/src/executables/elm_ats_api/ats_variables.hh +++ b/src/executables/elm_ats_api/ats_variables.hh @@ -2,7 +2,7 @@ namespace ATS { namespace ELM { enum class VarID : int { - NULL = 0, + NO_VARIABLE = 0, BASE_POROSITY = 1, HYDRAULIC_CONDUCTIVITY = 2, CLAPP_HORNBERGER_B = 3, @@ -13,7 +13,7 @@ enum class VarID : int { SURFACE_WATER_CONTENT = 8, WATER_CONTENT = 9, PRESSURE = 10, - SURFACE_WATER_SOURCE = 11, + GROSS_SURFACE_WATER_SOURCE = 11, POTENTIAL_EVAPORATION = 12, POTENTIAL_TRANSPIRATION = 13, EVAPORATION = 14, diff --git a/src/executables/elm_ats_api/elm_ats_api.cc b/src/executables/elm_ats_api/elm_ats_api.cc index 4eba893df..c3f18c6fd 100644 --- a/src/executables/elm_ats_api/elm_ats_api.cc +++ b/src/executables/elm_ats_api/elm_ats_api.cc @@ -22,6 +22,8 @@ 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) { @@ -49,25 +51,11 @@ void ats_delete(ELM_ATSDriver_ptr ats) } -// 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, - int * const nlevgrnd) + int * const nlevgrnd, + double * const dzs) { ATS::ELM_ATSDriver::MeshInfo info = reinterpret_cast(ats)->getMeshInfo(); *ncols_local = info.ncols_local; @@ -97,45 +85,45 @@ void ats_advance(ELM_ATSDriver_ptr ats, reinterpret_cast(ats)->advance(dt, checkpoint, visualize); } + // call driver advance_test() void ats_advance_test(ELM_ATSDriver_ptr ats) { - reinterpret_cast(ats)->advance(dt, checkpoint, visualize); + reinterpret_cast(ats)->advanceTest(); } - // // Memory movement/data passing // ----------------------------------------------------------------------------- -// double ats_get_scalar(ELM_ATSDriver_ptr ats, int scalar_id) -// { -// return reinterpret_cast(ats)->getScalar(reinterpret_cast(scalar_id)); -// } +double ats_get_scalar(ELM_ATSDriver_ptr ats, int var_id) +{ + return reinterpret_cast(ats)->getScalar(static_cast(var_id)); +} -// void ats_set_scalar(ELM_ATSDriver_ptr ats, int scalar_id, double in) -// { -// reinterpret_cast(ats)->setScalar(reinterpret_cast(scalar_id), in); -// } +void ats_set_scalar(ELM_ATSDriver_ptr ats, int var_id, double in) +{ + reinterpret_cast(ats)->setScalar(static_cast(var_id), in); +} void ats_get_field(ELM_ATSDriver_ptr ats, int var_id, double * const in) { - reinterpret_cast(ats)->getField(reinterpret_cast(var_id), in); + reinterpret_cast(ats)->getField(static_cast(var_id), in); } double const * ats_get_field_ptr(ELM_ATSDriver_ptr ats, int var_id) { - return reinterpret_cast(ats)->getFieldPtr(reinterpret_cast(var_id)); + return reinterpret_cast(ats)->getFieldPtr(static_cast(var_id)); } double * ats_get_field_ptr_w(ELM_ATSDriver_ptr ats, int var_id) { - return reinterpret_cast(ats)->getFieldPtrW(reinterpret_cast(var_id)); + return reinterpret_cast(ats)->getFieldPtrW(static_cast(var_id)); } void ats_set_field(ELM_ATSDriver_ptr ats, int var_id, double const * const in) { - reinterpret_cast(ats)->SetField(reinterpret_cast(var_id), in); + reinterpret_cast(ats)->setField(static_cast(var_id), in); } diff --git a/src/executables/elm_ats_api/elm_ats_driver.cc b/src/executables/elm_ats_api/elm_ats_driver.cc index d20db03d3..89052118a 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.cc +++ b/src/executables/elm_ats_api/elm_ats_driver.cc @@ -21,7 +21,7 @@ #include "CompositeVector.hh" #include "IO.hh" #include "UnstructuredObservations.hh" -#include "pk_helpers.hh" +#include "PK_Helpers.hh" #include "elm_ats_driver.hh" namespace ATS { @@ -80,27 +80,27 @@ ELM_ATSDriver::ELM_ATSDriver(const Teuchos::RCP& plist, // potential sources gross_water_source_key_ = Keys::readKey(*plist_, domain_surf_, "gross water source", "gross_water_source"); - key_map_[ELM::SURFACE_WATER_SOURCE] = gross_water_source_key_; + key_map_[ELM::VarID::GROSS_SURFACE_WATER_SOURCE] = gross_water_source_key_; pot_evap_key_ = Keys::readKey(*plist_, domain_surf_, "potential evaporation mps", "potential_evaporation_mps"); - key_map_[ELM::POTENTIAL_EVAPORATION] = pot_evap_key_; + key_map_[ELM::VarID::POTENTIAL_EVAPORATION] = pot_evap_key_; pot_trans_key_ = Keys::readKey(*plist_, domain_surf_, "potential transpiration mps", "potential_transpiration_mps"); - key_map_[ELM::POTENTIAL_TRANSPIRATION] = pot_trans_key_; + key_map_[ELM::VarID::POTENTIAL_TRANSPIRATION] = pot_trans_key_; // water state surf_wc_key_ = Keys::readKey(*plist_, domain_surf_, "water content", "water_content"); - key_map_[ELM::SURFACE_WATER_CONTENT] = surface_wc_key_; + key_map_[ELM::VarID::SURFACE_WATER_CONTENT] = surf_wc_key_; wc_key_ = Keys::readKey(*plist_, domain_subsurf_, "water content", "water_content"); - key_map_[ELM::WATER_CONTENT] = wc_key_; + key_map_[ELM::VarID::WATER_CONTENT] = wc_key_; // actual water fluxes evap_key_ = Keys::readKey(*plist_, domain_surf_, "evaporation", "evaporation"); - key_map_[ELM::EVAPORATION] = evap_key_; + key_map_[ELM::VarID::EVAPORATION] = evap_key_; col_trans_key_ = Keys::readKey(*plist_, domain_surf_, "total transpiration", "total_transpiration"); - key_map_[ELM::TRANSPIRATION] = col_trans_key_; + key_map_[ELM::VarID::TRANSPIRATION] = col_trans_key_; col_baseflow_key_ = Keys::readKey(*plist_, domain_surf_, "column total baseflow", "baseflow"); - key_map_[ELM::BASEFLOW] = col_baseflow_key_; + key_map_[ELM::VarID::BASEFLOW] = col_baseflow_key_; col_runoff_key_ = Keys::readKey(*plist_, domain_surf_, "runoff generation", "runoff_generation"); - key_map_[ELM::RUNOFF] = col_runoff_key_; + key_map_[ELM::VarID::RUNOFF] = col_runoff_key_; // 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"); @@ -123,51 +123,49 @@ void ELM_ATSDriver::setup() { // potential fluxes (ELM -> ATS) - requireAtNext(gross_water_source_key_, Amanzi::Tags::NEXT, *S_, gross_water_source_key_) + ATS::requireEvaluatorAtNext(gross_water_source_key_, Amanzi::Tags::NEXT, *S_, gross_water_source_key_) .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - requireAtNext(pot_evap_key_, Amanzi::Tags::NEXT, *S_, pot_evap_key_) + ATS::requireEvaluatorAtNext(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_) + ATS::requireEvaluatorAtNext(pot_trans_key_, Amanzi::Tags::NEXT, *S_, pot_trans_key_) .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - requireAtNext(surf_cv_key_, Amanzi::Tags::NEXT, *S_) + ATS::requireEvaluatorAtNext(surf_cv_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - requireAtNext(cv_key_, Amanzi::Tags::NEXT, *S_) + ATS::requireEvaluatorAtNext(cv_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); - requireAtNext(surf_mol_dens_key_, Amanzi::Tags::NEXT, *S_) + ATS::requireEvaluatorAtNext(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_) + ATS::requireEvaluatorAtNext(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_) + ATS::requireEvaluatorAtNext(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_) + ATS::requireEvaluatorAtNext(subsurf_mass_dens_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); - requireAtNext(evap_key_, Amanzi::Tags::NEXT, *S_) + ATS::requireEvaluatorAtNext(evap_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - requireAtNext(col_trans_key_, Amanzi::Tags::NEXT, *S_) + ATS::requireEvaluatorAtNext(col_trans_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - requireAtNext(col_baseflow_key_, Amanzi::Tags::NEXT, *S_, col_baseflow_key_) + ATS::requireEvaluatorAtNext(col_baseflow_key_, Amanzi::Tags::NEXT, *S_, col_baseflow_key_) .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); - requireAtNext(col_runoff_key_, Amanzi::Tags::NEXT, *S_, col_runoff_key_) + ATS::requireEvaluatorAtNext(col_runoff_key_, Amanzi::Tags::NEXT, *S_, col_runoff_key_) .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); Coordinator::setup(); - requireAtNext(pres_key_, Amanzi::Tags::NEXT, *S_) + ATS::requireEvaluatorAtNext(pres_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); - requireAtNext(surf_pres_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); } // // dz & depth are currently ignored -- presumed identical between ATS & ELM // -MeshInfo ELM_ATSDriver::getMeshInfo() +ELM_ATSDriver::MeshInfo ELM_ATSDriver::getMeshInfo() { - MeshInfo info; + 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_; @@ -175,12 +173,10 @@ MeshInfo ELM_ATSDriver::getMeshInfo() } -void ELM_ATSDriver::initialize(double t, - double const * const elm_water_content, - double const * const elm_surf_water_content) +void ELM_ATSDriver::initialize() { // Assign start time to ATS - t0_ = t; + t0_ = S_->get_time(); // initialize ATS data, commit initial conditions Coordinator::initialize(); @@ -191,8 +187,8 @@ void ELM_ATSDriver::initialize(double t, initZero_(pot_evap_key_); // no current evaluators for these, treat as primary - initZero_(baseflow_key_); - initZero_(runoff_key_); + initZero_(col_baseflow_key_); + initZero_(col_runoff_key_); // visualization at IC -- TODO remove this or place behind flag visualize(); @@ -223,7 +219,7 @@ void ELM_ATSDriver::advance(double dt, bool do_chkp, bool do_vis) // 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 @@ -252,9 +248,10 @@ void ELM_ATSDriver::initZero_(const Key& key) -void ELM_ATSDriver::copyToSurf_(double const * const in, const Key& key, Key owner) +void ELM_ATSDriver::copyToSurf_(double const * const in, const Key& key) { - if (owner.empty()) owner = key; + Key owner = S_->GetRecord(key, Amanzi::Tags::NEXT).owner(); + // surf maps directly into columns auto& vec = *S_->GetW(key, Amanzi::Tags::NEXT, owner) .ViewComponent("cell", false); @@ -270,9 +267,9 @@ void ELM_ATSDriver::copyToSurf_(double const * const in, const Key& key, Key own // 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) +void ELM_ATSDriver::copyToSub_(double const * const in, const Key& key) { - if (owner.empty()) owner = key; + Key owner = S_->GetRecord(key, Amanzi::Tags::NEXT).owner(); auto& vec = *S_->GetW(key, Amanzi::Tags::NEXT, owner) .ViewComponent("cell", false); @@ -316,7 +313,7 @@ void ELM_ATSDriver::copyFromSub_(double * const out, const Key& key) const } -void ELM_ATSDriver::getField(const VarID& var_id, double * const out) +void ELM_ATSDriver::getField(const ELM::VarID& var_id, double * const out) { switch(var_id) { case ELM::VarID::SURFACE_WATER_CONTENT: getField_(out); break; @@ -326,11 +323,12 @@ void ELM_ATSDriver::getField(const VarID& var_id, double * const out) case ELM::VarID::RUNOFF: getField_(out); break; case ELM::VarID::BASEFLOW: getField_(out); break; default: throw std::runtime_error("Ungettable variable"); + } } -void ELM_ATSDriver::setField(const VarID& var_id, double const * const out) +void ELM_ATSDriver::setField(const ELM::VarID& var_id, double const * const out) { switch(var_id) { case ELM::VarID::BASE_POROSITY: setField_(out); break; @@ -342,17 +340,18 @@ void ELM_ATSDriver::setField(const VarID& var_id, double const * const out) case ELM::VarID::ROOT_FRACTION: setField_(out); break; case ELM::VarID::SURFACE_WATER_CONTENT: setField_(out); break; case ELM::VarID::WATER_CONTENT: setField_(out); break; - case ELM::VarID::SURFACE_WATER_SOURCE: setField_(out); break; + case ELM::VarID::GROSS_SURFACE_WATER_SOURCE: setField_(out); break; case ELM::VarID::POTENTIAL_EVAPORATION: setField_(out); break; case ELM::VarID::POTENTIAL_TRANSPIRATION: setField_(out); break; default: throw std::runtime_error("Unsettable variable"); + } } template void ELM_ATSDriver::getField_(double * const out) { - Key field_key = keys_[var_id]; + Key field_key = key_map_[var_id]; if (Keys::getDomain(field_key) == domain_subsurf_) { copyFromSub_(out, field_key); } else { @@ -364,7 +363,7 @@ void ELM_ATSDriver::getField_(double * const out) template void ELM_ATSDriver::setField_(double const * const out) { - Key field_key = keys_[var_id]; + Key field_key = key_map_[var_id]; if (Keys::getDomain(field_key) == domain_subsurf_) { copyToSub_(out, field_key); } else { diff --git a/src/executables/elm_ats_api/elm_ats_driver.hh b/src/executables/elm_ats_api/elm_ats_driver.hh index f6d14bb0a..b5bd1f9e4 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.hh +++ b/src/executables/elm_ats_api/elm_ats_driver.hh @@ -16,6 +16,8 @@ #include "Key.hh" #include "coordinator.hh" +#include "ats_variables.hh" + namespace ATS { using namespace Amanzi; @@ -43,20 +45,20 @@ class ELM_ATSDriver : public Coordinator { void advanceTest(); void finalize(); - void setScalar(const ScalarID& scalar_id, double in); - double getScalar(const ScalarID& scalar_id); + void setScalar(const ELM::VarID& scalar_id, double in); + double getScalar(const ELM::VarID& scalar_id); - void setField(const VarID& var_id, double const * const in); - void getField(const VarID& var_id, double * const out); - double const * getFieldPtr(const VarID& var_id); - double * getFieldPtrW(const VarID& var_id); + void setField(const ELM::VarID& var_id, double const * const in); + void getField(const ELM::VarID& var_id, double * const out); + double const * getFieldPtr(const ELM::VarID& var_id); + double * getFieldPtrW(const ELM::VarID& var_id); private: - template + template void setField_(double const * const in); - template + template void getField_(double * const out); void copyToSurf_(double const * const in, const Key& key); @@ -83,6 +85,7 @@ class ELM_ATSDriver : public Coordinator { Key pot_evap_key_; Key pot_trans_key_; + Key pres_key_; Key wc_key_; Key surf_wc_key_; @@ -99,6 +102,8 @@ class ELM_ATSDriver : public Coordinator { Key surf_cv_key_; Key cv_key_; + std::map key_map_; + }; From f56a6b79170874840787cd5ff07e88970a3381a8 Mon Sep 17 00:00:00 2001 From: Ethan Coon Date: Tue, 14 Oct 2025 12:00:46 -0600 Subject: [PATCH 05/16] cleanup as things compiile --- src/executables/elm_ats_api/elm_ats_api.h | 4 +-- src/executables/elm_ats_api/elm_ats_driver.cc | 30 +++++++++++++++++++ 2 files changed, 32 insertions(+), 2 deletions(-) diff --git a/src/executables/elm_ats_api/elm_ats_api.h b/src/executables/elm_ats_api/elm_ats_api.h index 03db2e838..8039175eb 100644 --- a/src/executables/elm_ats_api/elm_ats_api.h +++ b/src/executables/elm_ats_api/elm_ats_api.h @@ -75,8 +75,8 @@ void ats_advance_test(ELM_ATSDriver_ptr ats); // // Memory movement/data passing // ----------------------------------------------------------------------------- -// double ats_get_scalar(ELM_ATSDriver_ptr ats, int scalar_id); -// void ats_set_scalar(ELM_ATSDriver_ptr ats, int scalar_id, double in); +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); diff --git a/src/executables/elm_ats_api/elm_ats_driver.cc b/src/executables/elm_ats_api/elm_ats_driver.cc index 89052118a..7af5c70ea 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.cc +++ b/src/executables/elm_ats_api/elm_ats_driver.cc @@ -313,6 +313,19 @@ void ELM_ATSDriver::copyFromSub_(double * const out, const Key& key) const } +void ELM_ATSDriver::setScalar(const ELM::VarID& scalar_id, double in) +{ + Key scalar_key = key_map_[scalar_id]; + S_->GetW(scalar_key, Tags::DEFAULT, scalar_key) = in; +} + +double ELM_ATSDriver::getScalar(const ELM::VarID& scalar_id) +{ + Key scalar_key = key_map_[scalar_id]; + return S_->Get(scalar_key, Tags::DEFAULT); +} + + void ELM_ATSDriver::getField(const ELM::VarID& var_id, double * const out) { switch(var_id) { @@ -348,6 +361,23 @@ void ELM_ATSDriver::setField(const ELM::VarID& var_id, double const * const out) } +double const * +ELM_ATSDriver::getFieldPtr(const ELM::VarID& var_id) +{ + Key var_key = key_map_[var_id]; + return &(*S_->Get(var_key, Tags::NEXT).ViewComponent("cell", false))[0][0]; +} + +double * +ELM_ATSDriver::getFieldPtrW(const ELM::VarID& var_id) +{ + Key var_key = key_map_[var_id]; + Key owner = S_->GetRecord(var_key, Amanzi::Tags::NEXT).owner(); + return &(*S_->GetW(var_key, Tags::NEXT, owner).ViewComponent("cell", false))[0][0]; +} + + + template void ELM_ATSDriver::getField_(double * const out) { From 1a4d5deca8ab5083a2e1f4190a84c815d33b61dc Mon Sep 17 00:00:00 2001 From: Ethan Coon Date: Mon, 20 Oct 2025 14:08:13 -0600 Subject: [PATCH 06/16] fixes problem with signaling vs quiet NaN, splits parseParameterList() from setup() in coordinator --- src/executables/ats_driver.cc | 1 + src/executables/coordinator.cc | 15 ++- src/executables/coordinator.hh | 1 + src/executables/elm_ats_api/elm_ats_api.cc | 7 ++ src/executables/elm_ats_api/elm_ats_api.h | 6 ++ src/executables/elm_ats_api/elm_ats_driver.cc | 34 ++++--- src/executables/elm_ats_api/elm_ats_driver.hh | 1 + .../land_cover/LandCover.cc | 96 +++++++++++++------ 8 files changed, 115 insertions(+), 46 deletions(-) diff --git a/src/executables/ats_driver.cc b/src/executables/ats_driver.cc index 498569c18..4609a3a15 100644 --- a/src/executables/ats_driver.cc +++ b/src/executables/ats_driver.cc @@ -69,6 +69,7 @@ ATSDriver::cycle_driver() << "Beginning setup stage..." << std::endl << std::flush; } + parseParameterList(); setup(); if (vo_->os_OK(Teuchos::VERB_LOW)) { *vo_->os() << " ... completed: "; diff --git a/src/executables/coordinator.cc b/src/executables/coordinator.cc index b2023c8e6..f643a175c 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,10 +273,18 @@ Coordinator::Coordinator(const Teuchos::RCP& plist, } +void +Coordinator::parseParameterList() +{ + Teuchos::TimeMonitor monitor(*timers_.at("2a: parseParameterList")); + 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"); @@ -288,8 +297,6 @@ Coordinator::setup() // 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/elm_ats_api.cc b/src/executables/elm_ats_api/elm_ats_api.cc index c3f18c6fd..389f12d56 100644 --- a/src/executables/elm_ats_api/elm_ats_api.cc +++ b/src/executables/elm_ats_api/elm_ats_api.cc @@ -51,6 +51,13 @@ void ats_delete(ELM_ATSDriver_ptr ats) } +// call driver setup() +void ats_parse_parameter_list(ELM_ATSDriver_ptr ats) +{ + reinterpret_cast(ats)->parseParameterList(); +} + + void ats_get_mesh_info(ELM_ATSDriver_ptr ats, int * const ncols_local, int * const ncols_global, diff --git a/src/executables/elm_ats_api/elm_ats_api.h b/src/executables/elm_ats_api/elm_ats_api.h index 8039175eb..9215bf99f 100644 --- a/src/executables/elm_ats_api/elm_ats_api.h +++ b/src/executables/elm_ats_api/elm_ats_api.h @@ -42,6 +42,12 @@ void ats_delete(ELM_ATSDriver_ptr ats); // Control // ----------------------------------------------------------------------------- +// +// simulation setup +// +void ats_parse_parameter_list(ELM_ATSDriver_ptr ats); + + // // These quantities should be compared against ELM to ensure consistent setup // diff --git a/src/executables/elm_ats_api/elm_ats_driver.cc b/src/executables/elm_ats_api/elm_ats_driver.cc index 7af5c70ea..ad2f7689c 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.cc +++ b/src/executables/elm_ats_api/elm_ats_driver.cc @@ -69,7 +69,13 @@ ELM_ATSDriver::ELM_ATSDriver(const Teuchos::RCP& plist, // -- set default verbosity level to no output // -- TODO make the verbosity level an input argument VerboseObject::global_default_level = Teuchos::VERB_NONE; +} + +void +ELM_ATSDriver::parseParameterList() +{ + // parse my parameter list // domain names domain_subsurf_ = Keys::readDomain(*plist_, "domain", "domain"); domain_surf_ = Keys::readDomainHint(*plist_, domain_subsurf_, "subsurface", "surface"); @@ -116,20 +122,24 @@ ELM_ATSDriver::ELM_ATSDriver(const Teuchos::RCP& plist, ncolumns_ = mesh_surf_->getNumEntities(AmanziMesh::CELL, AmanziMesh::Parallel_kind::OWNED); auto col_zero = mesh_subsurf_->columns.getCells(0); ncells_per_col_ = col_zero.size(); -} - -void -ELM_ATSDriver::setup() -{ + // require my primaries // potential fluxes (ELM -> ATS) ATS::requireEvaluatorAtNext(gross_water_source_key_, Amanzi::Tags::NEXT, *S_, gross_water_source_key_) - .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); ATS::requireEvaluatorAtNext(pot_evap_key_, Amanzi::Tags::NEXT, *S_, pot_evap_key_) - .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); ATS::requireEvaluatorAtNext(pot_trans_key_, Amanzi::Tags::NEXT, *S_, pot_trans_key_) - .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); + + Coordinator::parseParameterList(); +} + +void +ELM_ATSDriver::setup() +{ + // and now require our output variables ATS::requireEvaluatorAtNext(surf_cv_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); ATS::requireEvaluatorAtNext(cv_key_, Amanzi::Tags::NEXT, *S_) @@ -153,10 +163,12 @@ ELM_ATSDriver::setup() ATS::requireEvaluatorAtNext(col_runoff_key_, Amanzi::Tags::NEXT, *S_, col_runoff_key_) .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); - Coordinator::setup(); + // ATS::requireEvaluatorAtNext(pres_key_, Amanzi::Tags::NEXT, *S_) + // .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); - ATS::requireEvaluatorAtNext(pres_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_subsurf_)->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(); } diff --git a/src/executables/elm_ats_api/elm_ats_driver.hh b/src/executables/elm_ats_api/elm_ats_driver.hh index b5bd1f9e4..5a03268d1 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.hh +++ b/src/executables/elm_ats_api/elm_ats_driver.hh @@ -39,6 +39,7 @@ class ELM_ATSDriver : public Coordinator { int npfts = 17); MeshInfo getMeshInfo(); + void parseParameterList(); void setup(); void initialize(); void advance(double dt, bool checkpoint, bool vis); 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) {} From 330881925d8c1c048fc5fd4fb032c1b49e672281 Mon Sep 17 00:00:00 2001 From: Ethan Coon Date: Tue, 21 Oct 2025 13:32:40 -0600 Subject: [PATCH 07/16] better logging -- Case 0 --- src/executables/ats_driver.cc | 1 - src/executables/elm_ats_api/elm_ats_driver.cc | 114 ++++++++++++++---- 2 files changed, 88 insertions(+), 27 deletions(-) diff --git a/src/executables/ats_driver.cc b/src/executables/ats_driver.cc index 4609a3a15..48a286348 100644 --- a/src/executables/ats_driver.cc +++ b/src/executables/ats_driver.cc @@ -183,7 +183,6 @@ ATSDriver::cycle_driver() reportOneTimer_("4: solve"); } - // finalizing simulation if (vo_->os_OK(Teuchos::VERB_LOW)) { *vo_->os() << "================================================================================" diff --git a/src/executables/elm_ats_api/elm_ats_driver.cc b/src/executables/elm_ats_api/elm_ats_driver.cc index ad2f7689c..dbd842929 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.cc +++ b/src/executables/elm_ats_api/elm_ats_driver.cc @@ -68,13 +68,19 @@ ELM_ATSDriver::ELM_ATSDriver(const Teuchos::RCP& plist, { // -- 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; } 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"); @@ -84,6 +90,8 @@ ELM_ATSDriver::parseParameterList() mesh_subsurf_ = S_->GetMesh(domain_subsurf_); mesh_surf_ = S_->GetMesh(domain_surf_); + key_map_[ELM::VarID::TIME] = "time"; + // potential sources gross_water_source_key_ = Keys::readKey(*plist_, domain_surf_, "gross water source", "gross_water_source"); key_map_[ELM::VarID::GROSS_SURFACE_WATER_SOURCE] = gross_water_source_key_; @@ -133,12 +141,23 @@ ELM_ATSDriver::parseParameterList() .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); Coordinator::parseParameterList(); + if (vo_->os_OK(Teuchos::VERB_LOW)) { + *vo_->os() << " ... completed: "; + reportOneTimer_("2a: parseParameterList"); + } } void ELM_ATSDriver::setup() { + if (vo_->os_OK(Teuchos::VERB_LOW)) { + *vo_->os() << "================================================================================" + << std::endl + << "Beginning setup stage..." << std::endl + << std::flush; + } + // and now require our output variables ATS::requireEvaluatorAtNext(surf_cv_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); @@ -169,6 +188,10 @@ ELM_ATSDriver::setup() // 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"); + } } @@ -187,8 +210,12 @@ ELM_ATSDriver::MeshInfo ELM_ATSDriver::getMeshInfo() void ELM_ATSDriver::initialize() { - // Assign start time to ATS - t0_ = S_->get_time(); + if (vo_->os_OK(Teuchos::VERB_LOW)) { + *vo_->os() << "================================================================================" + << std::endl + << "Beginning initialize stage..." << std::endl + << std::flush; + } // initialize ATS data, commit initial conditions Coordinator::initialize(); @@ -205,28 +232,56 @@ void ELM_ATSDriver::initialize() // visualization at IC -- TODO remove this or place behind flag visualize(); checkpoint(); + + if (vo_->os_OK(Teuchos::VERB_LOW)) { + *vo_->os() << " ... completed: "; + reportOneTimer_("3: initialize"); + } } void ELM_ATSDriver::advance(double dt, bool do_chkp, bool do_vis) { - S_->Assign("dt", Amanzi::Tags::DEFAULT, "dt", dt); - S_->advance_time(Amanzi::Tags::NEXT, dt); - - // solve model for a single timestep - bool fail = Coordinator::advance(); - if (fail) { - Errors::Message msg("ELM_ATSDriver: advance(dt) failed. Make ATS subcycle for proper ELM use."); - Exceptions::amanzi_throw(msg); - } + { + 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); + + // solve model for a single timestep + bool fail = Coordinator::advance(); + if (fail) { + Errors::Message msg("ELM_ATSDriver: advance(dt) failed. Make ATS subcycle for proper ELM use."); + Exceptions::amanzi_throw(msg); + } - // make observations, vis, and checkpoints - for (const auto& obs : observations_) obs->MakeObservations(S_.ptr()); + S_->set_time(Amanzi::Tags::CURRENT, S_->get_time(Amanzi::Tags::NEXT)); + S_->advance_cycle(); - // vis/checkpoint if EITHER ATS or ELM request it - if (do_vis && !visualize()) visualize(true); - if (do_chkp && !checkpoint()) checkpoint(true); + // vis/checkpoint if EITHER ATS or ELM request it + if (do_vis && !visualize()) visualize(true); + if (do_chkp && !checkpoint()) checkpoint(true); + observe(); + } + if (vo_->os_OK(Teuchos::VERB_LOW)) { + *vo_->os() << " ... completed: "; + reportOneTimer_("4: solve"); + } } @@ -243,10 +298,17 @@ void ELM_ATSDriver::advanceTest() void ELM_ATSDriver::finalize() { - WriteStateStatistics(*S_, *vo_); - report_memory(); - Teuchos::TimeMonitor::summarize(*vo_->os()); + if (vo_->os_OK(Teuchos::VERB_LOW)) { + *vo_->os() << "================================================================================" + << std::endl + << "Beginning finalize stage..." << std::endl + << std::flush; + } Coordinator::finalize(); + if (vo_->os_OK(Teuchos::VERB_LOW)) { + *vo_->os() << " ... completed: "; + reportOneTimer_("5: finalize"); + } } @@ -327,13 +389,13 @@ void ELM_ATSDriver::copyFromSub_(double * const out, const Key& key) const void ELM_ATSDriver::setScalar(const ELM::VarID& scalar_id, double in) { - Key scalar_key = key_map_[scalar_id]; + Key scalar_key = key_map_.at(scalar_id); S_->GetW(scalar_key, Tags::DEFAULT, scalar_key) = in; } double ELM_ATSDriver::getScalar(const ELM::VarID& scalar_id) { - Key scalar_key = key_map_[scalar_id]; + Key scalar_key = key_map_.at(scalar_id); return S_->Get(scalar_key, Tags::DEFAULT); } @@ -376,14 +438,14 @@ void ELM_ATSDriver::setField(const ELM::VarID& var_id, double const * const out) double const * ELM_ATSDriver::getFieldPtr(const ELM::VarID& var_id) { - Key var_key = key_map_[var_id]; + Key var_key = key_map_.at(var_id); return &(*S_->Get(var_key, Tags::NEXT).ViewComponent("cell", false))[0][0]; } double * ELM_ATSDriver::getFieldPtrW(const ELM::VarID& var_id) { - Key var_key = key_map_[var_id]; + Key var_key = key_map_.at(var_id); Key owner = S_->GetRecord(var_key, Amanzi::Tags::NEXT).owner(); return &(*S_->GetW(var_key, Tags::NEXT, owner).ViewComponent("cell", false))[0][0]; } @@ -393,7 +455,7 @@ ELM_ATSDriver::getFieldPtrW(const ELM::VarID& var_id) template void ELM_ATSDriver::getField_(double * const out) { - Key field_key = key_map_[var_id]; + Key field_key = key_map_.at(var_id); if (Keys::getDomain(field_key) == domain_subsurf_) { copyFromSub_(out, field_key); } else { @@ -405,7 +467,7 @@ void ELM_ATSDriver::getField_(double * const out) template void ELM_ATSDriver::setField_(double const * const out) { - Key field_key = key_map_[var_id]; + Key field_key = key_map_.at(var_id); if (Keys::getDomain(field_key) == domain_subsurf_) { copyToSub_(out, field_key); } else { From 6ead7f8fa7660fbb27c8f42087719e5690c251ab Mon Sep 17 00:00:00 2001 From: Ethan Coon Date: Tue, 21 Oct 2025 22:03:37 -0600 Subject: [PATCH 08/16] adds logfile passing, adds saturation/ponded depth rather than water content, removes unused keys --- src/executables/elm_ats_api/ats_variables.F90 | 3 + src/executables/elm_ats_api/ats_variables.hh | 5 +- src/executables/elm_ats_api/elm_ats_api.cc | 25 +- src/executables/elm_ats_api/elm_ats_api.h | 5 +- src/executables/elm_ats_api/elm_ats_driver.cc | 273 ++++++------------ src/executables/elm_ats_api/elm_ats_driver.hh | 39 +-- 6 files changed, 128 insertions(+), 222 deletions(-) diff --git a/src/executables/elm_ats_api/ats_variables.F90 b/src/executables/elm_ats_api/ats_variables.F90 index ac99d294e..992c8df8b 100644 --- a/src/executables/elm_ats_api/ats_variables.F90 +++ b/src/executables/elm_ats_api/ats_variables.F90 @@ -23,6 +23,9 @@ module ats_variables 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() diff --git a/src/executables/elm_ats_api/ats_variables.hh b/src/executables/elm_ats_api/ats_variables.hh index d8a4c1155..c71e726e6 100644 --- a/src/executables/elm_ats_api/ats_variables.hh +++ b/src/executables/elm_ats_api/ats_variables.hh @@ -21,7 +21,10 @@ enum class VarID : int { RUNOFF = 16, BASEFLOW = 17, ELEVATION = 18, - TIME = 19 + TIME = 19, + INITIAL_WATER_CONTENT = 20, + SATURATION_LIQUID = 21, + PONDED_DEPTH = 22 }; } // namespace ELM diff --git a/src/executables/elm_ats_api/elm_ats_api.cc b/src/executables/elm_ats_api/elm_ats_api.cc index 389f12d56..3f271fc04 100644 --- a/src/executables/elm_ats_api/elm_ats_api.cc +++ b/src/executables/elm_ats_api/elm_ats_api.cc @@ -25,14 +25,16 @@ 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) +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)); + return reinterpret_cast(ATS::createELM_ATSDriver(f_comm, input_filename, logfile_filename)); } @@ -62,12 +64,19 @@ void ats_get_mesh_info(ELM_ATSDriver_ptr ats, int * const ncols_local, int * const ncols_global, int * const nlevgrnd, - double * const dzs) + double * const dzs, + double * const areas) { 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]; + } + for (int i = 0; i != info.ncols_local; ++i) { + areas[i] = info.areas[i]; + } } @@ -112,11 +121,6 @@ void ats_set_scalar(ELM_ATSDriver_ptr ats, int var_id, double in) reinterpret_cast(ats)->setScalar(static_cast(var_id), in); } -void ats_get_field(ELM_ATSDriver_ptr ats, int var_id, double * const in) -{ - reinterpret_cast(ats)->getField(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)); @@ -128,11 +132,6 @@ double * ats_get_field_ptr_w(ELM_ATSDriver_ptr ats, int var_id) return reinterpret_cast(ats)->getFieldPtrW(static_cast(var_id)); } -void ats_set_field(ELM_ATSDriver_ptr ats, int var_id, double const * const in) -{ - reinterpret_cast(ats)->setField(static_cast(var_id), in); -} - #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 9215bf99f..d4d9af497 100644 --- a/src/executables/elm_ats_api/elm_ats_api.h +++ b/src/executables/elm_ats_api/elm_ats_api.h @@ -32,7 +32,7 @@ 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); @@ -55,7 +55,8 @@ void ats_get_mesh_info(ELM_ATSDriver_ptr ats, int * const ncols_local, int * const ncols_global, int * const nlevgrnd, - double * const dzs); + double * const dzs, + double * const areas); // // simulation setup diff --git a/src/executables/elm_ats_api/elm_ats_driver.cc b/src/executables/elm_ats_api/elm_ats_driver.cc index dbd842929..aea28430b 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.cc +++ b/src/executables/elm_ats_api/elm_ats_driver.cc @@ -27,7 +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); @@ -38,6 +38,7 @@ 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()) { @@ -52,7 +53,7 @@ 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); } @@ -60,6 +61,7 @@ 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), @@ -69,6 +71,7 @@ ELM_ATSDriver::ELM_ATSDriver(const Teuchos::RCP& plist, // -- set default verbosity level to no output // -- TODO make the verbosity level an input argument VerboseObject::global_default_level = Teuchos::VERB_HIGH; + VerboseObject::global_logfile = logfile; } @@ -90,41 +93,47 @@ ELM_ATSDriver::parseParameterList() mesh_subsurf_ = S_->GetMesh(domain_subsurf_); mesh_surf_ = S_->GetMesh(domain_surf_); - key_map_[ELM::VarID::TIME] = "time"; + key_map_[ELM::VarID::TIME] = {"time", Tags::NEXT}; + + // parameters + poro_key_ = Keys::readKey(*plist_, domain_subsurf_, "porosity", "porosity"); + key_map_[ELM::VarID::EFFECTIVE_POROSITY] = { poro_key_, Tags::NEXT }; // potential sources gross_water_source_key_ = Keys::readKey(*plist_, domain_surf_, "gross water source", "gross_water_source"); - key_map_[ELM::VarID::GROSS_SURFACE_WATER_SOURCE] = gross_water_source_key_; + key_map_[ELM::VarID::GROSS_SURFACE_WATER_SOURCE] = { gross_water_source_key_, Tags::NEXT }; pot_evap_key_ = Keys::readKey(*plist_, domain_surf_, "potential evaporation mps", "potential_evaporation_mps"); - key_map_[ELM::VarID::POTENTIAL_EVAPORATION] = pot_evap_key_; + key_map_[ELM::VarID::POTENTIAL_EVAPORATION] = { pot_evap_key_, Tags::NEXT }; pot_trans_key_ = Keys::readKey(*plist_, domain_surf_, "potential transpiration mps", "potential_transpiration_mps"); - key_map_[ELM::VarID::POTENTIAL_TRANSPIRATION] = pot_trans_key_; + key_map_[ELM::VarID::POTENTIAL_TRANSPIRATION] = { pot_trans_key_, Tags::NEXT }; // water state - surf_wc_key_ = Keys::readKey(*plist_, domain_surf_, "water content", "water_content"); - key_map_[ELM::VarID::SURFACE_WATER_CONTENT] = surf_wc_key_; - wc_key_ = Keys::readKey(*plist_, domain_subsurf_, "water content", "water_content"); - key_map_[ELM::VarID::WATER_CONTENT] = wc_key_; + pres_key_ = Keys::readKey(*plist_, domain_subsurf_, "pressure", "pressure"); + key_map_[ELM::VarID::PRESSURE] = { pres_key_, Tags::NEXT }; + sat_key_ = Keys::readKey(*plist_, domain_subsurf_, "saturation liquid", "saturation_liquid"); + key_map_[ELM::VarID::SATURATION_LIQUID] = { sat_key_, Tags::NEXT }; + pd_key_ = Keys::readKey(*plist_, domain_surf_, "ponded depth", "ponded_depth"); + key_map_[ELM::VarID::PONDED_DEPTH] = { pd_key_, Tags::NEXT }; // actual water fluxes evap_key_ = Keys::readKey(*plist_, domain_surf_, "evaporation", "evaporation"); - key_map_[ELM::VarID::EVAPORATION] = evap_key_; + key_map_[ELM::VarID::EVAPORATION] = { evap_key_, Tags::NEXT }; col_trans_key_ = Keys::readKey(*plist_, domain_surf_, "total transpiration", "total_transpiration"); - key_map_[ELM::VarID::TRANSPIRATION] = col_trans_key_; + key_map_[ELM::VarID::TRANSPIRATION] = { col_trans_key_, Tags::NEXT }; col_baseflow_key_ = Keys::readKey(*plist_, domain_surf_, "column total baseflow", "baseflow"); - key_map_[ELM::VarID::BASEFLOW] = col_baseflow_key_; + key_map_[ELM::VarID::BASEFLOW] = { col_baseflow_key_, Tags::NEXT }; col_runoff_key_ = Keys::readKey(*plist_, domain_surf_, "runoff generation", "runoff_generation"); - key_map_[ELM::VarID::RUNOFF] = col_runoff_key_; + 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"); + // // 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"); // cell vol keys - surf_cv_key_ = Keys::getKey(domain_surf_, "cell_volume"); - cv_key_ = Keys::getKey(domain_subsurf_, "cell_volume"); + // surf_cv_key_ = Keys::getKey(domain_surf_, "cell_volume"); + // cv_key_ = Keys::getKey(domain_subsurf_, "cell_volume"); // -- mesh info ncolumns_ = mesh_surf_->getNumEntities(AmanziMesh::CELL, AmanziMesh::Parallel_kind::OWNED); @@ -132,6 +141,10 @@ ELM_ATSDriver::parseParameterList() ncells_per_col_ = col_zero.size(); // require my primaries + // parameters set by ELM + ATS::requireEvaluatorAtNext(poro_key_, Amanzi::Tags::NEXT, *S_, poro_key_) + .SetMesh(mesh_subsurf_)->SetComponent("cell", AmanziMesh::CELL, 1); + // potential fluxes (ELM -> ATS) ATS::requireEvaluatorAtNext(gross_water_source_key_, Amanzi::Tags::NEXT, *S_, gross_water_source_key_) .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); @@ -159,18 +172,26 @@ ELM_ATSDriver::setup() } // and now require our output variables - ATS::requireEvaluatorAtNext(surf_cv_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - ATS::requireEvaluatorAtNext(cv_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + // ATS::requireEvaluatorAtNext(surf_cv_key_, Amanzi::Tags::NEXT, *S_) + // .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + // ATS::requireEvaluatorAtNext(cv_key_, Amanzi::Tags::NEXT, *S_) + // .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); - ATS::requireEvaluatorAtNext(surf_mol_dens_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - ATS::requireEvaluatorAtNext(surf_mass_dens_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - ATS::requireEvaluatorAtNext(subsurf_mol_dens_key_, Amanzi::Tags::NEXT, *S_) - .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); - ATS::requireEvaluatorAtNext(subsurf_mass_dens_key_, Amanzi::Tags::NEXT, *S_) + // ATS::requireEvaluatorAtNext(surf_mol_dens_key_, Amanzi::Tags::NEXT, *S_) + // .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + // ATS::requireEvaluatorAtNext(surf_mass_dens_key_, Amanzi::Tags::NEXT, *S_) + // .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + // ATS::requireEvaluatorAtNext(subsurf_mol_dens_key_, Amanzi::Tags::NEXT, *S_) + // .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + // ATS::requireEvaluatorAtNext(subsurf_mass_dens_key_, Amanzi::Tags::NEXT, *S_) + // .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + + // output variables + ATS::requireEvaluatorAtNext(pd_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + ATS::requireEvaluatorAtNext(sat_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + ATS::requireEvaluatorAtNext(pres_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); ATS::requireEvaluatorAtNext(evap_key_, Amanzi::Tags::NEXT, *S_) @@ -182,8 +203,6 @@ ELM_ATSDriver::setup() ATS::requireEvaluatorAtNext(col_runoff_key_, Amanzi::Tags::NEXT, *S_, col_runoff_key_) .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); - // ATS::requireEvaluatorAtNext(pres_key_, Amanzi::Tags::NEXT, *S_) - // .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); // This must be last always -- it allocates memory calling State::setup, so // all other setup must be done. @@ -204,6 +223,20 @@ ELM_ATSDriver::MeshInfo ELM_ATSDriver::getMeshInfo() 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) { + info.dzs[i] = mesh_subsurf_->getCellVolume(i) / surface_area; + } + + // surface area + info.areas.resize(ncolumns_); + for (int i = 0; i != ncolumns_; ++i) { + info.areas[i] = mesh_surf_->getCellVolume(i); + } return info; } @@ -217,17 +250,20 @@ void ELM_ATSDriver::initialize() << std::flush; } - // initialize ATS data, commit initial conditions - Coordinator::initialize(); + // 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_(gross_water_source_key_); - initZero_(pot_trans_key_); - initZero_(pot_evap_key_); + initValue_(gross_water_source_key_); + initValue_(pot_trans_key_); + initValue_(pot_evap_key_); // no current evaluators for these, treat as primary - initZero_(col_baseflow_key_); - initZero_(col_runoff_key_); + 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(); @@ -312,170 +348,43 @@ void ELM_ATSDriver::finalize() } -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 = S_->GetRecord(key, Amanzi::Tags::NEXT).owner(); - - // 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 = S_->GetRecord(key, Amanzi::Tags::NEXT).owner(); - 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]; - } - } - - changedEvaluatorPrimary(key, Amanzi::Tags::NEXT, *S_); -} - - -void ELM_ATSDriver::copyFromSurf_(double * const out, const Key& key) const +void ELM_ATSDriver::initValue_(const Key& key, double value) { - 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 = mesh_subsurf_->columns.getCells(i); - for (int j=0; j!=ncells_per_col_; ++j) { - out[j * ncolumns_ + i] = vec[0][cells[j]]; - } - } + auto& vec = S_->GetW(key, Tags::NEXT, key); + vec.PutScalar(value); + S_->GetRecordW(key, Tags::NEXT, key).set_initialized(); } void ELM_ATSDriver::setScalar(const ELM::VarID& scalar_id, double in) { - Key scalar_key = key_map_.at(scalar_id); - S_->GetW(scalar_key, Tags::DEFAULT, scalar_key) = in; + KeyTag scalar_key = key_map_.at(scalar_id); + S_->GetW(scalar_key.first, scalar_key.second, scalar_key.first) = in; } double ELM_ATSDriver::getScalar(const ELM::VarID& scalar_id) { - Key scalar_key = key_map_.at(scalar_id); - return S_->Get(scalar_key, Tags::DEFAULT); -} - - -void ELM_ATSDriver::getField(const ELM::VarID& var_id, double * const out) -{ - switch(var_id) { - case ELM::VarID::SURFACE_WATER_CONTENT: getField_(out); break; - case ELM::VarID::WATER_CONTENT: getField_(out); break; - case ELM::VarID::EVAPORATION: getField_(out); break; - case ELM::VarID::TRANSPIRATION: getField_(out); break; - case ELM::VarID::RUNOFF: getField_(out); break; - case ELM::VarID::BASEFLOW: getField_(out); break; - default: throw std::runtime_error("Ungettable variable"); - } -} - - - -void ELM_ATSDriver::setField(const ELM::VarID& var_id, double const * const out) -{ - switch(var_id) { - case ELM::VarID::BASE_POROSITY: setField_(out); break; - case ELM::VarID::HYDRAULIC_CONDUCTIVITY: setField_(out); break; - case ELM::VarID::CLAPP_HORNBERGER_B: setField_(out); break; - case ELM::VarID::CLAPP_HORNBERGER_PSI_SAT: setField_(out); break; - case ELM::VarID::RESIDUAL_SATURATION: setField_(out); break; - case ELM::VarID::EFFECTIVE_POROSITY: setField_(out); break; - case ELM::VarID::ROOT_FRACTION: setField_(out); break; - case ELM::VarID::SURFACE_WATER_CONTENT: setField_(out); break; - case ELM::VarID::WATER_CONTENT: setField_(out); break; - case ELM::VarID::GROSS_SURFACE_WATER_SOURCE: setField_(out); break; - case ELM::VarID::POTENTIAL_EVAPORATION: setField_(out); break; - case ELM::VarID::POTENTIAL_TRANSPIRATION: setField_(out); break; - default: throw std::runtime_error("Unsettable variable"); - } + KeyTag scalar_key = key_map_.at(scalar_id); + return S_->Get(scalar_key.first, scalar_key.second); } double const * ELM_ATSDriver::getFieldPtr(const ELM::VarID& var_id) { - Key var_key = key_map_.at(var_id); - return &(*S_->Get(var_key, Tags::NEXT).ViewComponent("cell", false))[0][0]; + Amanzi::KeyTag var_key = key_map_.at(var_id); + return &(*S_->Get(var_key.first, var_key.second).ViewComponent("cell", false))[0][0]; } double * ELM_ATSDriver::getFieldPtrW(const ELM::VarID& var_id) { - Key var_key = key_map_.at(var_id); - Key owner = S_->GetRecord(var_key, Amanzi::Tags::NEXT).owner(); - return &(*S_->GetW(var_key, Tags::NEXT, owner).ViewComponent("cell", false))[0][0]; -} - - + Amanzi::KeyTag var_key = key_map_.at(var_id); + Amanzi::changedEvaluatorPrimary(var_key.first, var_key.second, *S_); -template -void ELM_ATSDriver::getField_(double * const out) -{ - Key field_key = key_map_.at(var_id); - if (Keys::getDomain(field_key) == domain_subsurf_) { - copyFromSub_(out, field_key); - } else { - copyFromSurf_(out, field_key); - } + 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]; } -template -void ELM_ATSDriver::setField_(double const * const out) -{ - Key field_key = key_map_.at(var_id); - if (Keys::getDomain(field_key) == domain_subsurf_) { - copyToSub_(out, field_key); - } else { - copyToSurf_(out, field_key); - } -} - - - - } // namespace ATS diff --git a/src/executables/elm_ats_api/elm_ats_driver.hh b/src/executables/elm_ats_api/elm_ats_driver.hh index 5a03268d1..8454a2dce 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.hh +++ b/src/executables/elm_ats_api/elm_ats_driver.hh @@ -30,12 +30,15 @@ class ELM_ATSDriver : public Coordinator { int ncols_local; int ncols_global; int nlevgrnd; + std::vector dzs; + std::vector areas; }; 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); MeshInfo getMeshInfo(); @@ -49,25 +52,11 @@ class ELM_ATSDriver : public Coordinator { void setScalar(const ELM::VarID& scalar_id, double in); double getScalar(const ELM::VarID& scalar_id); - void setField(const ELM::VarID& var_id, double const * const in); - void getField(const ELM::VarID& var_id, double * const out); double const * getFieldPtr(const ELM::VarID& var_id); double * getFieldPtrW(const ELM::VarID& var_id); private: - template - void setField_(double const * const in); - - - template - void getField_(double * const out); - - void copyToSurf_(double const * const in, const Key& key); - void copyToSub_(double const * const in, const Key& key); - void copyFromSurf_(double * const out, const Key& key) const; - void copyFromSub_(double * const out, const Key& key) const; - - void initZero_(const Key& key); + void initValue_(const Key& key, double value = 0.); private: Teuchos::RCP elm_list_; @@ -82,28 +71,30 @@ class ELM_ATSDriver : public Coordinator { Teuchos::RCP mesh_subsurf_; Teuchos::RCP mesh_surf_; + Key poro_key_; + Key gross_water_source_key_; Key pot_evap_key_; Key pot_trans_key_; Key pres_key_; - Key wc_key_; - Key surf_wc_key_; + Key sat_key_; + Key pd_key_; Key evap_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 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_; + // Key surf_cv_key_; + // Key cv_key_; - std::map key_map_; + std::map key_map_; }; From c01be822f23303517fedfa871f3932596831b96a Mon Sep 17 00:00:00 2001 From: Ethan Coon Date: Tue, 21 Oct 2025 22:14:57 -0600 Subject: [PATCH 09/16] missed API change --- src/executables/elm_ats_api/elm_ats_driver.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/executables/elm_ats_api/elm_ats_driver.hh b/src/executables/elm_ats_api/elm_ats_driver.hh index 8454a2dce..05ee30091 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.hh +++ b/src/executables/elm_ats_api/elm_ats_driver.hh @@ -103,7 +103,7 @@ class ELM_ATSDriver : public Coordinator { // 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); +createELM_ATSDriver(MPI_Fint *f_comm, const char *infile, const char *logfile, int npfts=17); } // namespace ATS From 6d2666963423c2d912b8d33352889cbee0e4ce73 Mon Sep 17 00:00:00 2001 From: Ethan Coon Date: Tue, 28 Oct 2025 20:47:51 -0600 Subject: [PATCH 10/16] adds divergence of fluxes evaluator for discharge, baseflow --- src/constitutive_relations/CMakeLists.txt | 5 ++ .../surface_subsurface_fluxes/CMakeLists.txt | 1 + .../flux_divergence_evaluator.cc | 65 +++++++++++++++++++ .../flux_divergence_evaluator.hh | 59 +++++++++++++++++ .../flux_divergence_evaluator_reg.hh | 20 ++++++ src/executables/elm_ats_api/elm_ats_driver.cc | 16 +++-- src/pks/mpc/mpc_coupled_water.cc | 13 ++-- 7 files changed, 168 insertions(+), 11 deletions(-) create mode 100644 src/constitutive_relations/surface_subsurface_fluxes/flux_divergence_evaluator.cc create mode 100644 src/constitutive_relations/surface_subsurface_fluxes/flux_divergence_evaluator.hh create mode 100644 src/constitutive_relations/surface_subsurface_fluxes/flux_divergence_evaluator_reg.hh 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/elm_ats_api/elm_ats_driver.cc b/src/executables/elm_ats_api/elm_ats_driver.cc index aea28430b..d3879b082 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.cc +++ b/src/executables/elm_ats_api/elm_ats_driver.cc @@ -120,9 +120,9 @@ ELM_ATSDriver::parseParameterList() key_map_[ELM::VarID::EVAPORATION] = { evap_key_, Tags::NEXT }; col_trans_key_ = Keys::readKey(*plist_, domain_surf_, "total transpiration", "total_transpiration"); key_map_[ELM::VarID::TRANSPIRATION] = { col_trans_key_, Tags::NEXT }; - col_baseflow_key_ = Keys::readKey(*plist_, domain_surf_, "column total baseflow", "baseflow"); + col_baseflow_key_ = Keys::readKey(*plist_, domain_surf_, "baseflow generation", "baseflow_mps"); key_map_[ELM::VarID::BASEFLOW] = { col_baseflow_key_, Tags::NEXT }; - col_runoff_key_ = Keys::readKey(*plist_, domain_surf_, "runoff generation", "runoff_generation"); + col_runoff_key_ = Keys::readKey(*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 @@ -198,10 +198,10 @@ ELM_ATSDriver::setup() .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); ATS::requireEvaluatorAtNext(col_trans_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - ATS::requireEvaluatorAtNext(col_baseflow_key_, Amanzi::Tags::NEXT, *S_, col_baseflow_key_) - .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); - ATS::requireEvaluatorAtNext(col_runoff_key_, Amanzi::Tags::NEXT, *S_, col_runoff_key_) - .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); + ATS::requireEvaluatorAtNext(col_baseflow_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + ATS::requireEvaluatorAtNext(col_runoff_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 @@ -373,6 +373,10 @@ double const * ELM_ATSDriver::getFieldPtr(const ELM::VarID& var_id) { 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_)); + } + return &(*S_->Get(var_key.first, var_key.second).ViewComponent("cell", false))[0][0]; } 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); From 541fe173210caa1d35f2de6be0d56a7f55774636 Mon Sep 17 00:00:00 2001 From: Ethan Coon Date: Mon, 3 Nov 2025 19:18:39 -0700 Subject: [PATCH 11/16] continued progress --- src/executables/elm_ats_api/ats_variables.hh | 3 +- src/executables/elm_ats_api/elm_ats_driver.cc | 48 ++++++++++--------- src/executables/elm_ats_api/elm_ats_driver.hh | 3 +- src/pks/flow/overland_pressure_pk.cc | 1 + 4 files changed, 31 insertions(+), 24 deletions(-) diff --git a/src/executables/elm_ats_api/ats_variables.hh b/src/executables/elm_ats_api/ats_variables.hh index c71e726e6..d69bddcd5 100644 --- a/src/executables/elm_ats_api/ats_variables.hh +++ b/src/executables/elm_ats_api/ats_variables.hh @@ -24,7 +24,8 @@ enum class VarID : int { TIME = 19, INITIAL_WATER_CONTENT = 20, SATURATION_LIQUID = 21, - PONDED_DEPTH = 22 + PONDED_DEPTH = 22, + DEPTH_TO_WATER_TABLE = 23 }; } // namespace ELM diff --git a/src/executables/elm_ats_api/elm_ats_driver.cc b/src/executables/elm_ats_api/elm_ats_driver.cc index d3879b082..c6c7484a9 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.cc +++ b/src/executables/elm_ats_api/elm_ats_driver.cc @@ -66,7 +66,8 @@ ELM_ATSDriver::ELM_ATSDriver(const Teuchos::RCP& plist, : Coordinator(plist, wallclock_timer, teuchos_comm, comm), npfts_(npfts), ncolumns_(-1), - ncells_per_col_(-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 @@ -86,8 +87,8 @@ ELM_ATSDriver::parseParameterList() } // 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 mesh_subsurf_ = S_->GetMesh(domain_subsurf_); @@ -96,40 +97,42 @@ ELM_ATSDriver::parseParameterList() key_map_[ELM::VarID::TIME] = {"time", Tags::NEXT}; // parameters - poro_key_ = Keys::readKey(*plist_, domain_subsurf_, "porosity", "porosity"); + poro_key_ = Keys::readKey(*elm_plist_, domain_subsurf_, "porosity", "porosity"); key_map_[ELM::VarID::EFFECTIVE_POROSITY] = { poro_key_, Tags::NEXT }; // potential sources - gross_water_source_key_ = Keys::readKey(*plist_, domain_surf_, "gross water source", "gross_water_source"); + 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(*plist_, domain_surf_, "potential evaporation mps", "potential_evaporation_mps"); + 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(*plist_, domain_surf_, "potential transpiration mps", "potential_transpiration_mps"); + 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 - pres_key_ = Keys::readKey(*plist_, domain_subsurf_, "pressure", "pressure"); + pres_key_ = Keys::readKey(*elm_plist_, domain_subsurf_, "pressure", "pressure"); key_map_[ELM::VarID::PRESSURE] = { pres_key_, Tags::NEXT }; - sat_key_ = Keys::readKey(*plist_, domain_subsurf_, "saturation liquid", "saturation_liquid"); + sat_key_ = Keys::readKey(*elm_plist_, domain_subsurf_, "saturation liquid", "saturation_liquid"); key_map_[ELM::VarID::SATURATION_LIQUID] = { sat_key_, Tags::NEXT }; - pd_key_ = Keys::readKey(*plist_, domain_surf_, "ponded depth", "ponded_depth"); + 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 }; // actual water fluxes - evap_key_ = Keys::readKey(*plist_, domain_surf_, "evaporation", "evaporation"); + evap_key_ = Keys::readKey(*elm_plist_, domain_surf_, "evaporation", "evaporation"); key_map_[ELM::VarID::EVAPORATION] = { evap_key_, Tags::NEXT }; - col_trans_key_ = Keys::readKey(*plist_, domain_surf_, "total transpiration", "total_transpiration"); + 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(*plist_, domain_surf_, "baseflow generation", "baseflow_mps"); + 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(*plist_, domain_surf_, "runoff generation", "runoff_generation_mps"); + 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"); + // surf_mol_dens_key_ = Keys::readKey(*elm_plist_, domain_surf_, "surface molar density", "molar_density_liquid"); + // surf_mass_dens_key_ = Keys::readKey(*elm_plist_, domain_surf_, "surface mass density", "mass_density_liquid"); + // subsurf_mol_dens_key_ = Keys::readKey(*elm_plist_, domain_subsurf_, "molar density", "molar_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"); @@ -189,6 +192,8 @@ ELM_ATSDriver::setup() // output variables ATS::requireEvaluatorAtNext(pd_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + ATS::requireEvaluatorAtNext(zwt_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); ATS::requireEvaluatorAtNext(sat_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); ATS::requireEvaluatorAtNext(pres_key_, Amanzi::Tags::NEXT, *S_) @@ -277,7 +282,7 @@ void ELM_ATSDriver::initialize() -void ELM_ATSDriver::advance(double dt, bool do_chkp, bool do_vis) +void ELM_ATSDriver::advance(double dt, bool force_chkp, bool force_vis) { { Teuchos::TimeMonitor timer(*timers_.at("4: solve")); @@ -309,9 +314,8 @@ void ELM_ATSDriver::advance(double dt, bool do_chkp, bool do_vis) S_->set_time(Amanzi::Tags::CURRENT, S_->get_time(Amanzi::Tags::NEXT)); S_->advance_cycle(); - // vis/checkpoint if EITHER ATS or ELM request it - if (do_vis && !visualize()) visualize(true); - if (do_chkp && !checkpoint()) checkpoint(true); + visualize(force_vis); + checkpoint(force_chkp); observe(); } if (vo_->os_OK(Teuchos::VERB_LOW)) { diff --git a/src/executables/elm_ats_api/elm_ats_driver.hh b/src/executables/elm_ats_api/elm_ats_driver.hh index 05ee30091..2c89bfae0 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.hh +++ b/src/executables/elm_ats_api/elm_ats_driver.hh @@ -59,7 +59,7 @@ class ELM_ATSDriver : public Coordinator { void initValue_(const Key& key, double value = 0.); private: - Teuchos::RCP elm_list_; + Teuchos::RCP elm_plist_; int ncolumns_; int ncells_per_col_; @@ -80,6 +80,7 @@ class ELM_ATSDriver : public Coordinator { Key pres_key_; Key sat_key_; Key pd_key_; + Key zwt_key_; Key evap_key_; Key col_trans_key_; 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_); } From 5e174b99a09e15c1a64977a720996bfc8f7a1192 Mon Sep 17 00:00:00 2001 From: Ethan Coon Date: Wed, 5 Nov 2025 09:15:18 -0700 Subject: [PATCH 12/16] updates debugging passing WC_old to ATS from ELM to allow e.g. freezing of water in ELM --- src/executables/elm_ats_api/ats_variables.hh | 4 +- src/executables/elm_ats_api/elm_ats_driver.cc | 123 ++++++++++++------ src/executables/elm_ats_api/elm_ats_driver.hh | 11 +- 3 files changed, 95 insertions(+), 43 deletions(-) diff --git a/src/executables/elm_ats_api/ats_variables.hh b/src/executables/elm_ats_api/ats_variables.hh index d69bddcd5..f2c13a466 100644 --- a/src/executables/elm_ats_api/ats_variables.hh +++ b/src/executables/elm_ats_api/ats_variables.hh @@ -25,7 +25,9 @@ enum class VarID : int { INITIAL_WATER_CONTENT = 20, SATURATION_LIQUID = 21, PONDED_DEPTH = 22, - DEPTH_TO_WATER_TABLE = 23 + DEPTH_TO_WATER_TABLE = 23, + WATER_CONTENT_OLD = 24, + SURFACE_WATER_CONTENT_OLD = 25 }; } // namespace ELM diff --git a/src/executables/elm_ats_api/elm_ats_driver.cc b/src/executables/elm_ats_api/elm_ats_driver.cc index c6c7484a9..e579296f4 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.cc +++ b/src/executables/elm_ats_api/elm_ats_driver.cc @@ -118,6 +118,11 @@ ELM_ATSDriver::parseParameterList() 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 }; + // actual water fluxes evap_key_ = Keys::readKey(*elm_plist_, domain_surf_, "evaporation", "evaporation"); key_map_[ELM::VarID::EVAPORATION] = { evap_key_, Tags::NEXT }; @@ -128,33 +133,29 @@ ELM_ATSDriver::parseParameterList() 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(*elm_plist_, domain_surf_, "surface molar density", "molar_density_liquid"); + // keys for fields used to convert ELM units to ATS units + 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_mol_dens_key_ = Keys::readKey(*elm_plist_, domain_subsurf_, "molar density", "molar_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"); + surf_cv_key_ = Keys::getKey(domain_surf_, "cell_volume"); + cv_key_ = Keys::getKey(domain_subsurf_, "cell_volume"); // -- 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(); - // require my primaries + // require my primary variables // parameters set by ELM - ATS::requireEvaluatorAtNext(poro_key_, Amanzi::Tags::NEXT, *S_, poro_key_) - .SetMesh(mesh_subsurf_)->SetComponent("cell", AmanziMesh::CELL, 1); + requireEvaluatorAtNext(poro_key_, Amanzi::Tags::NEXT, *S_, poro_key_); // potential fluxes (ELM -> ATS) - ATS::requireEvaluatorAtNext(gross_water_source_key_, Amanzi::Tags::NEXT, *S_, gross_water_source_key_) - .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); - ATS::requireEvaluatorAtNext(pot_evap_key_, Amanzi::Tags::NEXT, *S_, pot_evap_key_) - .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); - ATS::requireEvaluatorAtNext(pot_trans_key_, Amanzi::Tags::NEXT, *S_, pot_trans_key_) - .SetMesh(mesh_surf_)->SetComponent("cell", AmanziMesh::CELL, 1); + 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)) { @@ -174,38 +175,54 @@ ELM_ATSDriver::setup() << std::flush; } - // and now require our output variables - // ATS::requireEvaluatorAtNext(surf_cv_key_, Amanzi::Tags::NEXT, *S_) - // .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - // ATS::requireEvaluatorAtNext(cv_key_, Amanzi::Tags::NEXT, *S_) - // .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); - - // ATS::requireEvaluatorAtNext(surf_mol_dens_key_, Amanzi::Tags::NEXT, *S_) - // .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - // ATS::requireEvaluatorAtNext(surf_mass_dens_key_, Amanzi::Tags::NEXT, *S_) - // .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - // ATS::requireEvaluatorAtNext(subsurf_mol_dens_key_, Amanzi::Tags::NEXT, *S_) - // .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); - // ATS::requireEvaluatorAtNext(subsurf_mass_dens_key_, Amanzi::Tags::NEXT, *S_) - // .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); - - // output variables - ATS::requireEvaluatorAtNext(pd_key_, Amanzi::Tags::NEXT, *S_) + // variables needed for unit changes + requireEvaluatorAtNext(surf_cv_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); + requireEvaluatorAtNext(cv_key_, Amanzi::Tags::NEXT, *S_) + .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); + + requireEvaluatorAtNext(surf_mol_dens_key_, Amanzi::Tags::NEXT, *S_) + .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); - ATS::requireEvaluatorAtNext(zwt_key_, Amanzi::Tags::NEXT, *S_) + requireEvaluatorAtNext(zwt_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - ATS::requireEvaluatorAtNext(sat_key_, Amanzi::Tags::NEXT, *S_) + requireEvaluatorAtNext(sat_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); - ATS::requireEvaluatorAtNext(pres_key_, Amanzi::Tags::NEXT, *S_) + requireEvaluatorAtNext(pres_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); - ATS::requireEvaluatorAtNext(evap_key_, Amanzi::Tags::NEXT, *S_) + requireEvaluatorAtNext(evap_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - ATS::requireEvaluatorAtNext(col_trans_key_, Amanzi::Tags::NEXT, *S_) + requireEvaluatorAtNext(col_trans_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - ATS::requireEvaluatorAtNext(col_baseflow_key_, Amanzi::Tags::NEXT, *S_) + requireEvaluatorAtNext(col_baseflow_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); - ATS::requireEvaluatorAtNext(col_runoff_key_, Amanzi::Tags::NEXT, *S_) + requireEvaluatorAtNext(col_runoff_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1); @@ -304,6 +321,36 @@ void ELM_ATSDriver::advance(double dt, bool force_chkp, bool force_vis) 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]; + } + } + + // -- 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 bool fail = Coordinator::advance(); if (fail) { diff --git a/src/executables/elm_ats_api/elm_ats_driver.hh b/src/executables/elm_ats_api/elm_ats_driver.hh index 2c89bfae0..0ec5ce88c 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.hh +++ b/src/executables/elm_ats_api/elm_ats_driver.hh @@ -77,6 +77,9 @@ class ELM_ATSDriver : public Coordinator { Key pot_evap_key_; Key pot_trans_key_; + Key wc_key_; + Key surf_wc_key_; + Key pres_key_; Key sat_key_; Key pd_key_; @@ -87,13 +90,13 @@ class ELM_ATSDriver : public Coordinator { Key col_baseflow_key_; Key col_runoff_key_; - // Key surf_mol_dens_key_; + Key surf_mol_dens_key_; + Key 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_; + Key surf_cv_key_; + Key cv_key_; std::map key_map_; From 8777622ebd57f3b189a53fc0bf8def82102bd506 Mon Sep 17 00:00:00 2001 From: Ethan Coon Date: Tue, 2 Dec 2025 16:13:20 -0700 Subject: [PATCH 13/16] fixes bug in timer names --- src/executables/ats_driver.cc | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/executables/ats_driver.cc b/src/executables/ats_driver.cc index 48a286348..733b69bc8 100644 --- a/src/executables/ats_driver.cc +++ b/src/executables/ats_driver.cc @@ -70,10 +70,15 @@ ATSDriver::cycle_driver() << 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"); } // From daecf316a2b3436235520337f8a394d6095142d1 Mon Sep 17 00:00:00 2001 From: Ethan Coon Date: Tue, 2 Dec 2025 17:07:15 -0700 Subject: [PATCH 14/16] fixes bug -- time may be needed sooner now --- src/executables/coordinator.cc | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/executables/coordinator.cc b/src/executables/coordinator.cc index f643a175c..60310319c 100644 --- a/src/executables/coordinator.cc +++ b/src/executables/coordinator.cc @@ -277,6 +277,11 @@ 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(); } @@ -290,10 +295,6 @@ Coordinator::setup() 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 From 148f5a83626f8a0752111bca2e2231d508481f8c Mon Sep 17 00:00:00 2001 From: Ethan Coon Date: Tue, 9 Dec 2025 16:26:03 -0700 Subject: [PATCH 15/16] pass lat and lon back to ELM --- src/executables/elm_ats_api/elm_ats_api.cc | 6 +++- src/executables/elm_ats_api/elm_ats_api.h | 4 ++- src/executables/elm_ats_api/elm_ats_driver.cc | 28 +++++++++++++++++++ src/executables/elm_ats_api/elm_ats_driver.hh | 4 +++ 4 files changed, 40 insertions(+), 2 deletions(-) diff --git a/src/executables/elm_ats_api/elm_ats_api.cc b/src/executables/elm_ats_api/elm_ats_api.cc index 3f271fc04..0a01f59d9 100644 --- a/src/executables/elm_ats_api/elm_ats_api.cc +++ b/src/executables/elm_ats_api/elm_ats_api.cc @@ -65,7 +65,9 @@ void ats_get_mesh_info(ELM_ATSDriver_ptr ats, int * const ncols_global, int * const nlevgrnd, double * const dzs, - double * const areas) + double * const areas, + double * const lat, + double * const lon) { ATS::ELM_ATSDriver::MeshInfo info = reinterpret_cast(ats)->getMeshInfo(); *ncols_local = info.ncols_local; @@ -76,6 +78,8 @@ void ats_get_mesh_info(ELM_ATSDriver_ptr ats, } for (int i = 0; i != info.ncols_local; ++i) { areas[i] = info.areas[i]; + lat[i] = info.latitudes[i]; + lon[i] = info.longitudes[i]; } } diff --git a/src/executables/elm_ats_api/elm_ats_api.h b/src/executables/elm_ats_api/elm_ats_api.h index d4d9af497..0e9493d0a 100644 --- a/src/executables/elm_ats_api/elm_ats_api.h +++ b/src/executables/elm_ats_api/elm_ats_api.h @@ -56,7 +56,9 @@ void ats_get_mesh_info(ELM_ATSDriver_ptr ats, int * const ncols_global, int * const nlevgrnd, double * const dzs, - double * const areas); + double * const areas, + double * const lat, + double * const lon); // // simulation setup diff --git a/src/executables/elm_ats_api/elm_ats_driver.cc b/src/executables/elm_ats_api/elm_ats_driver.cc index e579296f4..6f9592e54 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.cc +++ b/src/executables/elm_ats_api/elm_ats_driver.cc @@ -123,6 +123,9 @@ ELM_ATSDriver::parseParameterList() 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 }; @@ -225,6 +228,13 @@ ELM_ATSDriver::setup() 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. @@ -259,6 +269,24 @@ ELM_ATSDriver::MeshInfo ELM_ATSDriver::getMeshInfo() 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; } diff --git a/src/executables/elm_ats_api/elm_ats_driver.hh b/src/executables/elm_ats_api/elm_ats_driver.hh index 0ec5ce88c..eab10650d 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.hh +++ b/src/executables/elm_ats_api/elm_ats_driver.hh @@ -32,6 +32,8 @@ class ELM_ATSDriver : public Coordinator { int nlevgrnd; std::vector dzs; std::vector areas; + std::vector latitudes; + std::vector longitudes; }; ELM_ATSDriver(const Teuchos::RCP& plist, @@ -97,6 +99,8 @@ class ELM_ATSDriver : public Coordinator { Key surf_cv_key_; Key cv_key_; + Key lat_key_; + Key lon_key_; std::map key_map_; From 1f1267e130777ae46c4b5876fc9258db804e96d8 Mon Sep 17 00:00:00 2001 From: Ethan Coon Date: Thu, 11 Dec 2025 17:05:24 -0700 Subject: [PATCH 16/16] adds a simpler getmeshinfo that simply gets ncolumns, which is used to allocate ELM space --- src/executables/elm_ats_api/elm_ats_api.cc | 6 ++++++ src/executables/elm_ats_api/elm_ats_api.h | 8 +++++++- src/executables/elm_ats_api/elm_ats_driver.cc | 16 ++++++++-------- src/executables/elm_ats_api/elm_ats_driver.hh | 3 ++- 4 files changed, 23 insertions(+), 10 deletions(-) diff --git a/src/executables/elm_ats_api/elm_ats_api.cc b/src/executables/elm_ats_api/elm_ats_api.cc index 0a01f59d9..9d2c38139 100644 --- a/src/executables/elm_ats_api/elm_ats_api.cc +++ b/src/executables/elm_ats_api/elm_ats_api.cc @@ -61,6 +61,12 @@ void ats_parse_parameter_list(ELM_ATSDriver_ptr ats) void ats_get_mesh_info(ELM_ATSDriver_ptr ats, + int * const ncols_local) +{ + *ncols_local = reinterpret_cast(ats)->ncolumns; +} + +void ats_get_mesh_info2(ELM_ATSDriver_ptr ats, int * const ncols_local, int * const ncols_global, int * const nlevgrnd, diff --git a/src/executables/elm_ats_api/elm_ats_api.h b/src/executables/elm_ats_api/elm_ats_api.h index 0e9493d0a..84b2ae359 100644 --- a/src/executables/elm_ats_api/elm_ats_api.h +++ b/src/executables/elm_ats_api/elm_ats_api.h @@ -49,9 +49,15 @@ void ats_parse_parameter_list(ELM_ATSDriver_ptr ats); // -// These quantities should be compared against ELM to ensure consistent setup +// This sets the decomposition // void ats_get_mesh_info(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, diff --git a/src/executables/elm_ats_api/elm_ats_driver.cc b/src/executables/elm_ats_api/elm_ats_driver.cc index 6f9592e54..55f89ca6a 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.cc +++ b/src/executables/elm_ats_api/elm_ats_driver.cc @@ -65,7 +65,7 @@ ELM_ATSDriver::ELM_ATSDriver(const Teuchos::RCP& plist, int npfts) : Coordinator(plist, wallclock_timer, teuchos_comm, comm), npfts_(npfts), - ncolumns_(-1), + ncolumns(-1), ncells_per_col_(-1), elm_plist_(Teuchos::sublist(Teuchos::sublist(plist, "cycle driver"), "ELM driver")) { @@ -147,7 +147,7 @@ ELM_ATSDriver::parseParameterList() cv_key_ = Keys::getKey(domain_subsurf_, "cell_volume"); // -- mesh info - ncolumns_ = mesh_surf_->getNumEntities(AmanziMesh::CELL, AmanziMesh::Parallel_kind::OWNED); + ncolumns = mesh_surf_->getNumEntities(AmanziMesh::CELL, AmanziMesh::Parallel_kind::OWNED); auto col_zero = mesh_subsurf_->columns.getCells(0); ncells_per_col_ = col_zero.size(); @@ -252,7 +252,7 @@ ELM_ATSDriver::setup() ELM_ATSDriver::MeshInfo ELM_ATSDriver::getMeshInfo() { ELM_ATSDriver::MeshInfo info; - info.ncols_local = ncolumns_; + info.ncols_local = ncolumns; info.ncols_global = mesh_surf_->getMap(AmanziMesh::Entity_kind::CELL, false).NumGlobalElements(); info.nlevgrnd = ncells_per_col_; @@ -265,14 +265,14 @@ ELM_ATSDriver::MeshInfo ELM_ATSDriver::getMeshInfo() } // surface area - info.areas.resize(ncolumns_); - for (int i = 0; i != ncolumns_; ++i) { + 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.); + 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"); @@ -282,7 +282,7 @@ ELM_ATSDriver::MeshInfo ELM_ATSDriver::getMeshInfo() const Epetra_MultiVector& lon = *S_->Get(lon_key_, Tags::NEXT) .ViewComponent("cell", false); - for (int i = 0; i != ncolumns_; ++i) { + for (int i = 0; i != ncolumns; ++i) { info.latitudes[i] = lat[0][i]; info.longitudes[i] = lon[0][i]; } diff --git a/src/executables/elm_ats_api/elm_ats_driver.hh b/src/executables/elm_ats_api/elm_ats_driver.hh index eab10650d..72bb8ce21 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.hh +++ b/src/executables/elm_ats_api/elm_ats_driver.hh @@ -57,13 +57,14 @@ class ELM_ATSDriver : public Coordinator { double const * getFieldPtr(const ELM::VarID& var_id); double * getFieldPtrW(const ELM::VarID& var_id); + int ncolumns; + private: void initValue_(const Key& key, double value = 0.); private: Teuchos::RCP elm_plist_; - int ncolumns_; int ncells_per_col_; int npfts_;