diff --git a/src/executables/CMakeLists.txt b/src/executables/CMakeLists.txt index f1cdec8c48..eafe916e40 100644 --- a/src/executables/CMakeLists.txt +++ b/src/executables/CMakeLists.txt @@ -121,6 +121,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 @@ -135,6 +139,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 22f7afa5af..d736899501 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) @@ -52,6 +59,13 @@ if (ENABLE_CLM) ) endif() +set(ENABLE_ELM true) +if (ENABLE_ELM) + 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 @@ -97,10 +111,17 @@ if (ENABLE_CLM) ) endif() +if (ENABLE_ELM) + list(APPEND ats_surface_balance_inc_files + ELMKernels/surface_balance_ELMKernels.hh + ) +endif() + set(ats_surface_balance_link_libs ${Teuchos_LIBRARIES} ${Epetra_LIBRARIES} ${CLM_LIBRARIES} + ${ELM_LIBRARIES} error_handling atk mesh @@ -147,6 +168,13 @@ if (ENABLE_CLM) ) endif() +if (ENABLE_ELM) + 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 0000000000..21c9ccf899 --- /dev/null +++ b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.cc @@ -0,0 +1,462 @@ +/* + 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" + +#define NUM_LC_CLASSES 18 + +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]"); + + // dummy ncols for now + int ncol = 1; + // Set up the ELMKernels interface + elm_ = Teuchos::rcp(new ELM::ELMInterface(ncol)); +} + +// main methods +// -- Setup data. +void +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_) + .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_); +} + +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) +{ + elm_->setup(); + // 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_; + + 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); + 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 0000000000..56bd24ba6a --- /dev/null +++ b/src/pks/surface_balance/ELMKernels/surface_balance_ELMKernels.hh @@ -0,0 +1,110 @@ +/* -*- 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" +#include "elm_kokkos_interface.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_; + + Teuchos::RCP elm_{nullptr}; + + private: + // factory registration + static RegisteredPKFactory reg_; +}; + +} // namespace SurfaceBalance +} // namespace ATS + +#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 0000000000..0238910ab2 --- /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