From cf220acd49f01b6facc2e260f8faf0e08cfba2c5 Mon Sep 17 00:00:00 2001 From: Ethan Coon Date: Fri, 10 Jan 2025 10:09:04 -0700 Subject: [PATCH] fixes problems in previous merge of elm_api --- src/executables/elm_ats_api/CMakeLists.txt | 32 +- src/executables/elm_ats_api/elm_ats_api.cc | 41 +-- src/executables/elm_ats_api/elm_ats_api.h | 29 +- src/executables/elm_ats_api/elm_ats_driver.cc | 333 ++++++++++-------- src/executables/elm_ats_api/elm_ats_driver.hh | 18 +- .../elm_ats_api/test/CMakeLists.txt | 18 - .../elm_ats_api/test/elm_ats_test.cc | 59 ++-- .../elm_ats_api/test/fortran/CMakeLists.txt | 17 - .../elm_ats_api/test/fortran/ats.f90 | 9 +- .../elm_ats_api/test/fortran/ats_mod.f90 | 13 +- .../elm_ats_api/test/fortran/test.f90 | 32 +- 11 files changed, 288 insertions(+), 313 deletions(-) delete mode 100644 src/executables/elm_ats_api/test/CMakeLists.txt delete mode 100644 src/executables/elm_ats_api/test/fortran/CMakeLists.txt diff --git a/src/executables/elm_ats_api/CMakeLists.txt b/src/executables/elm_ats_api/CMakeLists.txt index 37db82422..ec7af2623 100644 --- a/src/executables/elm_ats_api/CMakeLists.txt +++ b/src/executables/elm_ats_api/CMakeLists.txt @@ -12,6 +12,7 @@ include(FortranCInterface) FortranCInterface_VERIFY(CXX) include_directories(${ATS_SOURCE_DIR}/src/executables) +include_directories(${ATS_SOURCE_DIR}/src/executables/elm_ats_api) set(elm_ats_src_files ${ATS_SOURCE_DIR}/src/executables/coordinator.cc @@ -36,12 +37,27 @@ if (APPLE AND BUILD_SHARED_LIBS) set_target_properties(elm_ats PROPERTIES LINK_FLAGS "-Wl,-undefined,dynamic_lookup") endif() -# ## build static lib -# add_amanzi_library(elm_ats_static -# STATIC -# SOURCE ${elm_ats_src_files} -# HEADERS ${elm_ats_inc_files} -# LINK_LIBS ${fates_link_libs} ${tpl_link_libs} ${ats_link_libs} ${amanzi_link_libs}) -# re-enable later -##add_subdirectory(test) +# +# These used to work, but need an inputfile that I'm not sure where it is... --ETC +# +# if (BUILD_TESTS) + +# # Add UnitTest includes +# include_directories(${UnitTest_INCLUDE_DIRS}) +# include_directories(${EXECUTABLE_SOURCE_DIR}) + +# # test for ELM API +# add_amanzi_test(elm_ats_cc_test elm_ats_cc_test +# KIND int +# SOURCE test/elm_ats_test.cc +# LINK_LIBS elm_ats ${ats_link_libs} ${UnitTest_LIBRARIES}) + +# # test for ELM API in fortran +# add_amanzi_test(elm_ats_fort_test elm_ats_fort_test +# KIND int +# SOURCE test/fortran/ats_mod.f90 test/fortran/test.f90 +# LINK_LIBS elm_ats ${ats_link_libs} ${UnitTest_LIBRARIES}) + +# endif() + diff --git a/src/executables/elm_ats_api/elm_ats_api.cc b/src/executables/elm_ats_api/elm_ats_api.cc index c2bebe3d7..d06b24b3d 100644 --- a/src/executables/elm_ats_api/elm_ats_api.cc +++ b/src/executables/elm_ats_api/elm_ats_api.cc @@ -69,11 +69,11 @@ void ats_setup(ELM_ATSDriver_ptr ats) // call driver initialize() void ats_initialize(ELM_ATSDriver_ptr ats, - double const * const t, - double const * const patm, - double const * const soilp) + double const * const start_time, + double const * const soil_water_content, + double const * const soil_pressure) { - reinterpret_cast(ats)->initialize(*t, patm, soilp); + reinterpret_cast(ats)->initialize(*start_time, soil_water_content, soil_pressure); } @@ -90,15 +90,6 @@ void ats_set_soil_hydrologic_parameters(ELM_ATSDriver_ptr ats, } -void ats_set_veg_parameters(ELM_ATSDriver_ptr ats, - double const * const mafic_potential_full_turgor, - double const * const mafic_potential_wilt_point) -{ - reinterpret_cast(ats) - ->set_veg_parameters(mafic_potential_full_turgor, mafic_potential_wilt_point); -} - - void ats_set_soil_hydrologic_properties(ELM_ATSDriver_ptr ats, double const * const effective_porosity) { @@ -118,38 +109,36 @@ void ats_set_veg_properties(ELM_ATSDriver_ptr ats, // 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) + double const * const surface_source, + double const * const potential_evaporation, + double const * const potential_transpiration) { reinterpret_cast(ats) - ->set_potential_sources(surface_infiltration, surface_evaporation, subsurface_transpiration); + ->set_potential_sources(surface_source, potential_evaporation, potential_transpiration); } void ats_get_waterstate(ELM_ATSDriver_ptr ats, - double * const surface_ponded_depth, + double * const ponded_depth, double * const water_table_depth, - double * const soil_pressure, - double * const soil_psi, - double * const sat_liq, - double * const sat_ice) + double * const mass_water_content) { reinterpret_cast(ats) - ->get_waterstate(surface_ponded_depth, water_table_depth, soil_pressure, soil_psi, sat_liq, sat_ice); + ->get_waterstate(ponded_depth, water_table_depth, mass_water_content); } void ats_get_water_fluxes(ELM_ATSDriver_ptr ats, - double * const soil_infiltration, + double * const infiltration, double * const evaporation, double * const transpiration, + double * const root_flux, double * net_subsurface_fluxes, double * net_runon) { reinterpret_cast(ats) - ->get_water_fluxes(soil_infiltration, evaporation, transpiration, - net_subsurface_fluxes, net_runon); + ->get_water_fluxes(infiltration, evaporation, transpiration, + root_flux, 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..a4dde5b48 100644 --- a/src/executables/elm_ats_api/elm_ats_api.h +++ b/src/executables/elm_ats_api/elm_ats_api.h @@ -67,9 +67,9 @@ void ats_setup(ELM_ATSDriver_ptr ats); // call driver initialize() void ats_initialize(ELM_ATSDriver_ptr ats, - double const * const t, - double const * const patm, - double const * const soilp); + double const * const start_time, + double const * const soil_water_content, + double const * const soil_pressure); // set material parameters, which are constant in time void ats_set_soil_hydrologic_parameters(ELM_ATSDriver_ptr ats, @@ -79,11 +79,6 @@ void ats_set_soil_hydrologic_parameters(ELM_ATSDriver_ptr ats, double const * const clapp_horn_smpsat, double const * const clapp_horn_sr); -// set veg parameters, which are constant in time -void ats_set_veg_parameters(ELM_ATSDriver_ptr ats, - double const * const mafic_potential_full_turgor, - double const * const mafic_potential_wilt_point); - // // Called prior to timestep advance // ----------------------------------------------------------------------------- @@ -96,12 +91,10 @@ 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); + double const * const surface_source, + double const * const potential_evaporation, + double const * const potential_transpiration); // @@ -111,17 +104,15 @@ void ats_set_sources(ELM_ATSDriver_ptr ats, // 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 ponded_depth, double * const water_table_depth, - double * const soil_pressure, - double * const soil_psi, - double * const sat_liq, - double * const sat_ice); + double * const mass_water_content); void ats_get_water_fluxes(ELM_ATSDriver_ptr ats, - double * const soil_infiltration, + double * const infiltration, double * const evaporation, double * const transpiration, + double * const root_flux, double * net_subsurface_fluxes, double * net_runon); diff --git a/src/executables/elm_ats_api/elm_ats_driver.cc b/src/executables/elm_ats_api/elm_ats_driver.cc index 137850bb6..b909e5413 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.cc +++ b/src/executables/elm_ats_api/elm_ats_driver.cc @@ -2,6 +2,7 @@ #include #include #include +#include #include "errors.hh" #include "dbc.hh" @@ -59,6 +60,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, @@ -109,7 +111,7 @@ ELM_ATSDriver::ELM_ATSDriver(const Teuchos::RCP& plist, 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"); + ss_flux_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"); @@ -229,7 +231,7 @@ ELM_ATSDriver::setup() requireAtNext(trans_key_, Amanzi::Tags::NEXT, *S_) .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1); - requireAtNext(infilt_key_, Amanzi::Tags::NEXT, *S_) + requireAtNext(pot_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); @@ -239,7 +241,6 @@ ELM_ATSDriver::setup() } - // // dz & depth are currently ignored -- presumed identical between ATS & ELM // @@ -296,15 +297,28 @@ void ELM_ATSDriver::initialize(double t, // initialize ATS data, commit initial conditions Coordinator::initialize(); - // initialize pressure field + // initialize pressure via + // constant depth from surface water table, or + // from ELM's water content + // TODO connect this to an input parameter ELM_ATSDriver::init_pressure_from_wc_(elm_water_content); + //ELM_ATSDriver::init_pressure_from_wt_(1.0); + + // 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_); // set as zero and tag as initialized initZero_(root_frac_key_); initZero_(pot_infilt_key_); initZero_(pot_trans_key_); initZero_(pot_evap_key_); - initZero_(infilt_key_); + initZero_(ss_flux_key_); initZero_(trans_key_); initZero_(evap_key_); initZero_(total_trans_key_); @@ -314,79 +328,6 @@ void ELM_ATSDriver::initialize(double t, 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) { @@ -484,13 +425,6 @@ void ELM_ATSDriver::set_soil_hydrologic_parameters(double const * const base_por } -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 @@ -504,35 +438,37 @@ void ELM_ATSDriver::set_veg_properties(double const * const rooting_fraction) void -ELM_ATSDriver::set_potential_sources(double const * const elm_surface_input, +ELM_ATSDriver::set_potential_sources(double const * const elm_surface_source, double const * const elm_evaporation, double const * const elm_transpiration) { - // ELM's infiltration, evaporation, and transpiration are in units of mm s^-1 + // infiltration, evaporation, and transpiration + // ELM's fluxes are in units of mm s^-1 (or kg m^-2 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_) + + auto& input = *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_) + auto& evap = *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_) + auto& trans = *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_); + AMANZI_ASSERT(input.MyLength() == ncolumns_); + AMANZI_ASSERT(evap.MyLength() == ncolumns_); + AMANZI_ASSERT(trans.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); + const double m_per_mm{0.001}; 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; + const double factor = m_per_mm * surf_d[0][i]; + input[0][i] = elm_surface_source[i] * factor; + evap[0][i] = elm_evaporation[i] * factor; + trans[0][i] = elm_transpiration[i] * factor; } changedEvaluatorPrimary(pot_infilt_key_, Amanzi::Tags::NEXT, *S_); @@ -544,10 +480,7 @@ ELM_ATSDriver::set_potential_sources(double const * const elm_surface_input, 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? + double * const sat_liq) { // convert saturation into [kg/m2] from [-] S_->GetEvaluator(sat_key_, Amanzi::Tags::NEXT).Update(*S_, sat_key_); @@ -562,14 +495,14 @@ ELM_ATSDriver::get_waterstate(double * const ponded_depth, 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 + const auto dz = calcDZ_(); + 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; + sat_liq[j * ncolumns_ + i] = satl[0][cells[j]] * por[0][cells[j]] * dens[0][cells[j]] * dz[j]; } } @@ -584,28 +517,6 @@ ELM_ATSDriver::get_waterstate(double * const ponded_depth, // 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]]; -// } -// } } @@ -613,62 +524,190 @@ void ELM_ATSDriver::get_water_fluxes(double * const surf_subsurf_flx, double * const evaporation, double * const transpiration, + double * const root_flux, 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] + // transpiration [mol / m2 / s] [mmH2o/s] + // root_flux [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) + S_->GetEvaluator(ss_flux_key_, Amanzi::Tags::NEXT).Update(*S_, ss_flux_key_); + const auto& ss_flux = *S_->Get(ss_flux_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); + // Surface transpiration fluxes - at the leaves + S_->GetEvaluator(total_trans_key_, Amanzi::Tags::NEXT).Update(*S_, total_trans_key_); + const auto& trans = *S_->Get(total_trans_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 + // 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; + surf_subsurf_flx[i] = ss_flux[0][i] * mm_per_mol; evaporation[i] = evap[0][i] * mm_per_mol; + transpiration[i] = trans[0][i] * mm_per_mol; } - // Subsurface flux + // Subsurface transpiration fluxes - at the roots S_->GetEvaluator(trans_key_, Amanzi::Tags::NEXT).Update(*S_, "ELM"); - const auto& trans = *S_->Get(trans_key_, Amanzi::Tags::NEXT) + const auto& rootflux = *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); + const auto dz = calcDZ_(); - // convert mol/m3/s to mmH2O/s by integrating over dz - NO? - // treat the same as surface fluxes? + // convert mol/m3/s to mmH2O/s for ELM + // mol/m3/s * m3/mol * m * mm/m = mm/s for (int i=0; i!=ncolumns_; ++i) { - const auto& faces = mesh_subsurf_->columns.getFaces(i); - const auto& cells = mesh_subsurf_->columns.getCells(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; + double factor = dz[j] * 1000.0 / surfdens[0][i]; // use surface density at roots and leaves - test later! + root_flux[j * ncolumns_ + i] = rootflux[0][cells[j]] * factor; } } // unclear how to implement net_subsurface_fluxes and net_runon... --ETC } +// use prescribed WT depth to initialize terrain-following hydrostatic pressure field +void ELM_ATSDriver::init_pressure_from_wt_(double depth_to_wt) +{ + // atmospheric pressure and density-gravity factor + // hardwired for now + const double p_atm = 101325.0; + const double rho_g = 9806.65; + + // writable to pressure + auto& pres = *S_->GetW(pres_key_, Amanzi::Tags::NEXT, "flow") + .ViewComponent("cell", false); + + const auto cell_centers = calcCellDepths_(); + + // apply hydrostatic pressure on each column + for (int i=0; i!=ncolumns_; ++i) { + const auto cells = mesh_subsurf_->columns.getCells(i); + for (int j=0; j!=ncells_per_col_; ++j) { + pres[0][cells[j]] = p_atm + rho_g * (cell_centers[j] - depth_to_wt); + } + } +} + +// use incoming water content to initialize pressure field +void ELM_ATSDriver::init_pressure_from_wc_(double const * const elm_water_content) +{ + // atmospheric pressure and density-gravity factor + // hardwired for now + const double p_atm = 101325.0; + const double rho_g = 9806.65; + + // 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); + + // 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]; + + const auto dz = calcDZ_(); + + // 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 = 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 factor = 1 / (dz[j] * por[0][cells[j]] * mass_d[0][cells[j]]); + double satl = elm_water_content[j * ncolumns_ + i] * factor; + if (j > 9) satl = 1.0; // hardwire bottom 5 layers (ELM's bedrock) as fully saturated + if (satl < 1.0) { + pres[0][cells[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[j]; + pres[0][cells[j]] = p_atm + rho_g * (sat_depth - dz[j]/2); + } + } + } +} + + +// calculate dz for each cell in column +// only use this if all columns have identical vertical spacing +// returns a vector of length ncells_per_column +std::vector ELM_ATSDriver::calcDZ_() +{ + // use volume/area to calculate dz + 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); + + // assume all columns have the same dz spacing + const auto cells = mesh_subsurf_->columns.getCells(0); + std::vector dz; + for (int i=0; i ELM_ATSDriver::calcCellDepths_() +{ + const auto dz = calcDZ_(); // dz from cell volume/area + + std::vector cell_bot_depth(dz.size()); + std::partial_sum(dz.begin(), dz.end(), cell_bot_depth.begin()); // depth of bottom face + + std::vector cell_centers; + std::transform( + dz.begin(), dz.end(), + cell_bot_depth.begin(), + std::back_inserter(cell_centers), + [](const auto& delz, const auto& depth) + { return depth - delz/2; } + ); + + return cell_centers; // depth to cell center from surface +} + void ELM_ATSDriver::initZero_(const Key& key) { auto& vec = S_->GetW(key, Amanzi::Tags::NEXT, key); @@ -702,7 +741,7 @@ void ELM_ATSDriver::copyToSub_(double const * const in, const Key& key, Key owne .ViewComponent("cell", false); for (int i=0; i!=ncolumns_; ++i) { - const auto& cells = mesh_subsurf_->columns.getCells(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]; } diff --git a/src/executables/elm_ats_api/elm_ats_driver.hh b/src/executables/elm_ats_api/elm_ats_driver.hh index 5b76a564a..49cdc9e0e 100644 --- a/src/executables/elm_ats_api/elm_ats_driver.hh +++ b/src/executables/elm_ats_api/elm_ats_driver.hh @@ -81,29 +81,29 @@ public: 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_waterstate(double * const ponded_depth, + double * const wt_depth, + double * const sat_liq); void get_water_fluxes(double * const soil_infiltration, double * const evaporation, double * const transpiration, + double * const root_flux, double * const net_subsurface_fluxes, double * const net_runon); private: void init_pressure_from_wc_(double const * const elm_water_content); + void init_pressure_from_wt_(double depth_to_wt); + std::vector calcDZ_(); + std::vector calcCellDepths_(); void copyToSurf_(double const * const in, const Key& key, Key owner=""); void copyToSub_(double const * const in, const Key& key, Key owner=""); @@ -146,7 +146,7 @@ public: Key sat_key_; //Key sat_gas_key_; //Key sat_ice_key_; - Key infilt_key_; + Key ss_flux_key_; Key trans_key_; Key evap_key_; Key total_trans_key_; diff --git a/src/executables/elm_ats_api/test/CMakeLists.txt b/src/executables/elm_ats_api/test/CMakeLists.txt deleted file mode 100644 index 02b3b84de..000000000 --- a/src/executables/elm_ats_api/test/CMakeLists.txt +++ /dev/null @@ -1,18 +0,0 @@ -# -*- mode: cmake -*- -# -# ELM_ATS_API -# example/test drivers -# -# build both tests for now -# need to make configurable - -include_directories(${ATS_SOURCE_DIR}/src/executables) -include_directories(${ATS_SOURCE_DIR}/src/executables/elm_ats_api) - -add_amanzi_executable(elm_ats_cc_test - SOURCE elm_ats_test.cc - LINK_LIBS elm_ats ${fates_link_libs} ${tpl_link_libs} ${ats_link_libs} ${amanzi_link_libs} - OUTPUT_NAME elm_ats_cc_test - OUTPUT_DIRECTORY ${ATS_BINARY_DIR}) - -add_subdirectory(fortran) diff --git a/src/executables/elm_ats_api/test/elm_ats_test.cc b/src/executables/elm_ats_api/test/elm_ats_test.cc index 776591149..be0b56308 100644 --- a/src/executables/elm_ats_api/test/elm_ats_test.cc +++ b/src/executables/elm_ats_api/test/elm_ats_test.cc @@ -1,6 +1,5 @@ #include -#include #include #include @@ -75,53 +74,35 @@ int main(int argc, char *argv[]) if (input_filename.empty() && !opt_input_filename.empty()) input_filename = opt_input_filename; // dummy data - // 1 col, 100 cells double time = 0.; - int n = 1; - int m = 100; + int ncols = 5; + int ncells = 15; int ncols_local, ncols_global, ncells_per_col; - std::vector soil_infil(n, 0.0); - std::vector soil_evap(n, 0.0); - std::vector air_pres(n); - std::vector surf_pres(n); - std::vector elev(n); - std::vector surf_area_m2(n); - std::vector lat(n); - std::vector lon(n); - - std::vector dz(m); - std::vector depth(m); - std::vector root_tran(m); - std::vector soil_pres(m); - std::vector soil_pot(m); - std::vector satl(m); - std::vector sati(m); - std::vector pft_i(m); + std::vector soil_infil(ncols, 10.0); + std::vector soil_evap(ncols, 3.0); + std::vector root_tran(ncols, 6.0); + std::vector soil_pres(ncols*ncells); + std::vector satl(ncols*ncells); // dummy fortran comm MPI_Fint comm = 0; // test driver directly - auto driver = std::unique_ptr(ATS::createELM_ATSDriver(&comm, input_filename.data())); - driver->setup(); - //driver->get_mesh_info(ncols_local, ncols_global, lat.data(), lon.data(), elev.data(), surf_area_m2.data(), pft_i.data(), ncells_per_col, depth.data()); - driver->initialize(); - //driver->set_potential_sources(soil_infil.data(), soil_evap.data(), root_tran.data()); - driver->advance_test(); - //driver->get_waterstate(surf_pres.data(), soil_pres.data(), soil_pot.data(), satl.data(), sati.data()); - driver->finalize(); + auto driver = std::unique_ptr(ATS::createELM_ATSDriver(&comm, input_filename.data())); + driver->setup(); + driver->initialize(time, soil_pres.data(), satl.data()); + driver->advance_test(); + driver->finalize(); + std::cout << "DONE WITH DRIVER TEST" << std::endl; // test api - //auto driver_api = ats_create_c(&comm, input_filename.data()); - //ats_setup_c(driver_api); - //ats_get_mesh_info_c(driver_api, &ncols_local, &ncols_global, lat.data(), lon.data(), - // elev.data(), surf_area_m2.data(), pft_i.data(), &ncells_per_col, depth.data()); - //ats_initialize_c(driver_api, &time, air_pres.data(), soil_pres.data()); - //ats_set_potential_sources_c(driver_api, soil_infil.data(), soil_evap.data(), root_tran.data()); - //ats_advance_test_c(driver_api); - //ats_get_waterstate_c(driver_api, surf_pres.data(), soil_pres.data(), soil_pot.data(), satl.data(), sati.data()); - //ats_delete_c(driver_api); - std::cout << "DONE WITH ELM-ATS C++ TEST" << std::endl; + + auto driver_api = ats_create(&comm, input_filename.data()); + ats_setup(driver_api); + ats_initialize(driver_api, &time, soil_pres.data(), satl.data()); + ats_advance_test(driver_api); + ats_delete(driver_api); + std::cout << "DONE WITH API TEST" << std::endl; return 0; } diff --git a/src/executables/elm_ats_api/test/fortran/CMakeLists.txt b/src/executables/elm_ats_api/test/fortran/CMakeLists.txt deleted file mode 100644 index 1dedda87a..000000000 --- a/src/executables/elm_ats_api/test/fortran/CMakeLists.txt +++ /dev/null @@ -1,17 +0,0 @@ -# -*- mode: cmake -*- -# -# ELM_ATS_API -# example/test Fortran driver -# - -set(fort_src_files - ats_mod.f90 - test.f90 - ) - -add_amanzi_executable(fort_driver - SOURCE ${fort_src_files} - LINK_LIBS elm_ats - OUTPUT_NAME elm_ats_f90_test - OUTPUT_DIRECTORY ${ATS_BINARY_DIR}) - diff --git a/src/executables/elm_ats_api/test/fortran/ats.f90 b/src/executables/elm_ats_api/test/fortran/ats.f90 index 85891b5d8..90fd330f9 100644 --- a/src/executables/elm_ats_api/test/fortran/ats.f90 +++ b/src/executables/elm_ats_api/test/fortran/ats.f90 @@ -24,10 +24,13 @@ subroutine ats_setup_c(ats) bind(c, name="ats_setup") type(c_ptr), value :: ats end subroutine - subroutine ats_initialize_c(ats) bind(c, name="ats_initialize") + subroutine ats_initialize_c(ats, time, satur, soil_pres) bind(c, name="ats_initialize") use iso_c_binding implicit none type(c_ptr), value :: ats + real(c_double), intent(in) :: time + real(c_double), dimension(*), intent(in) :: satur + real(c_double), dimension(*), intent(in) :: soil_pres end subroutine subroutine ats_advance_test_c(ats) bind(c, name="ats_advance_test") @@ -43,7 +46,7 @@ subroutine ats_advance_c(ats, dt) bind(c, name="ats_advance") real(c_double), intent(in) :: dt end subroutine - subroutine ats_set_sources_c(ats, soil_infil, soil_evap, root_trans, ncols, ncells) & + subroutine ats_set_sources_c(ats, soil_infil, soil_evap, root_trans) & bind(c, name="ats_set_sources") use iso_c_binding implicit none @@ -51,8 +54,6 @@ subroutine ats_set_sources_c(ats, soil_infil, soil_evap, root_trans, ncols, ncel real(c_double), dimension(*), intent(in) :: soil_infil real(c_double), dimension(*), intent(in) :: soil_evap real(c_double), dimension(*), intent(in) :: root_trans - integer(c_int), intent(in) :: ncols - integer(c_int), intent(in) :: ncells end subroutine subroutine ats_get_waterstate_c(ats, surf_pres, soil_pres, satur, ncols, ncells) & diff --git a/src/executables/elm_ats_api/test/fortran/ats_mod.f90 b/src/executables/elm_ats_api/test/fortran/ats_mod.f90 index 01abed891..2b612fdc1 100644 --- a/src/executables/elm_ats_api/test/fortran/ats_mod.f90 +++ b/src/executables/elm_ats_api/test/fortran/ats_mod.f90 @@ -65,10 +65,13 @@ subroutine ats_setup(this) end subroutine !! Initialize ATS for this simulation - subroutine ats_initialize(this) + subroutine ats_initialize(this, time, satur, soil_pres) implicit none class(ats) :: this - call ats_initialize_c(this%ptr) + double precision, intent(in) :: time + double precision, dimension(*), intent(in) :: satur + double precision, dimension(*), intent(in) :: soil_pres + call ats_initialize_c(this%ptr, time, satur, soil_pres) end subroutine !! Advance ATS by dt @@ -88,15 +91,13 @@ subroutine ats_advance_test(this) !! WIP !! Sets ATS coupling variables infiltration, evaporation, and transpiration - subroutine ats_set_sources(this, soil_infil, soil_evap, root_trans, ncols, ncells) + subroutine ats_set_sources(this, soil_infil, soil_evap, root_trans) implicit none class(ats) :: this double precision, dimension(*), intent(in) :: soil_infil double precision, dimension(*), intent(in) :: soil_evap double precision, dimension(*), intent(in) :: root_trans - integer, intent(in) :: ncols - integer, intent(in) :: ncells - call ats_set_sources_c(this%ptr, soil_infil, soil_evap, root_trans, ncols, ncells) + call ats_set_sources_c(this%ptr, soil_infil, soil_evap, root_trans) end subroutine !! WIP diff --git a/src/executables/elm_ats_api/test/fortran/test.f90 b/src/executables/elm_ats_api/test/fortran/test.f90 index 3bfcadff3..91b492ee3 100644 --- a/src/executables/elm_ats_api/test/fortran/test.f90 +++ b/src/executables/elm_ats_api/test/fortran/test.f90 @@ -8,30 +8,22 @@ program elm_test integer :: ierror, i ! dummy data - ! 1 column, 100 cells - integer, parameter :: ncol = 1 - integer, parameter :: ncell = 100 + integer, parameter :: ncol = 5 + integer, parameter :: ncell = 15 + double precision, dimension(ncol*ncell) :: soil_pres + double precision, dimension(ncol*ncell) :: satur + double precision, dimension(ncol) :: tran double precision, dimension(ncol) :: infil double precision, dimension(ncol) :: evap - double precision, dimension(ncol) :: surf_pres - double precision, dimension(ncol) :: elev - double precision, dimension(ncol) :: surf_area_m2 - double precision, dimension(ncol) :: lat - double precision, dimension(ncol) :: lon - - double precision, dimension(ncell) :: dz - double precision, dimension(ncell) :: depth - double precision, dimension(ncell) :: soil_pres - double precision, dimension(ncell) :: satur - double precision, dimension(ncell) :: tran - integer :: ncols_local, ncols_global, ncells_per_col + double precision :: time character (len = 40) :: test_prefix infil(:) = 10.0 evap(:) = 3.0 + tran(:) = 6.0 - do i=1,ncell + do i=1,ncol tran(i) = (1.0 - (1.0/ncell)*i) end do @@ -42,6 +34,8 @@ program elm_test infile_name = trim(adjustl(infile_name(12:))) end if + time = 0.0 + call MPI_INIT(ierror) ! create an ATS driver object @@ -49,11 +43,9 @@ program elm_test ! call ATS methods call ats_driver%setup() - !!call ats_driver%get_mesh_info(ncols_local, ncols_global, ncells_per_col, dz, depth, elev, surf_area_m2, lat, lon) - call ats_driver%initialize() - !!call ats_driver%set_sources(infil, evap, tran, ncol, ncell) + call ats_driver%initialize(time, satur, soil_pres) + !call ats_driver%set_sources(infil, evap, tran) call ats_driver%advance_test() - !!call ats_driver%get_waterstate(surf_pres, soil_pres, satur, ncol, ncell) ! don't need to call ats_driver%delete now that final is used