From 091e979b2019d3506ce70eb20d32e7dabf27cba9 Mon Sep 17 00:00:00 2001 From: Joe Beisman Date: Mon, 3 Oct 2022 04:20:43 -0400 Subject: [PATCH 1/7] basic skeleton PK for ELMKernels coupling --- src/pks/surface_balance/CMakeLists.txt | 19 + .../ELMKernels/surface_balance_ELMKernels.cc | 452 ++++++++++++++++++ .../ELMKernels/surface_balance_ELMKernels.hh | 107 +++++ .../surface_balance_ELMKernels_reg.hh | 19 + 4 files changed, 597 insertions(+) create mode 100644 src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc create mode 100644 src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.hh create mode 100644 src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels_reg.hh diff --git a/src/pks/surface_balance/CMakeLists.txt b/src/pks/surface_balance/CMakeLists.txt index e8e1c1902..c16cd59dd 100644 --- a/src/pks/surface_balance/CMakeLists.txt +++ b/src/pks/surface_balance/CMakeLists.txt @@ -51,6 +51,12 @@ if (ENABLE_CLM) ) endif() +if (ENABLE_ELMKERNELS) + list(APPEND ats_surface_balance_src_files + ELMKernels/surface_balance_ELMKernels.cc + ) +endif() + set(ats_surface_balance_inc_files constitutive_relations/land_cover/seb_physics_defs.hh constitutive_relations/land_cover/seb_physics_funcs.hh @@ -95,6 +101,12 @@ if (ENABLE_CLM) ) endif() +if (ENABLE_ELMKERNELS) + list(APPEND ats_surface_balance_inc_files + ELMKernels/surface_balance_ELMKernels.hh + ) +endif() + set(ats_surface_balance_link_libs ${Teuchos_LIBRARIES} ${Epetra_LIBRARIES} @@ -145,6 +157,13 @@ if (ENABLE_CLM) ) endif() +if (ENABLE_ELMKERNELS) + register_evaluator_with_factory( + HEADERFILE ELMKernels/surface_balance_ELMKernels_reg.hh + LISTNAME ATS_SURFACE_BALANCE_REG + ) +endif() + register_evaluator_with_factory( HEADERFILE surface_balance_base_reg.hh LISTNAME ATS_SURFACE_BALANCE_REG diff --git a/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc new file mode 100644 index 000000000..7e9130fbc --- /dev/null +++ b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc @@ -0,0 +1,452 @@ +/* + 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. + +*/ + +#include + +#include "pk_helpers.hh" +#include "surface_balance_ELMKernels.hh" + +namespace Amanzi { +namespace SurfaceBalance { + +SurfaceBalanceELMKernels::SurfaceBalanceELMKernels(Teuchos::ParameterList& pk_tree, + const Teuchos::RCP& global_list, + const Teuchos::RCP& S, + const Teuchos::RCP& solution): + PK(pk_tree, global_list, S, solution), + PK_Physical_Default(pk_tree, global_list, S, solution) +{ + domain_ss_ = Keys::readDomainHint(*plist_, domain_, "surface", "subsurface"); + domain_snow_ = Keys::readDomainHint(*plist_, domain_, "surface", "snow"); + domain_can_ = Keys::readDomainHint(*plist_, domain_, "surface", "canopy"); + + // primary variables + key_ = Keys::readKey(*plist_, domain_snow_, "snow depth", key_); + surf_water_src_key_ = Keys::readKey(*plist_, domain_, "surface water source", "water_source"); + ss_water_src_key_ = Keys::readKey(*plist_, domain_ss_, "subsurface water source", "water_source"); + + // diagnostic keys + // evap_flux_key_ = Keys::readKey(*plist_, domain_, "evaporative flux", "evaporative_flux"); + qE_lh_key_ = Keys::readKey(*plist_, domain_, "latent heat of evaporation", "qE_latent_heat"); + qE_sh_key_ = Keys::readKey(*plist_, domain_, "sensible heat flux", "qE_sensible_heat"); + qE_lw_out_key_ = Keys::readKey(*plist_, domain_, "outgoing longwave radiation", "qE_lw_out"); + qE_cond_key_ = Keys::readKey(*plist_, domain_, "conducted energy flux", "qE_conducted"); + + snow_swe_key_ = Keys::readKey(*plist_, domain_snow_, "snow water equivalent", "water_equivalent"); + can_wc_key_ = Keys::readKey(*plist_, domain_can_, "canopy water content", "water_content"); + surf_temp_key_ = Keys::readKey(*plist_, domain_, "surface temperature", "temperature"); + soil_temp_key_ = Keys::readKey(*plist_, domain_ss_, "soil temperature", "temperature"); + can_temp_key_ = Keys::readKey(*plist_, domain_can_, "canopy temperature", "temperature"); + + // dependencies + // -- met data + met_sw_key_ = Keys::readKey(*plist_, domain_,"incoming shortwave radiation", "incoming_shortwave_radiation"); + met_lw_key_ = Keys::readKey(*plist_, domain_,"incoming longwave radiation", "incoming_longwave_radiation"); + met_air_temp_key_ = Keys::readKey(*plist_, domain_,"air temperature", "air_temperature"); + met_rel_hum_key_ = Keys::readKey(*plist_, domain_,"relative humidity", "relative_humidity"); + met_wind_speed_key_ = Keys::readKey(*plist_, domain_,"wind speed", "wind_speed"); + met_prain_key_ = Keys::readKey(*plist_, domain_,"precipitation rain", "precipitation_rain"); + met_psnow_key_ = Keys::readKey(*plist_, domain_snow_,"precipitation snow", "precipitation"); + + // soil state + pres_key_ = Keys::readKey(*plist_, domain_ss_, "pressure", "pressure"); + sl_key_ = Keys::readKey(*plist_, domain_ss_, "saturation_liquid", "saturation_liquid"); + + // soil properties + sand_frac_key_ = Keys::readKey(*plist_, domain_ss_, "sand fraction", "sand_fraction"); + silt_frac_key_ = Keys::readKey(*plist_, domain_ss_, "silt fraction", "silt_fraction"); + clay_frac_key_ = Keys::readKey(*plist_, domain_ss_, "clay fraction", "clay_fraction"); + poro_key_ = Keys::readKey(*plist_, domain_ss_, "porosity", "porosity"); + + // surface properties + color_index_key_ = Keys::readKey(*plist_, domain_ss_, "color index", "color_index"); + pft_index_key_ = Keys::readKey(*plist_, domain_ss_, "PFT index", "pft_index"); + + // ELMKernels timestep + dt_ = plist_->get("time step size [s]"); +} + +// main methods +// -- Setup data. +void +SurfaceBalanceELMKernels::Setup() +{ + PK_Physical_Default::Setup(); + auto subsurf_mesh = S_->GetMesh(domain_ss_); + + // requirements: primary variable + // -- snow depth + S_->Require(key_, tag_next_, name_) + .SetMesh(mesh_)->SetComponent("cell", AmanziMesh::CELL, 1); + + // requirements: other primary variables + // -- surface water source -- note we keep old and new in case of Crank-Nicholson Richards PK + S_->Require(surf_water_src_key_, tag_next_, name_) + .SetMesh(mesh_)->SetComponent("cell", AmanziMesh::CELL, 1); + requireEvaluatorPrimary(surf_water_src_key_, tag_next_, *S_); + + // -- subsurface water source -- note we keep old and new in case of Crank-Nicholson Richards PK + S_->Require(ss_water_src_key_, tag_next_, name_) + .SetMesh(subsurf_mesh)->SetComponent("cell", AmanziMesh::CELL, 1); + requireEvaluatorPrimary(ss_water_src_key_, tag_next_, *S_); + + // set requirements on dependencies + SetupDependencies_(tag_next_); + + // Set up the ELMKernels object + //ATS::ELMKernels::init(subsurf_mesh->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED), + // mesh_->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED), + // 2, mesh_->get_comm()->MyPID(), 3); +} + +void +SurfaceBalanceELMKernels::SetupDependencies_(const Tag& tag) +{ + auto subsurf_mesh = S_->GetMesh(domain_ss_); + auto snow_mesh = S_->GetMesh(domain_snow_); + auto can_mesh = S_->GetMesh(domain_can_); + + // requirements: energy balance diagnostic variables. Only at the new time. + // No evaluators for now? + // S_->Require(evap_flux_key_, tag_next_, name_) + // .SetMesh(mesh_)->SetComponent("cell", AmanziMesh::CELL, 1); + // S_->GetRecord(evap_flux_key_, tag_next_).set_io_checkpoint(false); + + S_->Require(qE_lh_key_, tag_next_, name_) + .SetMesh(mesh_)->SetComponent("cell", AmanziMesh::CELL, 1); + S_->GetRecordW(qE_lh_key_, tag_next_, name_).set_io_checkpoint(false); + + S_->Require(qE_sh_key_, tag_next_, name_) + .SetMesh(mesh_)->SetComponent("cell", AmanziMesh::CELL, 1); + S_->GetRecordW(qE_sh_key_, tag_next_, name_).set_io_checkpoint(false); + + S_->Require(qE_lw_out_key_, tag_next_, name_) + .SetMesh(mesh_)->SetComponent("cell", AmanziMesh::CELL, 1); + S_->GetRecordW(qE_lw_out_key_, tag_next_, name_).set_io_checkpoint(false); + + S_->Require(qE_cond_key_, tag_next_, name_) + .SetMesh(mesh_)->SetComponent("cell", AmanziMesh::CELL, 1); + S_->GetRecordW(qE_cond_key_, tag_next_, name_).set_io_checkpoint(false); + + // requirements: other diagnostics + S_->Require(snow_swe_key_, tag_next_, name_) + .SetMesh(snow_mesh)->SetComponent("cell", AmanziMesh::CELL, 1); + S_->GetRecordW(snow_swe_key_, tag_next_, name_).set_io_checkpoint(false); + + S_->Require(can_wc_key_, tag_next_, name_) + .SetMesh(can_mesh)->SetComponent("cell", AmanziMesh::CELL, 1); + S_->GetRecordW(can_wc_key_, tag_next_, name_).set_io_checkpoint(false); + + S_->Require(surf_temp_key_, tag_next_, name_) + .SetMesh(mesh_)->SetComponent("cell", AmanziMesh::CELL, 1); + S_->GetRecordW(surf_temp_key_, tag_next_, name_).set_io_checkpoint(false); + + S_->Require(soil_temp_key_, tag_next_, name_) + .SetMesh(subsurf_mesh)->SetComponent("cell", AmanziMesh::CELL, 1); + S_->GetRecordW(soil_temp_key_, tag_next_, name_).set_io_checkpoint(false); + + S_->Require(can_temp_key_, tag_next_, name_) + .SetMesh(can_mesh)->SetComponent("cell", AmanziMesh::CELL, 1); + S_->GetRecordW(can_temp_key_, tag_next_, name_).set_io_checkpoint(false); + + // requirements: independent variables (data from MET) + S_->RequireEvaluator(met_sw_key_, tag); + S_->Require(met_sw_key_, tag) + .SetMesh(mesh_)->AddComponent("cell", AmanziMesh::CELL, 1); + + S_->RequireEvaluator(met_lw_key_, tag); + S_->Require(met_lw_key_, tag) + .SetMesh(mesh_)->AddComponent("cell", AmanziMesh::CELL, 1); + + S_->RequireEvaluator(met_air_temp_key_, tag); + S_->Require(met_air_temp_key_, tag) + .SetMesh(mesh_)->AddComponent("cell", AmanziMesh::CELL, 1); + + S_->RequireEvaluator(met_rel_hum_key_, tag); + S_->Require(met_rel_hum_key_, tag) + .SetMesh(mesh_)->AddComponent("cell", AmanziMesh::CELL, 1); + + S_->RequireEvaluator(met_wind_speed_key_, tag); + S_->Require(met_wind_speed_key_, tag) + .SetMesh(mesh_)->AddComponent("cell", AmanziMesh::CELL, 1); + + S_->RequireEvaluator(met_prain_key_, tag); + S_->Require(met_prain_key_, tag) + .SetMesh(mesh_)->AddComponent("cell", AmanziMesh::CELL, 1); + + S_->RequireEvaluator(met_psnow_key_, tag); + S_->Require(met_psnow_key_, tag) + .SetMesh(snow_mesh)->AddComponent("cell", AmanziMesh::CELL, 1); + + // requirements: soil state + S_->RequireEvaluator(pres_key_, tag); + S_->Require(pres_key_, tag) + .SetMesh(subsurf_mesh)->AddComponent("cell", AmanziMesh::CELL, 1); + S_->RequireEvaluator(sl_key_, tag); + S_->Require(sl_key_, tag) + .SetMesh(subsurf_mesh)->AddComponent("cell", AmanziMesh::CELL, 1); + + // requirements: soil properties + S_->RequireEvaluator(poro_key_, tag); + S_->Require(poro_key_, tag) + .SetMesh(subsurf_mesh)->AddComponent("cell", AmanziMesh::CELL, 1); + S_->RequireEvaluator(sand_frac_key_, tag); + S_->Require(sand_frac_key_, tag) + .SetMesh(subsurf_mesh)->AddComponent("cell", AmanziMesh::CELL, 1); + S_->RequireEvaluator(silt_frac_key_, tag); + S_->Require(silt_frac_key_, tag) + .SetMesh(subsurf_mesh)->AddComponent("cell", AmanziMesh::CELL, 1); + S_->RequireEvaluator(clay_frac_key_, tag); + S_->Require(clay_frac_key_, tag) + .SetMesh(subsurf_mesh)->AddComponent("cell", AmanziMesh::CELL, 1); + + // requirements: surface properties + S_->RequireEvaluator(color_index_key_, tag); + S_->Require(color_index_key_, tag) + .SetMesh(mesh_)->AddComponent("cell", AmanziMesh::CELL, 1); + S_->RequireEvaluator(pft_index_key_, tag); + S_->Require(pft_index_key_, tag) + .SetMesh(mesh_)->AddComponent("cell", AmanziMesh::CELL, 1); +} + + +// -- Initialize owned (dependent) variables. +void +SurfaceBalanceELMKernels::Initialize() +{ + PK_Physical_Default::Initialize(); + InitializeELMKernels_(tag_next_); + InitializePrimaryVariables_(tag_next_); +} + + +void +SurfaceBalanceELMKernels::InitializeELMKernels_(const Tag& tag) +{ + // Initialize the ELMKernels instance + Teuchos::ParameterList& ic_list = plist_->sublist("initial condition"); + double snow_depth = ic_list.get("initial snow depth [m]"); + double temp = ic_list.get("initial temperature [K]"); + double year = ic_list.get("initial time [yr]"); + //ATS::ELMKernels::set_zero_time(year); + //ATS::ELMKernels::set_initial_state(temp, snow_depth); + + // lat/lon + auto latlon = plist_->get>("latitude,longitude [degrees]"); + int ncols = mesh_->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED); + double latlon_arr[ncols][2]; + for (int i=0; i!=ncols; ++i) { + latlon_arr[i][0] = latlon[0]; + latlon_arr[i][1] = latlon[1]; + } + + // soil properties + S_->GetEvaluator(sand_frac_key_, tag).Update(*S_, name_); + auto& sand = *S_->Get(sand_frac_key_, tag).ViewComponent("cell", false); + S_->GetEvaluator(silt_frac_key_, tag).Update(*S_, name_); + auto& silt = *S_->Get(silt_frac_key_, tag).ViewComponent("cell", false); + S_->GetEvaluator(clay_frac_key_, tag).Update(*S_, name_); + auto& clay = *S_->Get(clay_frac_key_, tag).ViewComponent("cell", false); + + S_->GetEvaluator(color_index_key_, tag).Update(*S_, name_); + auto& color_index_tmp = *S_->Get(color_index_key_, tag).ViewComponent("cell", false); + S_->GetEvaluator(pft_index_key_, tag).Update(*S_, name_); + auto& pft_index_tmp = *S_->Get(pft_index_key_, tag).ViewComponent("cell", false); + + std::vector color_index(ncols); + for (int i=0; i!=ncols; ++i) color_index[i] = std::round(color_index_tmp[0][i]); + + double pft_fraction[ncols][NUM_LC_CLASSES]; + for (int i=0; i!=ncols; ++i) { + for (int j=0; j!=NUM_LC_CLASSES; ++j) { + pft_fraction[i][j] = j == std::round(pft_index_tmp[0][i]) ? 1. : 0.; + } + } + //ATS::ELMKernels::set_ground_properties(&latlon_arr[0][0], sand, clay, color_index, &pft_fraction[0][0]); + + // ELMKernels setup stage + //ATS::ELMKernels::setup_begin(); + auto subsurf_mesh = S_->GetMesh(domain_ss_); + Epetra_MultiVector dz(subsurf_mesh->cell_map(false), 1); + for (int col=0; col!=ncols; ++col) { + auto& faces = subsurf_mesh->faces_of_column(col); + auto& cells = subsurf_mesh->cells_of_column(col); + for (int i=0; i!=cells.size(); ++i) { + dz[0][cells[i]] = subsurf_mesh->face_centroid(faces[i])[2] - + subsurf_mesh->face_centroid(faces[i+1])[2]; + AMANZI_ASSERT(dz[0][cells[i]] > 0.); + } + } + //ATS::ELMKernels::set_dz(dz); + //ATS::ELMKernels::set_et_controls(1, 2, 0.1, 1.0, 0.1); + //ATS::ELMKernels::setup_end(); + //ATS::ELMKernels::set_dz(dz); +} + + +void +SurfaceBalanceELMKernels::InitializePrimaryVariables_(const Tag& tag) +{ + // set as intialized the sources + S_->GetW(surf_water_src_key_, tag, name_).PutScalar(0.); + S_->GetRecordW(surf_water_src_key_, tag, name_).set_initialized(); + S_->GetW(ss_water_src_key_, tag, name_).PutScalar(0.); + S_->GetRecordW(ss_water_src_key_, tag, name_).set_initialized(); + + // set as intialized the diagnostics + // S_->GetRecordW(evap_flux_key_, tag, name_)->set_initialized(); + S_->GetRecordW(qE_lh_key_, tag, name_).set_initialized(); + S_->GetRecordW(qE_sh_key_, tag, name_).set_initialized(); + S_->GetRecordW(qE_lw_out_key_, tag, name_).set_initialized(); + S_->GetRecordW(qE_cond_key_, tag, name_).set_initialized(); + S_->GetRecordW(snow_swe_key_, tag, name_).set_initialized(); + S_->GetRecordW(can_wc_key_, tag, name_).set_initialized(); + S_->GetRecordW(surf_temp_key_, tag, name_).set_initialized(); + S_->GetRecordW(soil_temp_key_, tag, name_).set_initialized(); + S_->GetRecordW(can_temp_key_, tag, name_).set_initialized(); +} + + +bool +SurfaceBalanceELMKernels::AdvanceStep(double t_old, double t_new, bool reinit) +{ + Teuchos::OSTab tab = vo_->getOSTab(); + + bool debug = false; + Teuchos::RCP dcvo = Teuchos::null; + int rank = mesh_->get_comm()->MyPID(); + double dt = t_new -t_old; + AMANZI_ASSERT(std::abs(dt - dt_) < 1.e-4); + + if (vo_->os_OK(Teuchos::VERB_LOW)) + *vo_->os() << "----------------------------------------------------------------" << std::endl + << "Advancing: t0 = " << S_->get_time(tag_current_) + << " t1 = " << S_->get_time(tag_next_) << " h = " << dt << std::endl + << "----------------------------------------------------------------" << std::endl; + + Tag tag = tag_current_; + + // Set the state + S_->GetEvaluator(pres_key_, tag).Update(*S_, name_); + const Epetra_MultiVector& pres = *S_->Get(pres_key_, tag).ViewComponent("cell", false); + S_->GetEvaluator(poro_key_, tag).Update(*S_, name_); + const Epetra_MultiVector& poro = *S_->Get(poro_key_, tag).ViewComponent("cell", false); + S_->GetEvaluator(sl_key_, tag).Update(*S_, name_); + const Epetra_MultiVector& sl = *S_->Get(sl_key_, tag).ViewComponent("cell", false); + + double patm = S_->Get("atmospheric_pressure", Tags::DEFAULT); + + //ATS::ELMKernels::set_wc(poro, sl); + //ATS::ELMKernels::set_tksat_from_porosity(poro); + //ATS::ELMKernels::set_pressure(pres, patm); + + // set the forcing + S_->GetEvaluator(met_sw_key_, tag).Update(*S_, name_); + const Epetra_MultiVector& met_sw = *S_->Get(met_sw_key_, tag).ViewComponent("cell", false); + S_->GetEvaluator(met_lw_key_, tag).Update(*S_, name_); + const Epetra_MultiVector& met_lw = *S_->Get(met_lw_key_, tag).ViewComponent("cell", false); + S_->GetEvaluator(met_air_temp_key_, tag).Update(*S_, name_); + const Epetra_MultiVector& met_air_temp = *S_->Get(met_air_temp_key_, tag).ViewComponent("cell", false); + S_->GetEvaluator(met_rel_hum_key_, tag).Update(*S_, name_); + const Epetra_MultiVector& met_rel_hum = *S_->Get(met_rel_hum_key_, tag).ViewComponent("cell", false); + S_->GetEvaluator(met_wind_speed_key_, tag).Update(*S_, name_); + const Epetra_MultiVector& met_wind_speed = *S_->Get(met_wind_speed_key_, tag).ViewComponent("cell", false); + S_->GetEvaluator(met_prain_key_, tag).Update(*S_, name_); + const Epetra_MultiVector& met_prain = *S_->Get(met_prain_key_, tag).ViewComponent("cell", false); + S_->GetEvaluator(met_psnow_key_, tag).Update(*S_, name_); + const Epetra_MultiVector& met_psnow = *S_->Get(met_psnow_key_, tag).ViewComponent("cell", false); + + //ATS::ELMKernels::set_met_data(met_sw, met_lw, met_prain, met_psnow, met_air_temp, met_rel_hum, met_wind_speed, patm); + + // set the start time, endtime + //ATS::ELMKernels::advance_time(S_->get_cycle(tag), t_old, dt); // units in seconds + + // get diagnostics + Epetra_MultiVector& qE_lh = *S_->GetW(qE_lh_key_, tag, name_) + .ViewComponent("cell", false); + Epetra_MultiVector& qE_sh = *S_->GetW(qE_sh_key_, tag, name_) + .ViewComponent("cell", false); + Epetra_MultiVector& qE_lw_out = *S_->GetW(qE_lw_out_key_, tag, name_) + .ViewComponent("cell", false); + Epetra_MultiVector& qE_cond = *S_->GetW(qE_cond_key_, tag, name_) + .ViewComponent("cell", false); + //ATS::ELMKernels::get_ground_energy_fluxes(qE_lh, qE_sh, qE_lw_out, qE_cond); + + Epetra_MultiVector& snow_depth = *S_->GetW(key_, tag, name_) + .ViewComponent("cell", false); + Epetra_MultiVector& snow_swe = *S_->GetW(snow_swe_key_, tag, name_) + .ViewComponent("cell", false); + Epetra_MultiVector& can_wc = *S_->GetW(can_wc_key_, tag, name_) + .ViewComponent("cell", false); + Epetra_MultiVector& surf_temp = *S_->GetW(surf_temp_key_, tag, name_) + .ViewComponent("cell", false); + Epetra_MultiVector& soil_temp = *S_->GetW(soil_temp_key_, tag, name_) + .ViewComponent("cell", false); + Epetra_MultiVector& can_temp = *S_->GetW(can_temp_key_, tag, name_) + .ViewComponent("cell", false); + //ATS::ELMKernels::get_diagnostics(snow_swe, snow_depth, can_wc, surf_temp, can_temp, soil_temp); + changedEvaluatorPrimary(key_, tag, *S_); + + // get output + Epetra_MultiVector& surf_water_src = *S_->GetW(surf_water_src_key_, tag, name_) + .ViewComponent("cell", false); + Epetra_MultiVector& ss_water_src = *S_->GetW(ss_water_src_key_, tag, name_) + .ViewComponent("cell", false); + //ATS::ELMKernels::get_total_mass_fluxes(surf_water_src, ss_water_src); + changedEvaluatorPrimary(surf_water_src_key_, tag, *S_); + changedEvaluatorPrimary(ss_water_src_key_, tag, *S_); + + if (vo_->os_OK(Teuchos::VERB_HIGH)) { + std::vector vnames; + std::vector< Teuchos::Ptr > vecs; + + vnames.push_back("inc shortwave radiation [W/m^2]"); + vecs.push_back(S_->GetPtr(met_sw_key_, tag).ptr()); + vnames.push_back("inc longwave radiation [W/m^2]"); + vecs.push_back(S_->GetPtr(met_lw_key_, tag).ptr()); + vnames.push_back("inc latent heat [W/m^2]"); + vecs.push_back(S_->GetPtr(qE_lh_key_, tag).ptr()); + vnames.push_back("inc sensible heat [W/m^2]"); + vecs.push_back(S_->GetPtr(qE_sh_key_, tag).ptr()); + vnames.push_back("out longwave radiation [W/m^2]"); + vecs.push_back(S_->GetPtr(qE_lw_out_key_, tag).ptr()); + vnames.push_back("out conducted soil [W/m^2]"); + vecs.push_back(S_->GetPtr(qE_cond_key_, tag).ptr()); + + db_->WriteVectors(vnames, vecs, true); + db_->WriteDivider(); + + vnames.clear(); + vecs.clear(); + + vnames.push_back("surface water source [m/s]"); + vecs.push_back(S_->GetPtr(surf_water_src_key_, tag).ptr()); + db_->WriteVectors(vnames, vecs, true); + db_->WriteDivider(); + + vnames.clear(); + vecs.clear(); + + vnames.push_back("snow depth [m]"); + vecs.push_back(S_->GetPtr(key_, tag).ptr()); + vnames.push_back("snow swe [m]"); + vecs.push_back(S_->GetPtr(snow_swe_key_, tag).ptr()); + vnames.push_back("canopy storage [m]"); + vecs.push_back(S_->GetPtr(can_wc_key_, tag).ptr()); + vnames.push_back("skin temperature [K]"); + vecs.push_back(S_->GetPtr(surf_temp_key_, tag).ptr()); + vnames.push_back("leaf temperature [K]"); + vecs.push_back(S_->GetPtr(can_temp_key_, tag).ptr()); + db_->WriteVectors(vnames, vecs, true); + db_->WriteDivider(); + } + return false; +} + + +} // namespace +} // namespace diff --git a/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.hh b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.hh new file mode 100644 index 000000000..d1b9f29db --- /dev/null +++ b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.hh @@ -0,0 +1,107 @@ +/* -*- mode: c++; indent-tabs-mode: nil -*- */ + +/* ------------------------------------------------------------------------- + ATS + + License: see $ATS_DIR/COPYRIGHT + + PK utilizing surface mass and energy physics kernels from the + ELMKernels library + + ------------------------------------------------------------------------- */ + +#ifndef PK_SURFACE_BALANCE_ELMKERNELS_HH_ +#define PK_SURFACE_BALANCE_ELMKERNELS_HH_ + +#include "PK_Factory.hh" +#include "pk_physical_bdf_default.hh" + +namespace Amanzi { +namespace SurfaceBalance { + +class SurfaceBalanceELMKernels : public PK_Physical_Default { + +public: + + SurfaceBalanceELMKernels(Teuchos::ParameterList& pk_tree, + const Teuchos::RCP& global_list, + const Teuchos::RCP& S, + const Teuchos::RCP& solution); + + // main methods + // -- Setup data. + virtual void Setup() override; + + // -- Initialize owned (dependent) variables. + virtual void Initialize() override; + + + // // -- Commit any secondary (dependent) variables. + // virtual void CommitStep(double t_old, double t_new, ); + + // -- Calculate any diagnostics prior to doing vis + virtual void CalculateDiagnostics(const Tag& tag) override {} + + virtual void set_dt(double dt) override { + AMANZI_ASSERT(std::abs(dt - dt_) < 1.e-4); + } + virtual double get_dt() override { return dt_; } + + // Advance PK from time t_old to time t_new. True value of the last + // parameter indicates drastic change of boundary and/or source terms + // that may need PK's attention. + virtual bool AdvanceStep(double t_old, double t_new, bool reinit) override; + + protected: + void SetupDependencies_(const Tag& tag); + void InitializeELMKernels_(const Tag& tag); + void InitializePrimaryVariables_(const Tag& tag); + + protected: + Key domain_ss_; + Key domain_snow_; + Key domain_can_; + + double dt_; + + Key surf_water_src_key_; + Key ss_water_src_key_; + + Key met_sw_key_; + Key met_lw_key_; + Key met_air_temp_key_; + Key met_rel_hum_key_; + Key met_wind_speed_key_; + Key met_prain_key_; + Key met_psnow_key_; + + Key qE_lh_key_; + Key qE_sh_key_; + Key qE_lw_out_key_; + Key qE_cond_key_; + + Key snow_swe_key_; + Key can_wc_key_; + Key surf_temp_key_; + Key soil_temp_key_; + Key can_temp_key_; + + Key pres_key_; + Key poro_key_; + Key sl_key_; + + Key sand_frac_key_; + Key silt_frac_key_; + Key clay_frac_key_; + Key color_index_key_; + Key pft_index_key_; + + private: + // factory registration + static RegisteredPKFactory reg_; +}; + +} // namespace AmanziFlow +} // namespace Amanzi + +#endif diff --git a/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels_reg.hh b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels_reg.hh new file mode 100644 index 000000000..0238910ab --- /dev/null +++ b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels_reg.hh @@ -0,0 +1,19 @@ +/* -*- mode: c++; indent-tabs-mode: nil -*- */ +/* ------------------------------------------------------------------------- + * ATS + * + * License: see $ATS_DIR/COPYRIGHT + * + * ------------------------------------------------------------------------- */ + + +#include "surface_balance_ELMKernels.hh" + +namespace Amanzi { +namespace SurfaceBalance { + +RegisteredPKFactory +SurfaceBalanceELMKernels::reg_("surface balance ELMKernels"); + +} // namespace +} // namespace From 11907e930f3b09d8b7c043c404f410df2331d924 Mon Sep 17 00:00:00 2001 From: Joe Beisman Date: Thu, 10 Nov 2022 02:37:05 -0500 Subject: [PATCH 2/7] started elm kernels interface --- .../ELMKernels/ats_elm_kernels_interface.hh | 291 ++++++++++++++++++ .../ELMKernels/canopy_hydrology_evaluator.cc | 78 +++++ .../ELMKernels/canopy_hydrology_evaluator.hh | 53 ++++ 3 files changed, 422 insertions(+) create mode 100644 src/pks/surface_balance/ELMKernels/ats_elm_kernels_interface.hh create mode 100644 src/pks/surface_balance/ELMKernels/canopy_hydrology_evaluator.cc create mode 100644 src/pks/surface_balance/ELMKernels/canopy_hydrology_evaluator.hh diff --git a/src/pks/surface_balance/ELMKernels/ats_elm_kernels_interface.hh b/src/pks/surface_balance/ELMKernels/ats_elm_kernels_interface.hh new file mode 100644 index 000000000..9b0efe1ee --- /dev/null +++ b/src/pks/surface_balance/ELMKernels/ats_elm_kernels_interface.hh @@ -0,0 +1,291 @@ + +#ifndef ATS_CLM_INTERFACE_HH_ +#define ATS_CLM_INTERFACE_HH_ + +#include +#include + +#include "Epetra_MultiVector.h" + +#define NUM_LC_CLASSES 18 + +namespace ATS { +namespace ELMKernels { + + class ELMKernelDriver { + + std::shared_ptr ELMState_; + + // time-invariant, no spatial aspect + std::shared_ptr> snicar_data_; + std::shared_ptr> snw_rds_table_; + std::shared_ptr> pft_data_; + + // time variable, but no spatial aspect, unless domain very large + // holds 12 monthly values for aerosol dep + // interpolates input file using lat and lon + std::shared_ptr> aerosol_data_; + + // time and space variable - must be time-interpolated at every step + std::shared_ptr atm_forcing_; + std::shared_ptr> phen_data_; + + // derived from aerosol_data and snow physics + std::shared_ptr> aerosol_masses_; + std::shared_ptr> aerosol_concentrations_; + + // the path/to/file for phenology and atm forcing inputs + std::string fname_surfdata_; + std::string fname_forc_; + + // TODO better + // number of columns, cell per column + int ncols_{0}; + int cell_per_col{0}; + + // TODO should simplify this + // don't need 2D domain decomp with ATS + // will also likely have ATS do the file reads + DomainDecomposition<2> dd_; + + // for phenology input data - should hide this in phen_data + std::unordered_map host_phen_views_; + + ELM::Utils::Date start_; + + +// mapping between ATS and ELM +// if ATS list of local nodes for each rank isn't contiguously ordered +// use a std::map +// where map[0...ncols] = ATS node id +// so ELM_var(i) = ATS_var(map[i]) +// this should be fast to iterate? + +// SerialDenseVec? +// soil temp solve? +// mesh definition - pass in mesh_, or just vector of interfaces (zisoi) to get dz, zsoi, zisoi +// slope aspect + +// need: +// Land +// lat +// lon +// lat_r +// lon_r +// dewmx +// irrig_rate +// n_irrig_steps_left +// oldfflag +// veg_active +// do_capsnow +// topo_slope +// topo_std +// t_h2osfc +// tsoi +// altmax_indx +// altmax_lastyear_indx +// t10 +// t_veg +// dz | +// zsoi | - assigned together +// zisoi | +// h2osoi_ice | +// h2osoi_liq | - together? +// h2osoi_vol | + +// need to make phenology_data and helpers inot it's own call +// make AtmForc read start date from the files, instead of providing +// should make a decision on clock sync +// put coszen into functions, also max_dayl and dayl + + + + + public: + + int init_lsm(int ncols, const MPI_Comm& comm, std::string_view fsurfdata, + std::string_view fsnicar, std::string_view fforc, + std::string_view fparam, std::string_view faerosol, + std::string_view fsnowage) { + // assign filenames + fname_surfdata_{fsurfdata}; + fname_forc_{fforc}; + + ncols_ = ncols; + cell_per_col = nlevgrnd; // hardwired in various places within ELMKernels + // may break if not set at 15 + + // assert a bunch to make sure everything lines up + + // need to figure out mapping + // implement better way of defining local/global problem + // knowledge of global problem is currently required for + // pnetcdf usage - if ATS reads, we don't need it + dd_ = ELM::Utils::create_domain_decomposition_2D( + { ncols, ncols }, + { 1, 1 }, + { 0, 0 }); + dd_.comm = comm; + + // bloated state struct - will move some of this data + ELMState_ = std::make_shared(ncols); + + // time invariant constants + // without spatial awareness + snicar_data_ = std::make_shared>(); + snw_rds_table_ = std::make_shared>(); + pft_data_ = std::make_shared>(); + + // one year of aerosol deposition from nearest gridcell + // currently invariant after this read + aerosol_data_ = std::make_shared>(); + + ELM::initialize_kokkos_elm(*ELMState_, *snicar_data_, *snw_rds_table_, *pft_data_, + *aerosol_data_, dd_, fname_surfdata_, fname_param_, + fname_snicar_, fname_snowage_, fname_aerosol_); + return 0; + } + + + // initialization of state fields + // must deep copy all of these into kokkos + + // can also be provided on cell-by-cell basis, + // which will be necessary for restarts + int init_snow(double snow_depth, double snow_temp) { + return 0; + } + int init_snow(const Epetra_MultiVector& snow_depth) { + return 0; + } // snow temp needs to be taken care of here as well + + // should probably be run after snl calculated + int init_temperature(double const_temp) { + // deep copy to State view + assign(ELMState_->t_h2osfc, const_temp); + // copy constant value into subsurface only + // snow layers should be initialized as 0.0 + auto h_tsoi = Kokkos::create_mirror_view(ELMState_->t_soisno); + for (int i = 0; i < ncols_; ++i) { + for (int j = nlevgrnd; j < nlevsno + nlevgrnd; ++j) { + h_tsoi(i, j) = const_temp; + } + } + assign(ELMState_->t_soisno, h_tsoi); + return 0; + } + + // should make a decision on clock sync + int init_start_time(int year, int month, int day) { + start_ = ELM::Utils::Date(year, month, day); + return 0; + } + + // just pass in mesh_? or zisoi? + int init_mesh(const Epetra_MultiVector& dz) { + return 0; + } // should be constant with 1 column worth of dz, or zisoi + + + + + + + int set_energy_state(const Epetra_MultiVector& soil_temp) { + + // place soil_temp from ATS into t_soisno once mapping decided upon + return 0; + } + + int set_water_state(const Epetra_MultiVector& saturation, const Epetra_MultiVector& porosity) { + // place saturation from ATS into watsat once mapping decided upon + // h2osoi_liq, h2osoi_ice?, watsat + return 0; + } + + + int transfer_function_VG_to_CH() { + // we're ignoring pressure, so this probably isn't needed + return 0; + } + + + int setup_lsm() { + + int atm_nsteps = 101; // how many timesteps to read? + const auto fstart = ELM::Utils::Date(1985, 1, 1); // make atm_forcing read this - may already be done + + // need to interpolate forcing at every step and + // reread from file periodically - happens automatically when forcing step (atm_nsteps-1)/atm_nsteps has been consumed + atm_forcing_ = std::make_shared(fname_forc_, fstart, atm_nsteps, ncols); + + // phenology data manager + // make host mirrors - need to be persistent + phen_data_ = std::make_shared>(dd, ncells, 17); + host_phen_views_ = get_phen_host_views(*phen_data); + + // containers for aerosol deposition and concentration within snowpack layers + aerosol_masses_ = std::make_shared>(ncells); + aerosol_concentrations_ = std::make_shared>(ncells); + } + + + + int get_water_flux_diagnostic() { + // get some/all of the optional qflx outputs + // need to calc qflx_infiltration + // with qflx_top_soil, qflx_snomelt, etc + // from SoilHydrologyMod.F90 + // qflx_snow_grnd + // qflx_rain_grnd + // qflx_snow_melt + // qflx_evap_tot + // qflx_evap_soil + // qflx_evap_veg + // qflx_tran_veg + // qflx_irr + // irr_flag + return 0; + } + + int get_water_flux_total() { + // get the accumulated qflx outputs + // convert subsurf flux mm/s to 1/s -- mm to m 1e3 + // keep sfc flux in L/T + return 0; + } + + int get_energy_flux_diagnostic() { + // get some/all of the optional eflx outputs + //eflx_sh_grnd + //eflx_sh_snow + //eflx_sh_soil + //eflx_sh_h2osfc + return 0; + } + + int get_energy_flux_total() { + // get the accumulated eflx outputs + // eflx_lh_tot + // eflx_sh_tot + // eflx_lwrad_out - or eflx_lwrad_net?? + // eflx_soil + return 0; + } + + + int advance_time(double time, double dt) { + // call 2 functions - prepare init_timestep() and run_lsm_physics + return 0; + } + + + }; + + +} // namespace ATS +} // namespace ELMKernels + + + +#endif diff --git a/src/pks/surface_balance/ELMKernels/canopy_hydrology_evaluator.cc b/src/pks/surface_balance/ELMKernels/canopy_hydrology_evaluator.cc new file mode 100644 index 000000000..2993c5019 --- /dev/null +++ b/src/pks/surface_balance/ELMKernels/canopy_hydrology_evaluator.cc @@ -0,0 +1,78 @@ + +#include "Key.hh" + +namespace Amanzi { +namespace ELMKernels { + +CanopyHydrologyEvaluator::CanopyHydrologyEvaluator(Teuchos::ParameterList& plist) : + EvaluatorSecondaryMonotypeCV(plist) +{ + // determine the domain + Key akey = my_keys_.front().first; + domain_ = Keys::getDomain(akey); + akey = Keys::getVarName(akey); + domain_snow_ = Keys::readDomainHint(plist_, domain_, "surface", "snow"); + auto tag = my_keys_.front().second; + + // my keys + // -- sources + // gather all dependencies, assign any default parameters + // would prefer to do this as simply as possible, + // because the function signatures will change soon +} + +// Required methods from EvaluatorSecondaryMonotypeCV +void +CanopyHydrologyEvaluator::Evaluate_(const State& S, + const std::vector& results) +{ + auto mesh = S.GetMesh(domain_); + auto tag = my_keys_.front().second; + + // collect dependencies + + // collect output vecs + + // call canopy_hydrology kernels + ELM::kokkos_canopy_hydrology(*S, atm_forcing->forc_PREC, dtime, time_plus_half_dt) + +} + +void +CanopyHydrologyEvaluator::EvaluatePartialDerivative_(const State& S, + const Key& wrt_key, const Tag& wrt_tag, + const std::vector& results) {//nope} + + +// custom EC used to set subfield names +void +CanopyHydrologyEvaluator::EnsureCompatibility_Structure_(State& S) +{ + S.GetRecordSetW(my_keys_.front().first).set_subfieldnames( + {"bare", "water", "snow"}); +} + + +void CanopyHydrologyEvaluator::EnsureCompatibility_ToDeps_(State& S) +{ + // new state! + if (land_cover_.size() == 0) + land_cover_ = getLandCover(S.ICList().sublist("land cover types"), + {"albedo_ground", "emissivity_ground"}); + + for (auto dep : dependencies_) { + auto& fac = S.Require(dep.first, dep.second); + if (Keys::getDomain(dep.first) == domain_snow_) { + fac.SetMesh(S.GetMesh(domain_snow_)) + ->SetGhosted() + ->AddComponent("cell", AmanziMesh::CELL, 1); + } else { + fac.SetMesh(S.GetMesh(domain_)) + ->SetGhosted() + ->AddComponent("cell", AmanziMesh::CELL, 1); + } + } +} + +} // namespace Relations +} // namespace ELMKernels diff --git a/src/pks/surface_balance/ELMKernels/canopy_hydrology_evaluator.hh b/src/pks/surface_balance/ELMKernels/canopy_hydrology_evaluator.hh new file mode 100644 index 000000000..8e09e2e1c --- /dev/null +++ b/src/pks/surface_balance/ELMKernels/canopy_hydrology_evaluator.hh @@ -0,0 +1,53 @@ +#pragma once + +#include "Factory.hh" +#include "EvaluatorSecondaryMonotype.hh" +#include "LandCover.hh" + +namespace Amanzi { +namespace ELMKernels { + +class CanopyHydrologyEvaluator : public EvaluatorSecondaryMonotypeCV { + public: + explicit CanopyHydrologyEvaluator(Teuchos::ParameterList& plist); + CanopyHydrologyEvaluator(const CanopyHydrologyEvaluator& other) = default; + + virtual Teuchos::RCP Clone() const override { + return Teuchos::rcp(new CanopyHydrologyEvaluator(*this)); + } + + protected: + // custom EC used to set subfield names + virtual void EnsureCompatibility_Structure_(State& S) override; + + // custom EC used because deps have 1 component not 3 + virtual void EnsureCompatibility_ToDeps_(State& S) override; + + // Required methods from EvaluatorSecondaryMonotypeCV + virtual void Evaluate_(const State& S, + const std::vector& results) override; + + virtual void EvaluatePartialDerivative_(const State& S, + const Key& wrt_key, const Tag& wrt_tag, + const std::vector& results) override; + + protected: + Key domain_; + Key domain_snow_; + + Key albedo_key_, emissivity_key_; + Key snow_dens_key_; + Key unfrozen_fraction_key_; + + double a_water_, a_ice_; + double e_water_, e_ice_, e_snow_; + + // this is horrid, because this cannot yet live in state + // bring on new state! + + private: + static Utils::RegisteredFactory reg_; +}; + +} // namespace Relations +} // namespace ELMKernels From a326ed36c454435b2174e6cd0387a6711606619e Mon Sep 17 00:00:00 2001 From: Joe Beisman Date: Thu, 15 Dec 2022 09:34:30 -0500 Subject: [PATCH 3/7] hardwired to build with ELMKernels library --- src/executables/CMakeLists.txt | 5 +++++ src/pks/surface_balance/CMakeLists.txt | 15 ++++++++++++--- .../ELMKernels/surface_balance_ELMKernels.cc | 4 ++++ 3 files changed, 21 insertions(+), 3 deletions(-) diff --git a/src/executables/CMakeLists.txt b/src/executables/CMakeLists.txt index adf7e704f..68a62f573 100644 --- a/src/executables/CMakeLists.txt +++ b/src/executables/CMakeLists.txt @@ -122,6 +122,10 @@ set(ats_link_libs ) +set(ELM_LIBRARIES $ENV{ELM_DIR}/install/lib/libelm_kokkos.dylib) +set(ELM_INCLUDES_DIR $ENV{ELM_DIR}/install/include) +include_directories(${ELM_INCLUDES_DIR}) + # note, we can be inclusive here, because if they aren't enabled, # these won't be defined and will result in empty strings. set(tpl_link_libs @@ -136,6 +140,7 @@ set(tpl_link_libs ${HYPRE_LIBRARIES} ${HDF5_LIBRARIES} ${CLM_LIBRARIES} + ${ELM_LIBRARIES} ) add_amanzi_library(ats_executable diff --git a/src/pks/surface_balance/CMakeLists.txt b/src/pks/surface_balance/CMakeLists.txt index c16cd59dd..f643428f9 100644 --- a/src/pks/surface_balance/CMakeLists.txt +++ b/src/pks/surface_balance/CMakeLists.txt @@ -1,6 +1,13 @@ # -*- mode: cmake -*- # ATS Surface balance PKs describe Evaporation, energy fluxes from # long/showtwave radiation, precip, etc etc etc + +# temporary ELMKernels +set(ELM_LIBRARIES $ENV{ELM_DIR}/install/lib/libelm_kokkos.dylib) +set(ELM_INCLUDES_DIR $ENV{ELM_DIR}/install/include) +include_directories(${ELM_INCLUDES_DIR}) + + include_directories(${ATS_SOURCE_DIR}/src/pks) include_directories(${ATS_SOURCE_DIR}/src/constitutive_relations/surface_subsurface_fluxes) include_directories(${CMAKE_CURRENT_SOURCE_DIR}/constitutive_relations/land_cover) @@ -51,7 +58,8 @@ if (ENABLE_CLM) ) endif() -if (ENABLE_ELMKERNELS) +set(ENABLE_ELM true) +if (ENABLE_ELM) list(APPEND ats_surface_balance_src_files ELMKernels/surface_balance_ELMKernels.cc ) @@ -101,7 +109,7 @@ if (ENABLE_CLM) ) endif() -if (ENABLE_ELMKERNELS) +if (ENABLE_ELM) list(APPEND ats_surface_balance_inc_files ELMKernels/surface_balance_ELMKernels.hh ) @@ -111,6 +119,7 @@ set(ats_surface_balance_link_libs ${Teuchos_LIBRARIES} ${Epetra_LIBRARIES} ${CLM_LIBRARIES} + ${ELM_LIBRARIES} error_handling atk mesh @@ -157,7 +166,7 @@ if (ENABLE_CLM) ) endif() -if (ENABLE_ELMKERNELS) +if (ENABLE_ELM) register_evaluator_with_factory( HEADERFILE ELMKernels/surface_balance_ELMKernels_reg.hh LISTNAME ATS_SURFACE_BALANCE_REG diff --git a/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc index 7e9130fbc..08a8a2e7a 100644 --- a/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc +++ b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc @@ -10,6 +10,10 @@ #include "pk_helpers.hh" #include "surface_balance_ELMKernels.hh" +#include "canopy_fluxes_kokkos.hh" + +#define NUM_LC_CLASSES 18 + namespace Amanzi { namespace SurfaceBalance { From 0404040f02178ae410bae7f858a3068e4769f164 Mon Sep 17 00:00:00 2001 From: Joe Beisman Date: Thu, 15 Dec 2022 14:52:50 -0500 Subject: [PATCH 4/7] call ELM functionality from ATS --- .../ELMKernels/surface_balance_ELMKernels.cc | 10 ++++++---- .../ELMKernels/surface_balance_ELMKernels.hh | 3 +++ 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc index 08a8a2e7a..01c7ee437 100644 --- a/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc +++ b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc @@ -10,7 +10,7 @@ #include "pk_helpers.hh" #include "surface_balance_ELMKernels.hh" -#include "canopy_fluxes_kokkos.hh" +//#include "elm_kokkos_interface.hh" #define NUM_LC_CLASSES 18 @@ -102,9 +102,7 @@ SurfaceBalanceELMKernels::Setup() SetupDependencies_(tag_next_); // Set up the ELMKernels object - //ATS::ELMKernels::init(subsurf_mesh->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED), - // mesh_->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED), - // 2, mesh_->get_comm()->MyPID(), 3); + elm_ = std::make_unique(mesh_->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED)); } void @@ -231,6 +229,7 @@ SurfaceBalanceELMKernels::Initialize() void SurfaceBalanceELMKernels::InitializeELMKernels_(const Tag& tag) { + elm_->setup(); // Initialize the ELMKernels instance Teuchos::ParameterList& ic_list = plist_->sublist("initial condition"); double snow_depth = ic_list.get("initial snow depth [m]"); @@ -334,6 +333,9 @@ SurfaceBalanceELMKernels::AdvanceStep(double t_old, double t_new, bool reinit) Tag tag = tag_current_; + auto dummy_date = ELM::Utils::Date(2014, 1, 1); + elm_->advance(dummy_date, dt); + // Set the state S_->GetEvaluator(pres_key_, tag).Update(*S_, name_); const Epetra_MultiVector& pres = *S_->Get(pres_key_, tag).ViewComponent("cell", false); diff --git a/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.hh b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.hh index d1b9f29db..142b9cdbd 100644 --- a/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.hh +++ b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.hh @@ -15,6 +15,7 @@ #include "PK_Factory.hh" #include "pk_physical_bdf_default.hh" +#include "elm_kokkos_interface.hh" namespace Amanzi { namespace SurfaceBalance { @@ -96,6 +97,8 @@ public: Key color_index_key_; Key pft_index_key_; + std::unique_ptr elm_; + private: // factory registration static RegisteredPKFactory reg_; From 5867133e96bb7e53216433ce2574fbaaf6a12e56 Mon Sep 17 00:00:00 2001 From: Joe Beisman Date: Thu, 5 Jan 2023 15:03:11 -0500 Subject: [PATCH 5/7] progress toward ELM interface --- .../ELMKernels/ats_elm_interface.cc | 24 +++++++++++++++ .../ELMKernels/ats_elm_interface.hh | 30 +++++++++++++++++++ ...nels_interface.hh => interface_archive.hh} | 0 .../ELMKernels/surface_balance_ELMKernels.cc | 6 ++-- .../ELMKernels/surface_balance_ELMKernels.hh | 10 +++---- 5 files changed, 62 insertions(+), 8 deletions(-) create mode 100644 src/pks/surface_balance/ELMKernels/ats_elm_interface.cc create mode 100644 src/pks/surface_balance/ELMKernels/ats_elm_interface.hh rename src/pks/surface_balance/ELMKernels/{ats_elm_kernels_interface.hh => interface_archive.hh} (100%) diff --git a/src/pks/surface_balance/ELMKernels/ats_elm_interface.cc b/src/pks/surface_balance/ELMKernels/ats_elm_interface.cc new file mode 100644 index 000000000..37250b229 --- /dev/null +++ b/src/pks/surface_balance/ELMKernels/ats_elm_interface.cc @@ -0,0 +1,24 @@ +/* + 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. + +*/ + +#include + +#include "pk_helpers.hh" +#include "elm_kokkos_interface.hh" + +#define NUM_LC_CLASSES 18 + +namespace Amanzi { +namespace LandPhysics { + + Kernels::Kernels(Teuchos::ParameterList& pk_tree, + const Teuchos::RCP& global_list, + const Teuchos::RCP& S, + const Teuchos::RCP& solution): {} + +} // namespace LandPhysics +} // namespace ATS diff --git a/src/pks/surface_balance/ELMKernels/ats_elm_interface.hh b/src/pks/surface_balance/ELMKernels/ats_elm_interface.hh new file mode 100644 index 000000000..7a33104fa --- /dev/null +++ b/src/pks/surface_balance/ELMKernels/ats_elm_interface.hh @@ -0,0 +1,30 @@ + +#ifndef ATS_ELM_KERNELS_INTERFACE_HH_ +#define ATS_ELM_KERNELS_INTERFACE_HH_ + +#include +#include + +#include "Epetra_MultiVector.h" + +#define NUM_LC_CLASSES 18 + +namespace ATS { +namespace LandPhysics { + + class Kernels { + public: + Kernels(); + Kernels(Teuchos::ParameterList& pk_tree, + const Teuchos::RCP& global_list, + const Teuchos::RCP& S, + const Teuchos::RCP& solution); + +protected: + + Teuchos::RCP elm_kernels_; + + }; // Kernels + +} // namespace LandPhysics +} // namespace ATS diff --git a/src/pks/surface_balance/ELMKernels/ats_elm_kernels_interface.hh b/src/pks/surface_balance/ELMKernels/interface_archive.hh similarity index 100% rename from src/pks/surface_balance/ELMKernels/ats_elm_kernels_interface.hh rename to src/pks/surface_balance/ELMKernels/interface_archive.hh diff --git a/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc index 01c7ee437..0874c49ea 100644 --- a/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc +++ b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc @@ -15,7 +15,7 @@ #define NUM_LC_CLASSES 18 namespace Amanzi { -namespace SurfaceBalance { +namespace LandPhysics { SurfaceBalanceELMKernels::SurfaceBalanceELMKernels(Teuchos::ParameterList& pk_tree, const Teuchos::RCP& global_list, @@ -101,8 +101,8 @@ SurfaceBalanceELMKernels::Setup() // set requirements on dependencies SetupDependencies_(tag_next_); - // Set up the ELMKernels object - elm_ = std::make_unique(mesh_->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED)); + // Set up the ELMKernels interface + elm_interface_ = Teuchos::rcp(new Amanzi::LandPhysics::Kernels(engine_name, engine_inputfile)); } void diff --git a/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.hh b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.hh index 142b9cdbd..57eeac2e7 100644 --- a/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.hh +++ b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.hh @@ -15,10 +15,10 @@ #include "PK_Factory.hh" #include "pk_physical_bdf_default.hh" -#include "elm_kokkos_interface.hh" +#include "ats_elm_interface.hh" namespace Amanzi { -namespace SurfaceBalance { +namespace LandPhysics { class SurfaceBalanceELMKernels : public PK_Physical_Default { @@ -97,14 +97,14 @@ public: Key color_index_key_; Key pft_index_key_; - std::unique_ptr elm_; + Teuchos::RCP elm_interface{nullptr}; private: // factory registration static RegisteredPKFactory reg_; }; -} // namespace AmanziFlow -} // namespace Amanzi +} // namespace LandPhysics +} // namespace ATS #endif From 0fa679bacf7683dc34590da29ade60429a4e545e Mon Sep 17 00:00:00 2001 From: Joe Beisman Date: Mon, 23 Oct 2023 20:21:44 -0400 Subject: [PATCH 6/7] Deleted some unused code --- .../ELMKernels/ats_elm_interface.cc | 24 -- .../ELMKernels/ats_elm_interface.hh | 30 -- .../ELMKernels/canopy_hydrology_evaluator.cc | 78 ----- .../ELMKernels/canopy_hydrology_evaluator.hh | 53 ---- .../ELMKernels/interface_archive.hh | 291 ------------------ 5 files changed, 476 deletions(-) delete mode 100644 src/pks/surface_balance/ELMKernels/ats_elm_interface.cc delete mode 100644 src/pks/surface_balance/ELMKernels/ats_elm_interface.hh delete mode 100644 src/pks/surface_balance/ELMKernels/canopy_hydrology_evaluator.cc delete mode 100644 src/pks/surface_balance/ELMKernels/canopy_hydrology_evaluator.hh delete mode 100644 src/pks/surface_balance/ELMKernels/interface_archive.hh diff --git a/src/pks/surface_balance/ELMKernels/ats_elm_interface.cc b/src/pks/surface_balance/ELMKernels/ats_elm_interface.cc deleted file mode 100644 index 37250b229..000000000 --- a/src/pks/surface_balance/ELMKernels/ats_elm_interface.cc +++ /dev/null @@ -1,24 +0,0 @@ -/* - 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. - -*/ - -#include - -#include "pk_helpers.hh" -#include "elm_kokkos_interface.hh" - -#define NUM_LC_CLASSES 18 - -namespace Amanzi { -namespace LandPhysics { - - Kernels::Kernels(Teuchos::ParameterList& pk_tree, - const Teuchos::RCP& global_list, - const Teuchos::RCP& S, - const Teuchos::RCP& solution): {} - -} // namespace LandPhysics -} // namespace ATS diff --git a/src/pks/surface_balance/ELMKernels/ats_elm_interface.hh b/src/pks/surface_balance/ELMKernels/ats_elm_interface.hh deleted file mode 100644 index 7a33104fa..000000000 --- a/src/pks/surface_balance/ELMKernels/ats_elm_interface.hh +++ /dev/null @@ -1,30 +0,0 @@ - -#ifndef ATS_ELM_KERNELS_INTERFACE_HH_ -#define ATS_ELM_KERNELS_INTERFACE_HH_ - -#include -#include - -#include "Epetra_MultiVector.h" - -#define NUM_LC_CLASSES 18 - -namespace ATS { -namespace LandPhysics { - - class Kernels { - public: - Kernels(); - Kernels(Teuchos::ParameterList& pk_tree, - const Teuchos::RCP& global_list, - const Teuchos::RCP& S, - const Teuchos::RCP& solution); - -protected: - - Teuchos::RCP elm_kernels_; - - }; // Kernels - -} // namespace LandPhysics -} // namespace ATS diff --git a/src/pks/surface_balance/ELMKernels/canopy_hydrology_evaluator.cc b/src/pks/surface_balance/ELMKernels/canopy_hydrology_evaluator.cc deleted file mode 100644 index 2993c5019..000000000 --- a/src/pks/surface_balance/ELMKernels/canopy_hydrology_evaluator.cc +++ /dev/null @@ -1,78 +0,0 @@ - -#include "Key.hh" - -namespace Amanzi { -namespace ELMKernels { - -CanopyHydrologyEvaluator::CanopyHydrologyEvaluator(Teuchos::ParameterList& plist) : - EvaluatorSecondaryMonotypeCV(plist) -{ - // determine the domain - Key akey = my_keys_.front().first; - domain_ = Keys::getDomain(akey); - akey = Keys::getVarName(akey); - domain_snow_ = Keys::readDomainHint(plist_, domain_, "surface", "snow"); - auto tag = my_keys_.front().second; - - // my keys - // -- sources - // gather all dependencies, assign any default parameters - // would prefer to do this as simply as possible, - // because the function signatures will change soon -} - -// Required methods from EvaluatorSecondaryMonotypeCV -void -CanopyHydrologyEvaluator::Evaluate_(const State& S, - const std::vector& results) -{ - auto mesh = S.GetMesh(domain_); - auto tag = my_keys_.front().second; - - // collect dependencies - - // collect output vecs - - // call canopy_hydrology kernels - ELM::kokkos_canopy_hydrology(*S, atm_forcing->forc_PREC, dtime, time_plus_half_dt) - -} - -void -CanopyHydrologyEvaluator::EvaluatePartialDerivative_(const State& S, - const Key& wrt_key, const Tag& wrt_tag, - const std::vector& results) {//nope} - - -// custom EC used to set subfield names -void -CanopyHydrologyEvaluator::EnsureCompatibility_Structure_(State& S) -{ - S.GetRecordSetW(my_keys_.front().first).set_subfieldnames( - {"bare", "water", "snow"}); -} - - -void CanopyHydrologyEvaluator::EnsureCompatibility_ToDeps_(State& S) -{ - // new state! - if (land_cover_.size() == 0) - land_cover_ = getLandCover(S.ICList().sublist("land cover types"), - {"albedo_ground", "emissivity_ground"}); - - for (auto dep : dependencies_) { - auto& fac = S.Require(dep.first, dep.second); - if (Keys::getDomain(dep.first) == domain_snow_) { - fac.SetMesh(S.GetMesh(domain_snow_)) - ->SetGhosted() - ->AddComponent("cell", AmanziMesh::CELL, 1); - } else { - fac.SetMesh(S.GetMesh(domain_)) - ->SetGhosted() - ->AddComponent("cell", AmanziMesh::CELL, 1); - } - } -} - -} // namespace Relations -} // namespace ELMKernels diff --git a/src/pks/surface_balance/ELMKernels/canopy_hydrology_evaluator.hh b/src/pks/surface_balance/ELMKernels/canopy_hydrology_evaluator.hh deleted file mode 100644 index 8e09e2e1c..000000000 --- a/src/pks/surface_balance/ELMKernels/canopy_hydrology_evaluator.hh +++ /dev/null @@ -1,53 +0,0 @@ -#pragma once - -#include "Factory.hh" -#include "EvaluatorSecondaryMonotype.hh" -#include "LandCover.hh" - -namespace Amanzi { -namespace ELMKernels { - -class CanopyHydrologyEvaluator : public EvaluatorSecondaryMonotypeCV { - public: - explicit CanopyHydrologyEvaluator(Teuchos::ParameterList& plist); - CanopyHydrologyEvaluator(const CanopyHydrologyEvaluator& other) = default; - - virtual Teuchos::RCP Clone() const override { - return Teuchos::rcp(new CanopyHydrologyEvaluator(*this)); - } - - protected: - // custom EC used to set subfield names - virtual void EnsureCompatibility_Structure_(State& S) override; - - // custom EC used because deps have 1 component not 3 - virtual void EnsureCompatibility_ToDeps_(State& S) override; - - // Required methods from EvaluatorSecondaryMonotypeCV - virtual void Evaluate_(const State& S, - const std::vector& results) override; - - virtual void EvaluatePartialDerivative_(const State& S, - const Key& wrt_key, const Tag& wrt_tag, - const std::vector& results) override; - - protected: - Key domain_; - Key domain_snow_; - - Key albedo_key_, emissivity_key_; - Key snow_dens_key_; - Key unfrozen_fraction_key_; - - double a_water_, a_ice_; - double e_water_, e_ice_, e_snow_; - - // this is horrid, because this cannot yet live in state - // bring on new state! - - private: - static Utils::RegisteredFactory reg_; -}; - -} // namespace Relations -} // namespace ELMKernels diff --git a/src/pks/surface_balance/ELMKernels/interface_archive.hh b/src/pks/surface_balance/ELMKernels/interface_archive.hh deleted file mode 100644 index 9b0efe1ee..000000000 --- a/src/pks/surface_balance/ELMKernels/interface_archive.hh +++ /dev/null @@ -1,291 +0,0 @@ - -#ifndef ATS_CLM_INTERFACE_HH_ -#define ATS_CLM_INTERFACE_HH_ - -#include -#include - -#include "Epetra_MultiVector.h" - -#define NUM_LC_CLASSES 18 - -namespace ATS { -namespace ELMKernels { - - class ELMKernelDriver { - - std::shared_ptr ELMState_; - - // time-invariant, no spatial aspect - std::shared_ptr> snicar_data_; - std::shared_ptr> snw_rds_table_; - std::shared_ptr> pft_data_; - - // time variable, but no spatial aspect, unless domain very large - // holds 12 monthly values for aerosol dep - // interpolates input file using lat and lon - std::shared_ptr> aerosol_data_; - - // time and space variable - must be time-interpolated at every step - std::shared_ptr atm_forcing_; - std::shared_ptr> phen_data_; - - // derived from aerosol_data and snow physics - std::shared_ptr> aerosol_masses_; - std::shared_ptr> aerosol_concentrations_; - - // the path/to/file for phenology and atm forcing inputs - std::string fname_surfdata_; - std::string fname_forc_; - - // TODO better - // number of columns, cell per column - int ncols_{0}; - int cell_per_col{0}; - - // TODO should simplify this - // don't need 2D domain decomp with ATS - // will also likely have ATS do the file reads - DomainDecomposition<2> dd_; - - // for phenology input data - should hide this in phen_data - std::unordered_map host_phen_views_; - - ELM::Utils::Date start_; - - -// mapping between ATS and ELM -// if ATS list of local nodes for each rank isn't contiguously ordered -// use a std::map -// where map[0...ncols] = ATS node id -// so ELM_var(i) = ATS_var(map[i]) -// this should be fast to iterate? - -// SerialDenseVec? -// soil temp solve? -// mesh definition - pass in mesh_, or just vector of interfaces (zisoi) to get dz, zsoi, zisoi -// slope aspect - -// need: -// Land -// lat -// lon -// lat_r -// lon_r -// dewmx -// irrig_rate -// n_irrig_steps_left -// oldfflag -// veg_active -// do_capsnow -// topo_slope -// topo_std -// t_h2osfc -// tsoi -// altmax_indx -// altmax_lastyear_indx -// t10 -// t_veg -// dz | -// zsoi | - assigned together -// zisoi | -// h2osoi_ice | -// h2osoi_liq | - together? -// h2osoi_vol | - -// need to make phenology_data and helpers inot it's own call -// make AtmForc read start date from the files, instead of providing -// should make a decision on clock sync -// put coszen into functions, also max_dayl and dayl - - - - - public: - - int init_lsm(int ncols, const MPI_Comm& comm, std::string_view fsurfdata, - std::string_view fsnicar, std::string_view fforc, - std::string_view fparam, std::string_view faerosol, - std::string_view fsnowage) { - // assign filenames - fname_surfdata_{fsurfdata}; - fname_forc_{fforc}; - - ncols_ = ncols; - cell_per_col = nlevgrnd; // hardwired in various places within ELMKernels - // may break if not set at 15 - - // assert a bunch to make sure everything lines up - - // need to figure out mapping - // implement better way of defining local/global problem - // knowledge of global problem is currently required for - // pnetcdf usage - if ATS reads, we don't need it - dd_ = ELM::Utils::create_domain_decomposition_2D( - { ncols, ncols }, - { 1, 1 }, - { 0, 0 }); - dd_.comm = comm; - - // bloated state struct - will move some of this data - ELMState_ = std::make_shared(ncols); - - // time invariant constants - // without spatial awareness - snicar_data_ = std::make_shared>(); - snw_rds_table_ = std::make_shared>(); - pft_data_ = std::make_shared>(); - - // one year of aerosol deposition from nearest gridcell - // currently invariant after this read - aerosol_data_ = std::make_shared>(); - - ELM::initialize_kokkos_elm(*ELMState_, *snicar_data_, *snw_rds_table_, *pft_data_, - *aerosol_data_, dd_, fname_surfdata_, fname_param_, - fname_snicar_, fname_snowage_, fname_aerosol_); - return 0; - } - - - // initialization of state fields - // must deep copy all of these into kokkos - - // can also be provided on cell-by-cell basis, - // which will be necessary for restarts - int init_snow(double snow_depth, double snow_temp) { - return 0; - } - int init_snow(const Epetra_MultiVector& snow_depth) { - return 0; - } // snow temp needs to be taken care of here as well - - // should probably be run after snl calculated - int init_temperature(double const_temp) { - // deep copy to State view - assign(ELMState_->t_h2osfc, const_temp); - // copy constant value into subsurface only - // snow layers should be initialized as 0.0 - auto h_tsoi = Kokkos::create_mirror_view(ELMState_->t_soisno); - for (int i = 0; i < ncols_; ++i) { - for (int j = nlevgrnd; j < nlevsno + nlevgrnd; ++j) { - h_tsoi(i, j) = const_temp; - } - } - assign(ELMState_->t_soisno, h_tsoi); - return 0; - } - - // should make a decision on clock sync - int init_start_time(int year, int month, int day) { - start_ = ELM::Utils::Date(year, month, day); - return 0; - } - - // just pass in mesh_? or zisoi? - int init_mesh(const Epetra_MultiVector& dz) { - return 0; - } // should be constant with 1 column worth of dz, or zisoi - - - - - - - int set_energy_state(const Epetra_MultiVector& soil_temp) { - - // place soil_temp from ATS into t_soisno once mapping decided upon - return 0; - } - - int set_water_state(const Epetra_MultiVector& saturation, const Epetra_MultiVector& porosity) { - // place saturation from ATS into watsat once mapping decided upon - // h2osoi_liq, h2osoi_ice?, watsat - return 0; - } - - - int transfer_function_VG_to_CH() { - // we're ignoring pressure, so this probably isn't needed - return 0; - } - - - int setup_lsm() { - - int atm_nsteps = 101; // how many timesteps to read? - const auto fstart = ELM::Utils::Date(1985, 1, 1); // make atm_forcing read this - may already be done - - // need to interpolate forcing at every step and - // reread from file periodically - happens automatically when forcing step (atm_nsteps-1)/atm_nsteps has been consumed - atm_forcing_ = std::make_shared(fname_forc_, fstart, atm_nsteps, ncols); - - // phenology data manager - // make host mirrors - need to be persistent - phen_data_ = std::make_shared>(dd, ncells, 17); - host_phen_views_ = get_phen_host_views(*phen_data); - - // containers for aerosol deposition and concentration within snowpack layers - aerosol_masses_ = std::make_shared>(ncells); - aerosol_concentrations_ = std::make_shared>(ncells); - } - - - - int get_water_flux_diagnostic() { - // get some/all of the optional qflx outputs - // need to calc qflx_infiltration - // with qflx_top_soil, qflx_snomelt, etc - // from SoilHydrologyMod.F90 - // qflx_snow_grnd - // qflx_rain_grnd - // qflx_snow_melt - // qflx_evap_tot - // qflx_evap_soil - // qflx_evap_veg - // qflx_tran_veg - // qflx_irr - // irr_flag - return 0; - } - - int get_water_flux_total() { - // get the accumulated qflx outputs - // convert subsurf flux mm/s to 1/s -- mm to m 1e3 - // keep sfc flux in L/T - return 0; - } - - int get_energy_flux_diagnostic() { - // get some/all of the optional eflx outputs - //eflx_sh_grnd - //eflx_sh_snow - //eflx_sh_soil - //eflx_sh_h2osfc - return 0; - } - - int get_energy_flux_total() { - // get the accumulated eflx outputs - // eflx_lh_tot - // eflx_sh_tot - // eflx_lwrad_out - or eflx_lwrad_net?? - // eflx_soil - return 0; - } - - - int advance_time(double time, double dt) { - // call 2 functions - prepare init_timestep() and run_lsm_physics - return 0; - } - - - }; - - -} // namespace ATS -} // namespace ELMKernels - - - -#endif From e25ce0f40f2d80b62effa795ea5300f336fd99a0 Mon Sep 17 00:00:00 2001 From: Joe Beisman Date: Mon, 23 Oct 2023 20:28:18 -0400 Subject: [PATCH 7/7] ELMKernels surface balance PK compiles and links against ELMKernels library --- .../ELMKernels/surface_balance_ELMKernels.cc | 16 ++++++++++------ .../ELMKernels/surface_balance_ELMKernels.hh | 8 ++++---- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc index 0874c49ea..21c9ccf89 100644 --- a/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc +++ b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc @@ -10,12 +10,10 @@ #include "pk_helpers.hh" #include "surface_balance_ELMKernels.hh" -//#include "elm_kokkos_interface.hh" - #define NUM_LC_CLASSES 18 namespace Amanzi { -namespace LandPhysics { +namespace SurfaceBalance { SurfaceBalanceELMKernels::SurfaceBalanceELMKernels(Teuchos::ParameterList& pk_tree, const Teuchos::RCP& global_list, @@ -72,6 +70,11 @@ SurfaceBalanceELMKernels::SurfaceBalanceELMKernels(Teuchos::ParameterList& pk_tr // ELMKernels timestep dt_ = plist_->get("time step size [s]"); + + // dummy ncols for now + int ncol = 1; + // Set up the ELMKernels interface + elm_ = Teuchos::rcp(new ELM::ELMInterface(ncol)); } // main methods @@ -82,6 +85,10 @@ SurfaceBalanceELMKernels::Setup() PK_Physical_Default::Setup(); auto subsurf_mesh = S_->GetMesh(domain_ss_); + // Set up the ELMKernels interface + int ncols = subsurf_mesh->num_columns(); + elm_ = Teuchos::rcp(new ELM::ELMInterface(ncols)); + // requirements: primary variable // -- snow depth S_->Require(key_, tag_next_, name_) @@ -100,9 +107,6 @@ SurfaceBalanceELMKernels::Setup() // set requirements on dependencies SetupDependencies_(tag_next_); - - // Set up the ELMKernels interface - elm_interface_ = Teuchos::rcp(new Amanzi::LandPhysics::Kernels(engine_name, engine_inputfile)); } void diff --git a/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.hh b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.hh index 57eeac2e7..56bd24ba6 100644 --- a/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.hh +++ b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.hh @@ -15,10 +15,10 @@ #include "PK_Factory.hh" #include "pk_physical_bdf_default.hh" -#include "ats_elm_interface.hh" +#include "elm_kokkos_interface.hh" namespace Amanzi { -namespace LandPhysics { +namespace SurfaceBalance { class SurfaceBalanceELMKernels : public PK_Physical_Default { @@ -97,14 +97,14 @@ public: Key color_index_key_; Key pft_index_key_; - Teuchos::RCP elm_interface{nullptr}; + Teuchos::RCP elm_{nullptr}; private: // factory registration static RegisteredPKFactory reg_; }; -} // namespace LandPhysics +} // namespace SurfaceBalance } // namespace ATS #endif