diff --git a/.gitignore b/.gitignore index d4d11ca8..ddcbafa9 100755 --- a/.gitignore +++ b/.gitignore @@ -8,6 +8,8 @@ third_party/snowpack/plugins/SnowpackIO.cc *.bak CRHMr.log .Rhistory +**/*.swp +**/*.swo /.idea/ /.spack-env/ @@ -27,4 +29,4 @@ CRHMr.log /docs/api/* /docs/_build/* /docs/doxygen/* -/docs/_ext/__pycache__ \ No newline at end of file +/docs/_ext/__pycache__ diff --git a/CMakeLists.txt b/CMakeLists.txt index d73a4be4..c7c3a266 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required (VERSION 3.30) +cmake_minimum_required (VERSION 3.27) project (CHM-project C CXX Fortran) # MACOSX_RPATH is enabled by default. @@ -16,7 +16,7 @@ cmake_policy(SET CMP0144 NEW) # The NEW behavior is for find_package(Boost) to search for the upstream BoostConfig.cmake # https://cmake.org/cmake/help/latest/policy/CMP0167.html -cmake_policy(SET CMP0167 NEW) +#cmake_policy(SET CMP0167 NEW) # Needed in addition to the target cxx standard set later set(CMAKE_CXX_STANDARD 20) @@ -337,4 +337,4 @@ add_subdirectory(src) if(BUILD_DOCS) add_subdirectory(docs) -endif() \ No newline at end of file +endif() diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b0571273..4c6be3de 100755 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -42,6 +42,7 @@ set(MODULE_SRCS modules/Marsh_shading_iswr.cpp modules/Sicart_ilwr.cpp modules/PenmanMonteith_evaporation.cpp + modules/Evapotranspiration_All.cpp modules/interp_met/Thornton_p.cpp modules/Walcek_cloud.cpp modules/Harder_precip_phase.cpp @@ -50,7 +51,9 @@ set(MODULE_SRCS modules/Richard_albedo.cpp modules/snowpack.cpp modules/Gray_inf.cpp - modules/Simple_Canopy.cpp + modules/Infil_All.cpp + modules/soil_module.cpp + modules/Simple_Canopy.cpp modules/scale_wind_vert.cpp modules/sub_grid.cpp modules/snow_slide.cpp @@ -68,6 +71,21 @@ set(MODULE_SRCS CACHE INTERNAL "" FORCE ) +set(MODEL_SRCS + modules/EvapotranspirationModels/evapbase.cpp + modules/EvapotranspirationModels/PenmanMonteith.cpp + modules/EvapotranspirationModels/PriestleyTaylor.cpp + + modules/soil_submodules/soil_two_layer.cpp + modules/soil_submodules/soil_ET.cpp + modules/soil_submodules/soil_classes.cpp + + modules/soil_submodules/K_estimate.cpp + + modules/soil_submodules/XG_algorithm.cpp + + modules/infil_submodules/Crack.cpp + CACHE INTERNAL "" FORCE ) set(CHM_SRCS #main.cpp needs to be added below so we can re use CHM_SRCS in the gtest build @@ -77,6 +95,7 @@ set(CHM_SRCS metdata.cpp physics/Atmosphere.cpp + physics/Soil.cpp mesh/triangulation.cpp # mesh/ugrid.cpp @@ -109,7 +128,10 @@ set (HEADER_FILES modules/interp_met modules/snobal modules/testing - math + modules/EvapotranspirationModels/ + modules/soil_submodules/ + modules/infil_submodules/ + math libmaw interpolation timeseries @@ -227,7 +249,8 @@ add_executable( main.cpp # this needs to be here so we can reuse CHM_SRCS in the gtest build. but since it links it's own main, we cannot have this one. ${CHM_SRCS} ${FILTER_SRCS} - ${MODULE_SRCS} + ${MODEL_SRCS} + ${MODULE_SRCS} ) @@ -394,26 +417,31 @@ endif() if (BUILD_TESTS) message(STATUS "Tests enabled. Run with make check") - add_subdirectory (tests/googletest/googletest) + add_subdirectory (tests/googletest) set(TEST_SRCS - tests/test_station.cpp - tests/test_interpolation.cpp - tests/test_timeseries.cpp - tests/test_core.cpp - tests/test_variablestorage.cpp - tests/test_metdata.cpp - tests/test_netcdf.cpp + #tests/test_station.cpp + #tests/test_interpolation.cpp + #tests/test_timeseries.cpp + #tests/test_core.cpp + #tests/test_variablestorage.cpp + #tests/test_metdata.cpp + #tests/test_netcdf.cpp # test_mesh.cpp - tests/test_regexptokenizer.cpp + #tests/test_regexptokenizer.cpp # test_daily.cpp - tests/test_triangulation.cpp - tests/main.cpp -) + #tests/test_triangulation.cpp + tests/main.cpp + tests/submoduletests/test_crack.cpp + tests/submoduletests/test_ayers.cpp + tests/submoduletests/test_CSVreader.cpp + tests/submoduletests/test_daily_accumulator.cpp + ) set( GTEST_LINK gtest + gmock ) @@ -421,7 +449,7 @@ if (BUILD_TESTS) add_executable( runUnitTests ${CHM_SRCS} - ${FILTER_SRCS} + ${MODEL_SRCS} ${MODULE_SRCS} ${TEST_SRCS} ) diff --git a/src/logger.hpp b/src/logger.hpp index 5434dc17..eae135a9 100644 --- a/src/logger.hpp +++ b/src/logger.hpp @@ -41,9 +41,9 @@ }; //#ifdef USE_MPI -//#define MPI_RANK_DBG(RANK) if(_comm_world.rank() == RANK) \ -// {\ -// LOG_DEBUG << "\n\n-----------------------------\n\nI am PID " << getpid() <<", attach within 45s\n\n-----------------------------\n\n"; \ -// sleep(45);\ -// } -//#endif \ No newline at end of file +//#define MPI_RANK_DBG(RANK) if(_comm_world.rank() == RANK) +// { +// LOG_DEBUG << "\n\n-----------------------------\n\nI am PID " << getpid() <<", attach within 45s\n\n-----------------------------\n\n"; +// sleep(45); +// t } +//#endif diff --git a/src/mesh/triangulation.hpp b/src/mesh/triangulation.hpp index 705788b5..4164c645 100644 --- a/src/mesh/triangulation.hpp +++ b/src/mesh/triangulation.hpp @@ -324,6 +324,32 @@ class face : public Fb */ double veg_attribute(const std::string &variable); + /** + * Checks if the face has soil by loading the soil_storage_max parameter + * and then checking if it is greater than zero. + */ + bool has_soil(); + + /** + * Retrieves a soil attribute for the face. With a templated return type to allow for + * double, int or std::string types. Distributed variables will only be double/int. But + * look-up tables can point to strings or another type. + * */ + template + T soil_attribute(const std::string variable); + + /** + * This version of the function specifies that you need a look-up table. + * It is essentially identical to veg_attribute except it does not include the first if statement + * and it allows you to specify the look up map category name. For soils, it will always be soils. + * The obvious problem si that with veg_attribute you can keep a module unchanged and only change the + * config with or without a look-up. This version requires changing the module to stop using a look up. + * Consider changing in the future. + * This is only found in the one argument version of this function + * */ + template + T soil_attribute(const std::string variable,const std::string category); + /** * Sets the vector for the given variable. * Does not support timeseries output. @@ -1737,6 +1763,61 @@ double face::veg_attribute(const std::string &variable) return result; }; +template < class Gt, class Fb > +bool face::has_soil() +{ + if (has_parameter("soil_storage_max"_s) ) + { + double result = parameter("soil_storage_max"_s); + if (result >= 0.0) + return true; + } + + return false; +} + +template < class Gt, class Fb > +template< typename T > +T face::soil_attribute(const std::string variable) +{ + + T result{}; + // first see if we have a distributed map of this parameter + if(has_parameter(variable)) + result = parameter(variable); + else + { + CHM_THROW_EXCEPTION(module_error, "Parameter " + variable +" does not exist."); + } + return result; +}; + +template < class Gt, class Fb > +template< typename T > +T face::soil_attribute(const std::string variable,const std::string category) +{ + T result{}; + if (has_parameter(variable)) + { + int LC = parameter(variable); + auto param = _domain->_global->parameters; + try + { + result = param.get(category + "." + variable + "." + std::to_string(LC)); + } + catch(const boost::property_tree::ptree_bad_path& e) + { + CHM_THROW_EXCEPTION(module_error, "Parameter " + variable + " in category " + category +" does not exist."); + } + } + else + { + CHM_THROW_EXCEPTION(module_error, "Parameter " + variable + " in category " + category +" does not exist."); + } + return result; +}; + + template < class Gt, class Fb> Vector_3 face::face_vector(const std::string& variable) { diff --git a/src/modules/EvapotranspirationModels/PenmanMonteith.cpp b/src/modules/EvapotranspirationModels/PenmanMonteith.cpp new file mode 100644 index 00000000..81f1fe38 --- /dev/null +++ b/src/modules/EvapotranspirationModels/PenmanMonteith.cpp @@ -0,0 +1,114 @@ +#include "PenmanMonteith.hpp" + + +PenmanMonteith::PenmanMonteith(double& LAI, double& LAImax, double& veg_Ht, double& wind_height, double& stomatal_res_min, double& soil_d, double& F_to_g, const double Cp, const double K, const double tension, const double pore_sz, const double theta_pwp, const double phi) + : leaf_area_index(LAI), leaf_area_index_max(LAImax), Veg_height(veg_Ht), wind_measurement_height(wind_height), stomatal_resistance_min(stomatal_res_min), soil_depth(soil_d), Frac_to_ground(F_to_g), heat_capacity_air(Cp), kappa(K), air_entry_tension(tension), pore_size_dist(pore_sz), wilt_point(theta_pwp), porosity(phi) +{ + if (leaf_area_index == 0) + has_vegetation = false; + else + has_vegetation = true; + + + // TODO other checks on values might be good +} + +PenmanMonteith::~PenmanMonteith() +{ + //Do nothing +} + +void PenmanMonteith::CalcHeights() +{ + Z0 = Veg_height/7.6; + d = Veg_height*0.67; +} + +double PenmanMonteith::CalcAeroResistance(const PM_vars& var) +{ + if (wind_measurement_height - d > 0) + { + return pow( log((wind_measurement_height - d)/Z0),2) / + (pow(kappa,2) * var.wind_speed); + } + else + { + return 0; // I don't know if this is right, but it is at least... safe. + } +} + +double PenmanMonteith::CalcStomatalResistance(const PM_vars& var) +{ + double rcstar = stomatal_resistance_min; + + // In CRHM, the below calculation is an option, for now just use the minimum option. + if (has_vegetation) + { + + double LAI = Veg_height/2.0*leaf_area_index_max; //TODO ad hoc for test + rcstar = stomatal_resistance_min * leaf_area_index_max / LAI; + // rcstar = stomatal_resistance_min * leaf_area_index_max / leaf_area_index; TODO commented for test + } + // TODO check units. for example, short_wave_in - 1.5 is suspect + double f1 = 1.0; + if (var.short_wave_in > 0.0) + f1 = std::max(1.0, 500.0/var.short_wave_in - 1.5); +//max (1.0, 500.0/(var.short_wave_in - 1.5)); + + double f2 = std::max(1.0, 2.0 * (var.saturated_vapour_pressure - var.vapour_pressure) ); +// (1.0, 2.0 * (var.saturated_vapour_pressure - var.vapour_pressure) ); + + double p = air_entry_tension * pow(porosity / (var.soil_storage/soil_depth/1000.0 + wilt_point), pore_size_dist); + double f3 = std::max(1.0, p/40.0); + + double f4 = 1.0; + if (var.t < 5.0 || var.t > 40.0) + f4 = 5000/stomatal_resistance_min; + + if (var.short_wave_in <= 0) + return 5000; + else + { + return std::min(rcstar * f1 * f2 * f3 * f4, 5000.0); + } + +} + +void PenmanMonteith::CalcEvapT(var_base& basevar, model_output& output) +{ + + const PM_vars & var = static_cast(basevar); + PM_output &out = static_cast(output); + + double Q = var.all_wave_net * (1 - Frac_to_ground); + + if (IsFirstRun) + { + CalcHeights(); + IsFirstRun = false; + } + + double aero_resistance = CalcAeroResistance(var); + double stomatal_resistance = CalcStomatalResistance(var); + + double radiation = delta(var.t) * Q * 24; //24 --> hours/day, should be generalized or encoded in Q in a future version + // Actual units are mm/day + + double mass = AirDensity(var.t,var.vapour_pressure,var.P_atm) * heat_capacity_air * + (var.saturated_vapour_pressure - var.vapour_pressure)/ ( aero_resistance * lambda(var.t) * water_density ); + mass *= m_per_s_to_mm_per_day; + + out.ET = (radiation + mass) / ( delta(var.t) + gamma(var.P_atm, var.t, heat_capacity_air) * ( 1 + stomatal_resistance / aero_resistance )); + + out.ET *= 1.0 / 24.0; // As with the radiation variable, assumed the timestep is hours here. TODO for the future. + out.stomatal_resistance = stomatal_resistance; +} + +// TODO make delta, gamma, density fucntions + +double PenmanMonteith::AirDensity(double& t, double& ea, double& Pa) // atmospheric density (kg/m^3) +{ + const double R0 = 2870; + return (1E4*Pa /(R0*( 273.15 + t))*(1.0 - 0.379*(ea/Pa))); // +} + diff --git a/src/modules/EvapotranspirationModels/PenmanMonteith.hpp b/src/modules/EvapotranspirationModels/PenmanMonteith.hpp new file mode 100644 index 00000000..cddcaf3d --- /dev/null +++ b/src/modules/EvapotranspirationModels/PenmanMonteith.hpp @@ -0,0 +1,93 @@ +#pragma once + +#include "evapbase.hpp" +#include +#include +// Implementation: +// include the header +// std::unique_ptr mymodel = std::make_unqiue(); +// std::unique_ptr param = std::make_unique(); TODO may not need to make param a pointer because it doesnt need to be released on each timestep (in fact, it shouldn't) +// std::unique_ptr var = std::make_unique(); +// +// var must be remade because values are gathered from other modules. +// set param and var members. param never changes, var changes every timestep and should be released. + +struct PM_vars; +struct PM_output; + +class PenmanMonteith : public evapT_base +{ +public: + + PenmanMonteith(double& LAI, double& LAImax, double& veg_Ht, double& wind_height, + double& stomatal_res_min, double& soil_d, double& F_to_g, const double Cp, + const double K, const double tension, const double pore_sz, + const double theta_pwp, const double phi); + + + ~PenmanMonteith(void) override; // Deconstructor + + + void CalcEvapT(var_base& vars, model_output& output) override; + + // TODO These should be references + // double Veg_height; + // double Veg_height_max; + + double& leaf_area_index; + //double LAImin; + //double seasonal_growth; + double& leaf_area_index_max; + double& Veg_height; + double& wind_measurement_height; // This one might be uniform... + double& stomatal_resistance_min; // Also might be domain wide + double& soil_depth; + double& Frac_to_ground; + const double heat_capacity_air; + const double kappa; // also might be domain wide + const double air_entry_tension; + const double pore_size_dist; + const double wilt_point; + const double porosity; +private: + + // dont delete + void CalcHeights(void); + double CalcAeroResistance(const PM_vars& var); + double CalcStomatalResistance(const PM_vars& var); + double AirDensity(double& t, double& ea, double& Pa); + + double Z0; + double d; + bool has_vegetation; + bool IsFirstRun = true; + static constexpr double water_density = 1000; //kg/m^3 + static constexpr long m_per_s_to_mm_per_day = 1000 * 86400; +}; + + +struct PM_vars : public var_base +{ + double& wind_speed; + double& short_wave_in; + double& all_wave_net; + double& t; + double& soil_storage; + double& vapour_pressure; + double& saturated_vapour_pressure; + double& P_atm; + + PM_vars(double& U, double& Qsw, double& Qnet, double& temp, double& soil, double& ea, double& ea_star, double& P) + : wind_speed(U), short_wave_in(Qsw), all_wave_net(Qnet), t(temp), soil_storage(soil), vapour_pressure(ea), saturated_vapour_pressure(ea_star), P_atm(P) {}; // TODO Add a constructor, make these references. + // maybe not needed, remember it is normal to make a local copy of values on faces, which these are. +}; + + + // Make this a variable passed to evap? double ShortWave_in; + // Same as above double t; + + +struct PM_output : public model_output +{ + double stomatal_resistance; +}; diff --git a/src/modules/EvapotranspirationModels/PriestleyTaylor.cpp b/src/modules/EvapotranspirationModels/PriestleyTaylor.cpp new file mode 100644 index 00000000..1982c3e3 --- /dev/null +++ b/src/modules/EvapotranspirationModels/PriestleyTaylor.cpp @@ -0,0 +1,20 @@ +#include "PriestleyTaylor.hpp" + +PriestleyTaylor::PriestleyTaylor(const double& alpha_const, const double& Cp) : alpha(alpha_const), heat_capacity_air(Cp) +{ + +} + +PriestleyTaylor::~PriestleyTaylor() +{ + // Do nothing +} + +void PriestleyTaylor::CalcEvapT(var_base& basevar, model_output& output) +{ + const PT_vars& var = static_cast(basevar); + + double Q = var.all_wave_net * (1 - Frac_to_ground); + + output.ET = alpha * delta(var.air_temperature) * Q / (delta(var.air_temperature) + gamma(var.P_atm,var.air_temperature,heat_capacity_air) ); +} diff --git a/src/modules/EvapotranspirationModels/PriestleyTaylor.hpp b/src/modules/EvapotranspirationModels/PriestleyTaylor.hpp new file mode 100644 index 00000000..8af77242 --- /dev/null +++ b/src/modules/EvapotranspirationModels/PriestleyTaylor.hpp @@ -0,0 +1,49 @@ +#pragma once + +#include "evapbase.hpp" +#include + +// Implementation: +// include the header +// std::unique_ptr mymodel = std::make_unqiue(); +// std::unique_ptr param = std::make_unique(); TODO may not need to make param a pointer because it doesnt need to be released on each timestep (in fact, it shouldn't) +// std::unique_ptr var = std::make_unique(); +// +// var must be remade because values are gathered from other modules. +// set param and var members. param never changes, var changes every timestep and should be released. + + +class PriestleyTaylor : public evapT_base +{ +public: + + PriestleyTaylor(const double& alpha_const, const double& Cp); + + ~PriestleyTaylor(void) override; // Deconstructor + + + void CalcEvapT(var_base& vars, model_output& output) override; + + double Frac_to_ground; + const double& alpha; + const double heat_capacity_air; +private: + + +}; + + +struct PT_vars : public var_base +{ + double& all_wave_net; + double& P_atm; + double& air_temperature; + + PT_vars(double& Q, double& P, double& t) : all_wave_net(Q), P_atm(P), air_temperature(t) {}; +}; + + // Make this a variable passed to evap? double ShortWave_in; + // Same as above double t; + + + diff --git a/src/modules/EvapotranspirationModels/evapbase.cpp b/src/modules/EvapotranspirationModels/evapbase.cpp new file mode 100644 index 00000000..5b233ded --- /dev/null +++ b/src/modules/EvapotranspirationModels/evapbase.cpp @@ -0,0 +1,21 @@ +#include "EvapotranspirationModels/evapbase.hpp" + +double evapT_base::delta(const double& t) // Slope of sat vap p vs t, kPa/DEGREE_CELSIUS +{ + if (t > 0.0) + return(2504.0*exp(17.27 * t/(t+237.3)) / pow(t+237.3,2)); + else + return(3549.0*exp( 21.88 * t/(t+265.5)) / pow(t+265.5,2)); +} + +double evapT_base::lambda(const double& t) // Latent heat of vaporization (J/kg) +{ + // Equation 7-8 Dingman (2002) Second Edition + return (2.501 - 0.002361 * t) * 1e6; // original is MegaJoules/kg, 1e6 returns it to joules/kg +} + +double evapT_base::gamma(const double& Pa, const double& t, const double& c_a) // Psychrometric constant (kPa/DEGREE_CELSIUS) +{ + // Equation 7-13 Dingman Second Edition 2002 + return c_a * Pa / (0.622 * lambda(t)); // lambda (J/kg) +} diff --git a/src/modules/EvapotranspirationModels/evapbase.hpp b/src/modules/EvapotranspirationModels/evapbase.hpp new file mode 100644 index 00000000..4a65e2f2 --- /dev/null +++ b/src/modules/EvapotranspirationModels/evapbase.hpp @@ -0,0 +1,28 @@ +#pragma once + +#include +#include "logger.hpp" + +struct var_base +{ + +}; + +struct model_output +{ + double ET; + double rc; +}; + +class evapT_base +{ +public: + virtual ~evapT_base() = default; + + virtual void CalcEvapT(var_base& basevar, model_output& output) = 0; + + double delta(const double& t); + double gamma(const double& P_atm, const double& t, const double& c_a); + double lambda(const double& t); +}; + diff --git a/src/modules/Evapotranspiration_All.cpp b/src/modules/Evapotranspiration_All.cpp new file mode 100644 index 00000000..c8ed65b9 --- /dev/null +++ b/src/modules/Evapotranspiration_All.cpp @@ -0,0 +1,177 @@ +// +// Canadian Hydrological Model - The Canadian Hydrological Model (CHM) is a novel +// modular unstructured mesh based approach for hydrological modelling +// Copyright (C) 2018 Christopher Marsh +// +// This file is part of Canadian Hydrological Model. +// +// Canadian Hydrological Model is free software: you can redistribute it and/or +// modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// Canadian Hydrological Model is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with Canadian Hydrological Model. If not, see +// . +// + + +#include "Evapotranspiration_All.hpp" + +REGISTER_MODULE_CPP(Evapotranspiration_All); + +Evapotranspiration_All::Evapotranspiration_All(config_file cfg) + :module_base("Evapotranspiration_All", parallel::data, cfg) +{ + // TODO Constructor is not properly editted with all new inputs (see set vars function at the end) + depends("iswr"); + depends("netall"); + depends("P_atm"); + depends("ea"); + depends("t"); + depends("U_2m_above_srf"); // + depends("soil_storage"); // but is how albedo is used in CHM as of Sept, 2024 + + provides("ET"); + provides("stomatal_resistance"); + +} + +void Evapotranspiration_All::init(mesh& domain) +{ + alpha = cfg.get("alpha_PriestleyTaylor",1.26); + wind_height = cfg.get("wind_measurement_height",2.0); + stomatal_resistance_min = cfg.get("stomatal_resistance_min",62.0); + Frac_to_ground = cfg.get("Frac_to_ground",1.0); // iswr_subcanopy exists. + + SoilDataObj = std::make_unique(); + + for (size_t i = 0; i < domain->size_faces(); i++) + { + auto face = domain->face(i); + auto& d = face->make_module_data(ID); + + // Consider if an if statement is necessary. + + + d.soil_depth = face->soil_attribute("soil_depth"_s); + if (face->has_vegetation()) + { + d.LAI = face->veg_attribute("LAI"); + d.LAImax = face->veg_attribute("LAImax"); + d.vegetation_height = face->veg_attribute("CanopyHeight"); + } + else + { + d.LAI = 0.0; + d.LAImax = 0.0; + d.vegetation_height = 0.0; + } + // TODO Consider if we can have a single model object per triangle. + // Polymorpish would let this work well + // It wouldn't work if PT model is used for normal soils when they are saturated. + if (is_water(face)) + { + init_PriestleyTaylor(d,alpha); + } + // TODO add wetlands + else + { + init_PenmanMonteith(d, face, wind_height, stomatal_resistance_min, Frac_to_ground); + } + + } + +} + +void Evapotranspiration_All::run(mesh_elem& face) +{ + auto& d = face->get_module_data(ID); + + if (is_water(face)) + { + // Do PriestlyTaylor + PT_vars my_PT_vars = set_PriestleyTaylor_vars(face); + model_output output; + d.MyPriestleyTaylor->CalcEvapT(my_PT_vars,output); + + (*face)["ET"_s] = output.ET; + } + // TODO wetlands or saturated soils, need input from soils module. + // else if (If wetlands) + // { + // Also Do PriestlyTaylor + // } + else + { + // All members of PM_vars are references and must be set at initialization + // Therefore we need a copy of SVP to reference + // t is made its own copy to avoid dereferencing face for "t" twice + + double t = (*face)["t"_s]; + double SVP = Atmosphere::saturatedVapourPressure(t+273.15)/1000; // units of kelvin expected + double VP = (*face)["ea"_s]; + PM_vars my_PM_vars = set_PenmanMonteith_vars(face,t,SVP,VP); + PM_output output; + d.MyPenmanMonteith->CalcEvapT(my_PM_vars,output); + + (*face)["stomatal_resistance"_s] = output.stomatal_resistance; + + (*face)["ET"_s] = output.ET; + } + + // TODO total_ET, as well as PT ET and PM ET as separate. +} + +Evapotranspiration_All::~Evapotranspiration_All() +{ + +} + +void Evapotranspiration_All::init_PriestleyTaylor(Evapotranspiration_All::data& d,double& alpha) +{ + d.MyPriestleyTaylor = std::make_unique(alpha,Atmosphere::Cp); +} + +void Evapotranspiration_All::init_PenmanMonteith(Evapotranspiration_All::data& d,mesh_elem& face, double& wind_height, double& stomatal_resistance_min, double& Frac_to_ground) +{ + + const std::string soil_type = face->soil_attribute("soil_type"_s,"soils"); + + const double Cp = Atmosphere::Cp; + const double kappa = Atmosphere::kappa; + const double air_entry_tension = SoilDataObj->air_entry_tension(soil_type); + const double pore_size_dist = SoilDataObj->pore_size_dist(soil_type); + const double wilt_point = SoilDataObj->wilt_point(soil_type); + const double porosity = SoilDataObj->porosity(soil_type); + + double soil_storage_max = face->soil_attribute("soil_storage_max"_s); + // Leaf area index is not used if no vegetation, but LAI and LAImax are references in the PenmanMonteith model, therefore values are needed for initialization. It is ok if these values go out of scope as long as there is no vegetation. + + d.MyPenmanMonteith = std::make_unique(d.LAI, d.LAImax, d.vegetation_height, wind_height, + stomatal_resistance_min, d.soil_depth, Frac_to_ground, Cp, kappa, + air_entry_tension, pore_size_dist, wilt_point, porosity); + +} + +PM_vars Evapotranspiration_All::set_PenmanMonteith_vars(mesh_elem& face,double& t, double& saturated_vapour_pressure,double& vapour_pressure) +{ + PM_vars vars((*face)["U_2m_above_srf"_s],(*face)["iswr"_s],(*face)["netall"_s],t,(*face)["soil_storage"_s],vapour_pressure,saturated_vapour_pressure,(*face)["P_atm"_s]); + + return vars; +} + +PT_vars Evapotranspiration_All::set_PriestleyTaylor_vars(mesh_elem& face) +{ + // TODO P_atm, is a state variable and should have a copy local to the run function + // because it is calculated from the Atmosphere namespace, not done yet + PT_vars vars((*face)["netall"_s],(*face)["P_atm"_s],(*face)["t"_s]); + + return vars; +} diff --git a/src/modules/Evapotranspiration_All.hpp b/src/modules/Evapotranspiration_All.hpp new file mode 100644 index 00000000..1fe932a2 --- /dev/null +++ b/src/modules/Evapotranspiration_All.hpp @@ -0,0 +1,105 @@ +// +// Canadian Hydrological Model - The Canadian Hydrological Model (CHM) is a novel +// modular unstructured mesh based approach for hydrological modelling +// Copyright (C) 2018 Christopher Marsh +// +// This file is part of Canadian Hydrological Model. +// +// Canadian Hydrological Model is free software: you can redistribute it and/or +// modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// Canadian Hydrological Model is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with Canadian Hydrological Model. If not, see +// . +// + +#pragma once + +#include "logger.hpp" +#include "triangulation.hpp" +#include "module_base.hpp" +#include "Atmosphere.h" +#include "Soil.h" +#include +#include +#include +#include +#define _USE_MATH_DEFINES +#include +#include "EvapotranspirationModels/evapbase.hpp" +#include "EvapotranspirationModels/PenmanMonteith.hpp" +#include "EvapotranspirationModels/PriestleyTaylor.hpp" + + +/** + * \ingroup modules exp evap + * @{ + * \class Evapotranspiration_All + + * Calculates evapo-transpiration via Penman-Monteith. + * + * Not currently maintained. + * + * Depends: + * - Incoming shortwave radaition "iswr" [\f$ W \cdot m^{-2} \f$ ] + * - Incoming longwave radiation "ilwr" [\f$ W \cdot m^{-2} \f$ ] + * - Relative humidity "rh" [%] + * - Windspeed at 2m "U_2m_above_srf" [\f$ m \cdot s^{-1} \f$ ] + * + * Provides: + * - Evapotranspiration "ET" [\f$ mm \cdot dt^{-1} \f$ ] + * + * @} + */ +class Evapotranspiration_All : public module_base +{ +REGISTER_MODULE_HPP(Evapotranspiration_All); +public: + Evapotranspiration_All(config_file cfg); + ~Evapotranspiration_All(); + void init(mesh& domain); + void run(mesh_elem& face); + + class data : public face_info + { + public: + std::unique_ptr MyPenmanMonteith; + std::unique_ptr MyPriestleyTaylor; + + double LAI; + double LAImax; + double vegetation_height; + double soil_depth; + }; + +private: + + std::unique_ptr SoilDataObj; + + // PriestleyTaylor + double alpha; + // TODO add PT methods here + + void init_PriestleyTaylor(Evapotranspiration_All::data& d,double& alpha); + PT_vars set_PriestleyTaylor_vars(mesh_elem& face); + + // PenmanMonteith + double wind_height; + double stomatal_resistance_min; + double Frac_to_ground; + + void init_PenmanMonteith(Evapotranspiration_All::data& d, mesh_elem& face, double& wind_height, + double& stomatal_resistance_min, double& Frac_to_ground); + + PM_vars set_PenmanMonteith_vars(mesh_elem& face, double& t, + double& saturated_vapour_pressure, double& vapour_pressure); + +}; diff --git a/src/modules/Infil_All.cpp b/src/modules/Infil_All.cpp new file mode 100644 index 00000000..82aaebaf --- /dev/null +++ b/src/modules/Infil_All.cpp @@ -0,0 +1,391 @@ +// +// Canadian Hydrological Model - The Canadian Hydrological Model (CHM) is a novel +// modular unstructured mesh based approach for hydrological modelling +// Copyright (C) 2018 Christopher Marsh +// +// This file is part of Canadian Hydrological Model. +// +// Canadian Hydrological Model is free software: you can redistribute it and/or +// modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// Canadian Hydrological Model is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with Canadian Hydrological Model. If not, see +// . +// + +#include "Infil_All.hpp" +REGISTER_MODULE_CPP(Infil_All); + +Infil_All::Infil_All(config_file cfg) : module_base("Infil_All", parallel::data, cfg) +{ + + depends("swe"); + depends("snowmelt_int"); + depends("rainfall_int"); // NEW + depends("soil_storage_at_freeze"); // NEW, depends on Volumetric model, equivalent to fallstat in crhm + depends("soil_storage"); + depends("t"); + + provides("inf"); + provides("total_inf"); + provides("snowinf"); // NEW + provides("total_snowinf"); // NEW + provides("total_excess"); + provides("total_meltexcess"); // NEW + provides("runoff"); + provides("melt_runoff"); // NEW + provides("total_rain_on_snow"); // NEW + provides("rain_on_snow"); // NEW + provides("frozen"); + provides("major_melt_count"); +} + +Infil_All::~Infil_All() +{ + +} + +void Infil_All::init(mesh& domain) +{ + //store all of snobals global variables from this timestep to be used as ICs for the next timestep +#pragma omp parallel for + for (size_t i = 0; i < domain->size_faces(); i++) + { + auto face = domain->face(i); + auto& d = face->make_module_data(ID); + + // Triangle Totals + d.total_inf = 0.; + d.total_snowinf = 0.; // NEW + d.total_excess = 0.; + d.total_meltexcess = 0.; // NEW + d.total_rain_on_snow = 0.; // NEW + + // Triangle variables + d.crack_model_status.init(); + + d.soil_storage = face->soil_attribute("soil_storage"); + d.last_day = 0; + + // Model Parameters + infDays = cfg.get("max_inf_days",6); + min_swe_to_freeze = cfg.get("min_swe_to_freeze",25); + major = cfg.get("major",5); + AllowPriorInf = cfg.get("AllowPriorInf",true); + thaw_type = cfg.get("thaw_type",0); // Default is Ayers + + if (thaw_type == AYERS) + { + d.texture = face->soil_attribute("soil_texture","soils"); + d.ground_cover = face->soil_attribute("soil_groundcover","soils"); + } + else if (thaw_type == GREENAMPT) + { + d.soil_type = face->soil_attribute("soil_type","soils"); + d.ksaturated = SoilDataObj.saturated_conductivity(d.soil_type); + } + lenstemp = cfg.get("temperature_ice_lens",-10.0); + + + d.soil_storage_max = face->parameter("soil_storage_max"_s); + + + } +} +void Infil_All::run(mesh_elem &face) +{ + // TODO if its water it should probably take all rain as "infil", there will be no snowmelt... sorta, snow melting and leaking into the water under the ice in spring?? + if(is_water(face)) + { + set_all_nan_on_skip(face); + return; + } + + auto& d = face->get_module_data(ID); + + auto id = face->cell_local_id; + + // CRHM does total infil for snow and total infil separately, wonder if I should do this + double runoff = 0.0; + double melt_runoff = 0.0; + double inf = 0.0; + double snowinf = 0.0; + double rain_on_snow = 0.0; + + double snowmelt = (*face)["snowmelt_int"_s]; + double rainfall = (*face)["rainfall_int"_s]; // NEW + double swe = (*face)["swe"_s]; + double soil_storage_at_freeze = (*face)["soil_storage_at_freeze"_s]; + double airtemp = (*face)["t"_s]; + + if (thaw_type == GREENAMPT) + d.soil_storage = (*face)["soil_storage"_s]; + + + + if (swe > min_swe_to_freeze && !d.crack_model_status.frozen && is_new_day()) + { + d.crack_model_status.begin_freeze(); + d.crack_model_status.end_freeze_tomorrow = false; + } + + if (d.crack_model_status.frozen) // Gray's infiltration, 1985 + { + double steps_per_day = 86400.0 / global_param->dt(); + Crack crack(major, min_swe_to_freeze, infDays, + AllowPriorInf, lenstemp,steps_per_day,d.crack_model_status); + + crack.init_inputs(snowmelt, rainfall, swe, soil_storage_at_freeze, + airtemp, is_new_day()); + d.crack_model_status.daily_melt_total = snowmelt * steps_per_day; + crack.run(); + + runoff = crack.get_runoff() / steps_per_day; + melt_runoff = crack.get_melt_runoff() / steps_per_day; + inf = crack.get_inf() / steps_per_day; + snowinf = crack.get_snow_inf() / steps_per_day; + rain_on_snow = crack.get_rain_on_snow(); + + Increment_Totals(d,runoff,melt_runoff,inf,snowinf,rain_on_snow); + + if (is_new_day() && swe <= 0.0 && + d.crack_model_status.major_melt_count > 0) + { + d.crack_model_status.end_freeze(); + d.crack_model_status.end_freeze_tomorrow = true; + } + //if (swe <= 0.0 && d.crack_model_status.major_melt_count > 0) + // d.crack_model_status.end_freeze_tomorrow = true; + + //if (is_new_day()) + //{ + // d.last_day = global_param->day(); + //} + } + else if (thaw_type == AYERS) // if not frozen, do Ayers + { + Ayers ayers(rainfall, snowmelt, d.texture, d.ground_cover, SoilDataObj); + + ayers.run(); + + inf = ayers.get_inf(); + runoff = ayers.get_runoff(); + + //Below exists to match the function of CRHM + //snowmelt infiltration computed by crack continues for the next day after SWE = 0.0 + if (is_new_day()) + d.crack_model_status.end_freeze_tomorrow = false; + + if (d.crack_model_status.end_freeze_tomorrow) + { + double steps_per_day = 86400.0 / global_param->dt(); + snowinf = d.crack_model_status.current_snow_inf / + steps_per_day; + melt_runoff = d.crack_model_status.current_melt_runoff / + steps_per_day; + inf += snowinf - snowmelt; + runoff += melt_runoff; + } + else + snowinf = ayers.get_snow_inf(); + + // Increment totals + Increment_Totals(d,runoff,melt_runoff,inf,snowinf,rain_on_snow); + } + else if (thaw_type == GREENAMPT) // if not frozen, do GreenAmpt + { + d.GA_temp = std::make_unique(); + + if(rainfall > 0.0) { + d.GA_temp->intensity = convert_to_rate_hourly(rainfall); + + if(d.soil_type == "pavement"){ // TODO Not a real option, handle this + // ,this is a string handle pavement separately + runoff = rainfall; + } + else if(is_space_in_dry_soil(d.soil_storage,d.soil_storage_max,rainfall)){ + inf = rainfall; + } + else { + + //double F0 = d.soil_storage; + Initialize_GA_Variables(d); + + // TODO what about ponding + if (d.GA_temp->intensity > d.GA_temp->initial_rate) { // ponding is ongoing + + d.GA_temp->final_storage = d.GA_temp->initial_storage + rainfall; + + double ponding_time = global_param->dt() / 3600.0; + + find_final_storage(d,d.GA_temp,d.GA_temp->initial_storage,ponding_time); + + d.GA_temp->pond = rainfall - (d.GA_temp->final_storage - d.GA_temp->initial_storage); + + // TODO CRHM version calls find_final_storage again here, but appears unchanged. + // Note that during tests that this could cause a difference. + } + else { + + d.GA_temp->final_storage = d.GA_temp->initial_storage + rainfall; + d.GA_temp->final_rate = calc_GA_infiltration_rate(d,d.GA_temp,d.GA_temp->final_storage); //TODO calcf1 not a function anymore + + if (d.GA_temp->intensity > d.GA_temp->final_rate) { // ponding starts midway through the time step + initialize_ponding_vars(d,d.GA_temp); + + double ponding_time = global_param->dt() / 3600.0 - d.GA_temp->time_to_ponding; + + find_final_storage(d,d.GA_temp,d.GA_temp->storage_at_ponding,ponding_time); + + d.GA_temp->pond = rainfall - (d.GA_temp->final_storage - d.GA_temp->initial_storage); + } + + } + + + inf = d.GA_temp->final_storage - d.soil_storage; + if(d.GA_temp->pond > 0.0){ + runoff = d.GA_temp->pond; + } + } + + + // Increment totals + Increment_Totals(d,runoff,melt_runoff,inf,snowinf,rain_on_snow); + d.soil_storage += d.GA_temp->final_storage; + + melt_to_infil(inf,snowinf,snowmelt); + + d.GA_temp.reset(); + + } // if(net_rain[hh] + net_snow[hh] > 0.0) greenampt routine + } + + + + // set variables to face + (*face)["total_excess"_s]=d.total_excess; + (*face)["total_meltexcess"_s]=d.total_meltexcess; + (*face)["total_inf"_s]=d.total_inf; + (*face)["total_snowinf"_s]=d.total_snowinf; + (*face)["total_rain_on_snow"_s]=d.total_rain_on_snow; + + (*face)["runoff"_s]=runoff; + (*face)["inf"_s]=inf; + (*face)["rain_on_snow"_s]=rain_on_snow; + (*face)["snowinf"_s]=snowinf; + (*face)["melt_runoff"_s]=melt_runoff; + (*face)["frozen"_s]=static_cast(d.crack_model_status.frozen); + (*face)["major_melt_count"_s]=d.crack_model_status.major_melt_count; + + if (thaw_type == GREENAMPT) + (*face)["soil_storage"_s]=d.soil_storage; +} + +//General Functions +void Infil_All::Increment_Totals(Infil_All::data &d, double &runoff, double &melt_runoff, double &inf, double &snowinf, double &rain_on_snow) { + d.total_inf += inf; + d.total_excess += runoff; + d.total_snowinf += snowinf; + d.total_meltexcess += melt_runoff; + d.total_rain_on_snow += rain_on_snow; +} + +void Infil_All::melt_to_infil(double& inf,double& snowinf,double& snowmelt) +{ + // This function determines what to do with the excess melt after swe = 0 when the crack model shuts off + // Logan suggested setting it as runoff + if (snowmelt > 0.0) + { + inf += snowmelt; + snowinf += snowmelt; + } +}; + +// Crack Functions + + +bool Infil_All::is_new_day() +{ + // TODO This has hard coded elements, Chris suggested something different here: https://godbolt.org/z/3c51T1avT + int td = global_param->posix_time().time_of_day().total_seconds(); + int time_to_midnight = 86400 - td; + if (td >= 0 && td < global_param->dt()) //(time_to_midnight >= global_param->dt()) + { + return true; + } + else + return false; +}; + + +// Ayers + +// Green-Ampt Functions + +double Infil_All::convert_to_rate_hourly(double &rainfall) { + return rainfall / (global_param->dt() / 3600.0); +} + +bool Infil_All::is_space_in_dry_soil(double &moist, double &max, double &rainfall) { + return moist == 0.0 && max >= rainfall; +} + +void Infil_All::Initialize_GA_Variables(Infil_All::data &d) { + // This function requires d.soil_storage so the full object d must be passed. + // For simple reading, defined GA pointer to be consistent with other functions that use GA + // rather than GA_temp + // + // TODO Make this a constructor for the tempvars struct + std::unique_ptr &GA = d.GA_temp; + + GA->soil_storage_deficit = (1.0 - d.soil_storage/d.soil_storage_max); // TODO GA in Dingman is porosity - pore space filed + // Here: 1.0 means we've filled all the pores + // 0.4 - 0.2 = 0.2 (porosity) + // 1.0 - 0.5/1.0 = 0.5 (current) + // Is this a problem? + GA->initial_rate = calc_GA_infiltration_rate(d,GA,d.soil_storage); + GA->initial_storage = d.soil_storage; + GA->final_storage = GA->initial_storage; + GA->final_rate = GA->initial_rate; + GA->capillary_suction = SoilDataObj.capillary_suction(d.soil_type) + * GA->soil_storage_deficit; +} + +void Infil_All::initialize_ponding_vars(Infil_All::data& d,std::unique_ptr &GA) { + GA->storage_at_ponding = d.ksaturated * GA->capillary_suction / (GA->intensity - d.ksaturated); + GA->time_to_ponding = (GA->storage_at_ponding - GA->initial_storage)/GA->intensity; +} + +void Infil_All::find_final_storage(Infil_All::data& d,std::unique_ptr &GA, \ + double &initial_storage, double &dt) { + + double LastF1; + + do { + + LastF1 = GA->final_storage; + + GA->final_storage = initial_storage + d.ksaturated*dt + GA->capillary_suction * \ + log((GA->final_storage + GA->capillary_suction) \ + / (initial_storage + GA->capillary_suction)); + + } while(fabs(LastF1 - GA->final_storage) > 0.001); + +} + +double Infil_All::calc_GA_infiltration_rate(Infil_All::data& d,std::unique_ptr &GA, double &F){ + + return d.ksaturated*(GA->capillary_suction/F + 1.0); + +} + +//End Green-Ampt functions diff --git a/src/modules/Infil_All.hpp b/src/modules/Infil_All.hpp new file mode 100644 index 00000000..a6de5a34 --- /dev/null +++ b/src/modules/Infil_All.hpp @@ -0,0 +1,170 @@ +// +// Canadian Hydrological Model - The Canadian Hydrological Model (CHM) is a novel +// modular unstructured mesh based approach for hydrological modelling +// Copyright (C) 2018 Christopher Marsh +// +// This file is part of Canadian Hydrological Model. +// +// Canadian Hydrological Model is free software: you can redistribute it and/or +// modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// Canadian Hydrological Model is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with Canadian Hydrological Model. If not, see +// . +// + +#pragma once + +#include "logger.hpp" +#include "triangulation.hpp" +#include "module_base.hpp" +#include "TPSpline.hpp" +#include +#include "Soil.h" +#include "Crack.hpp" +#include "Ayers.hpp" +#include "boost/date_time/posix_time/posix_time_types.hpp" +/** + * \ingroup modules infil soils exp + * @{ + * \class Infil_All + * + * TODO some of this is correct, but likely will be outdated eventually. + * Estimates areal snowmelt infiltration into frozen soils for: + * a) Restricted - Water entry impeded by surface conditions + * b) Limited - Capiliary flow dominates and water flow influenced by soil physical properties + * c) Unlimited - Gravity flow dominates + * + * **Depends:** + * - Snow water equivalent "swe" [mm] + * - Snow melt for interval "snowmelt_int" [\f$mm \cdot dt^{-1}\f$] + * + * **Provides:** + * - Infiltration "inf" [\f$mm \cdot dt^{-1}\f$] + * - Total infiltration "total_inf" [mm] + * - Total infiltration excess "total_excess" [mm] + * - Total runoff "runoff" [mm] + * - Total soil storage "soil_storage" + * - Potential infiltration "potential_inf" + * - Opportunity time for infiltration to occur "opportunity_time" + * - Available storage for water of the soil "available_storage" + * + * \rst + * .. note:: + * Has hardcoded soil parameters that need to be read from the mesh parameters. + * + * \endrst + * + * **References:** + * - Gray, D., Toth, B., Zhao, L., Pomeroy, J., Granger, R. (2001). Estimating areal snowmelt infiltration into frozen soils + * Hydrological Processes 15(16), 3095-3111. https://dx.doi.org/10.1002/hyp.320 + * @} + */ +class Infil_All : public module_base +{ +REGISTER_MODULE_HPP(Infil_All) +public: + Infil_All(config_file cfg); + + ~Infil_All(); + + void run(mesh_elem &face); + void init(mesh& domain); + + class data : public face_info + { + public: + struct tempvars + { + double intensity; + double soil_storage_deficit; + double capillary_suction; + double initial_rate; + double initial_storage; + double final_rate; + double final_storage; + double pond; + double storage_at_ponding; + double time_to_ponding; + }; + double total_inf; + double total_excess; + double total_snowinf; + double total_meltexcess; + double total_rain_on_snow; + + std::string soil_type; + + // Crack + Crack::info crack_model_status; + int last_day; + bool end_freeze_tomorrow = false; // Delays end of freeze by one day so that we can get the melt distributed over the day. + + // Ayers + std::string texture; + std::string ground_cover; + + // GreenAmpt + double soil_storage; + double soil_storage_max; + double ksaturated; + std::unique_ptr GA_temp{nullptr}; + + }; + +private: + + const Soil::soils_na& SoilDataObj = Soil::get_soil_obj(); + + // General + bool is_new_day(void); + + // Crack + double major; + double min_swe_to_freeze; + unsigned int infDays; + bool AllowPriorInf; + double lenstemp; + + // General, thawed soil + enum ThawOptions { AYERS, GREENAMPT}; + unsigned int thaw_type; + + + // GreenAmpt + // double max_soil_storage; + // double soil_depth; + // double porosity; + enum GATable {PSI, KSAT, WILT, FCAP, PORG, PORE, AIENT, PORESZ, AVAIL}; // Used for mapping the soil table, PSI and KSAT are used, the others are unused but may but used in the future or other modules. + enum GAVars {TOTINF, RATEINF, SUCTION, THETA}; + //double soilproperties[][9]; + //double textureproperties[][6]; + + // TODO I should put all these functions on the data class + // General Functions + void Increment_Totals(data &d, double &runoff, double &melt_runoff, double &inf, double &snowinf, double &rain_on_snow); + void melt_to_infil(double& inf,double& snowinf,double& snowmelt); + + + // Green-Ampt Functions + double convert_to_rate_hourly(double &rainfall); + bool is_space_in_dry_soil(double &moist, double &max, double &rainfall); + void Initialize_GA_Variables(data &d); + void initialize_ponding_vars(data &d, std::unique_ptr &GA); + void find_final_storage(data &d, std::unique_ptr &GA, \ + double &initial_storage, double &dt); + double calc_GA_infiltration_rate(data &d, std::unique_ptr &GA, double &F); + + + + + +}; diff --git a/src/modules/infil_submodules/Ayers.hpp b/src/modules/infil_submodules/Ayers.hpp new file mode 100644 index 00000000..d9b3f134 --- /dev/null +++ b/src/modules/infil_submodules/Ayers.hpp @@ -0,0 +1,30 @@ +#include "submodule_base.hpp" + +template +class Ayers : public submodule_base +{ +public: + Ayers(const double& _rainfall,const double& _snowmelt, + const std::string& _texture,const std::string& _ground_cover, + const soil_data& _soil_data); + + ~Ayers() {}; + + virtual void run() override; + + const double& rainfall; + const double& snowmelt; + const std::string& texture; + const std::string& ground_cover; + const soil_data& soils; + + double runoff = 0.0; + double get_runoff() { return runoff; }; + double inf = 0.0; + double get_inf() { return inf; }; + double snow_inf = 0.0; + double get_snow_inf() { return snow_inf; }; + +}; + +#include "Ayers.ipp" diff --git a/src/modules/infil_submodules/Ayers.ipp b/src/modules/infil_submodules/Ayers.ipp new file mode 100644 index 00000000..8a40228a --- /dev/null +++ b/src/modules/infil_submodules/Ayers.ipp @@ -0,0 +1,35 @@ +#pragma once + +template +Ayers::Ayers(const double& _rainfall, + const double& _snowmelt, const std::string& _texture, + const std::string& _ground_cover, const soil_data& _soil_data) : + rainfall(_rainfall), + snowmelt(_snowmelt), + texture(_texture), + ground_cover(_ground_cover), + soils(_soil_data) {}; + +template +void Ayers::run() +{ + if (rainfall == 0.0 && snowmelt == 0.0) + return; + + if (rainfall > 0.0) + { + double maxinfil = (soils.*max_infil)(texture,ground_cover); + inf = std::min(maxinfil,rainfall); + runoff = rainfall - inf; + if (runoff < 1e-12) + runoff = 0.0; + + } + + if (snowmelt > 0.0) + { + inf += snowmelt; + snow_inf = snowmelt; + } + +} diff --git a/src/modules/infil_submodules/Crack.cpp b/src/modules/infil_submodules/Crack.cpp new file mode 100644 index 00000000..11a3ff71 --- /dev/null +++ b/src/modules/infil_submodules/Crack.cpp @@ -0,0 +1,183 @@ +#include "Crack.hpp" + +Crack::Crack(double& _major, double& _min_swe_to_freeze, unsigned int& _infDays, bool& _AllowPriorInf, double& _lenstemp,double& _steps_per_day, Crack::info& _d) : major(_major), min_swe_to_freeze(_min_swe_to_freeze), infDays(_infDays), AllowPriorInf(_AllowPriorInf), lenstemp(_lenstemp), steps_per_day(_steps_per_day), d(_d) +{ + +}; + +void Crack::init_inputs(double _snowmelt,double _rainfall, double _swe, double _soil_storage_at_freeze, double _airtemp, bool _newday) +{ + snowmelt = _snowmelt; + rainfall = _rainfall; + swe = _swe; + soil_storage_at_freeze = _soil_storage_at_freeze; + airtemp = _airtemp; + is_newday = _newday; +}; + +void Crack::run() +{ + + d.daily_rain_total += rainfall; + if (is_newday) // Gray's infiltration, 1985 + { + + + if (d.daily_melt_total > 0.0) + { + + if (soil_storage_at_freeze == 0) // Unlimited + { + inf += d.daily_melt_total; + d.major_melt_count = 1; + } + else if (soil_storage_at_freeze > 0 && soil_storage_at_freeze < 100) // Limited + { + + increment_major_count(); + + Check_for_ice_lens(); + + if (is_first_major()) + { + SPDLOG_DEBUG("First Major"); + Calc_Index(); + Calc_Actual_Inf(); + //increment_major_count(); + } + else if (is_limited_phase()) + { + SPDLOG_DEBUG("Limited Phase"); + Calc_Actual_Inf(); + + //increment_major_count(); + } + else if (is_prior_first_major() && AllowPriorInf) + { + SPDLOG_DEBUG("Prior"); + inf = d.daily_melt_total; + } + + } + else if (soil_storage_at_freeze == 100) // Restricted + { + inf = 0.; + d.major_melt_count = 1; + } + + + runoff = d.daily_melt_total - inf; + // melt_runoff and snowinf only track melt related quantities + // total runoff and infiltrated amounts from ANY source are stored in + // inf and runoff + + // This is a weird function, if there is any snowinf, then the rain on the snow also infiltrates + // this is ported directly from CRHM module crack + if (inf > 0.0) + { + inf += d.daily_rain_total; + } + else + { + runoff += d.daily_rain_total; + } + melt_runoff = runoff; + snow_inf = inf; + + + } + d.yesterday_melt = d.daily_melt_total; + d.current_inf = inf; + d.current_snow_inf = snow_inf; + d.current_runoff = runoff; + d.current_melt_runoff = melt_runoff; + + rain_on_snow = d.daily_rain_total; + d.daily_melt_total = 0.0; + d.daily_rain_total = 0.0; + + } + else + { + runoff = d.current_runoff; + melt_runoff = d.current_melt_runoff; + inf = d.current_inf; + snow_inf = d.current_snow_inf; + } + + if (!is_CRHM_compare_test) + d.daily_melt_total += snowmelt; + update_t_max(); +}; + + +void Crack::Calc_Index() { + d.index = 5 * (1 - soil_storage_at_freeze/100.0) * std::pow(swe,0.584); + // d.major_major_per_melt is obtained by dividing d.index by the + // total number of time steps to get to d.index + // This only works if 86400 / dt is a fraction which turns infDays into an integer + // Example: if dt is 345600 (4 days in seconds) and infDays is 6 days. + // the denominator is 1.5, which then requires 2 major melts for it to stop, not 1.5 + // this is actually OK behaviour, but difficult to understand. + + d.max_major_per_melt = d.index / infDays; + d.index = std::min(d.index / swe,1.0); + d.init_SWE = swe; +}; + +void Crack::Calc_Actual_Inf() { + inf = d.daily_melt_total * d.index; + if (inf > d.max_major_per_melt && is_major_melt()) { + inf = d.max_major_per_melt; + } +}; + + +void Crack::update_t_max() +{ + if (is_newday) + d.tmax = airtemp; + else + d.tmax = std::max(d.tmax,airtemp); +}; + +void Crack::Check_for_ice_lens() +{ + if (d.major_melt_count > 0 && d.tmax < lenstemp) + { + SPDLOG_DEBUG("Ice lens found"); + d.major_melt_count = infDays + 4; + } + d.tmax = 0.0; +}; + +bool Crack::is_first_major() +{ + // TODO This isn't really `is_first_major()`, only the part to the left of OR. Shoudl split this into two functions + return (is_major_melt() && swe >= d.init_SWE && + (is_prior_first_major() || is_limited_phase()));//( (d.major_melt_count == 0) & (is_major_melt()) ) || ( (swe >= d.init_SWE) & (is_limited_phase())); +}; + +bool Crack::is_major_melt() +{ + return d.daily_melt_total > major; +}; + +void Crack::increment_major_count() +{ + if (is_major_melt()) + { + d.major_melt_count++; + } +}; + + +bool Crack::is_limited_phase() +{ + return d.major_melt_count > 0 && d.major_melt_count <= infDays; +}; + +bool Crack::is_prior_first_major() +{ + return d.major_melt_count == 0; +}; diff --git a/src/modules/infil_submodules/Crack.hpp b/src/modules/infil_submodules/Crack.hpp new file mode 100644 index 00000000..383767b6 --- /dev/null +++ b/src/modules/infil_submodules/Crack.hpp @@ -0,0 +1,115 @@ + +#include "submodule_base.hpp" + +class Crack : submodule_base +{ +public: + struct info; + Crack(double& _major,double& _min_swe_to_freeze,unsigned int& _infDays,bool& _AllowPriorInf,double& _lenstemp,double& _steps_per_day,info& _d); + ~Crack() {}; + + virtual void run() override; + void init_inputs(double _snowmelt,double _rainfall, double _swe, double _soil_storage_at_freeze, double _airtemp, bool _newday); + + double runoff = 0.0; + double get_runoff() { return runoff; }; + double melt_runoff = 0.0; + double get_melt_runoff() { return melt_runoff; }; + double inf = 0.0; + double get_inf() { return inf; }; + double snow_inf = 0.0; + double get_snow_inf() { return snow_inf; }; + double rain_on_snow = 0.0; + double get_rain_on_snow() { return rain_on_snow; }; + + double snowmelt; + double rainfall; + double swe; + double soil_storage_at_freeze; + double airtemp; + bool is_newday; + bool is_CRHM_compare_test = false; + + const double& major; + const double& min_swe_to_freeze; + const unsigned int& infDays; + const bool& AllowPriorInf; + const double& lenstemp; + const double& steps_per_day; + + struct info + { + bool frozen; + unsigned int major_melt_count; + double index; + double max_major_per_melt; + double init_SWE; + double daily_melt_total; + double daily_rain_total; + bool current_day_is_major; + double tmax; + double current_inf; + double current_snow_inf; + double current_runoff; + double current_melt_runoff; + double yesterday_melt; + bool end_freeze_tomorrow; + + void init() + { + frozen = false; + major_melt_count = 0; + index = 0.0; + max_major_per_melt = 0.0; + init_SWE = 0.0; + daily_melt_total = 0.0; + daily_rain_total = 0.0; + current_day_is_major = false; + tmax = 0.0; + current_inf = 0.0; + current_snow_inf = 0.0; + current_runoff = 0.0; + current_melt_runoff = 0.0; + yesterday_melt = 0.0; + end_freeze_tomorrow = false; + }; + + void begin_freeze() + { + frozen = true; + index = 0.0; + max_major_per_melt = 0.0; + init_SWE = 0.0; + current_inf = 0.0; + current_runoff = 0.0; + yesterday_melt = 0.0; + end_freeze_tomorrow = false; + }; + + void end_freeze() + { + frozen = false; + major_melt_count = 0; + }; + + }; + info& d; + + void melt_to_infil(); + + // Crack Functions + void Calc_Index(); + void Calc_Actual_Inf(); + void Check_for_ice_lens(); + void increment_major_count(); + bool is_first_major(); + bool is_limited_phase(); + bool is_prior_first_major(); + bool is_major_melt(); + void daily_melt_increment(); + void update_t_max(); + + +}; + + diff --git a/src/modules/soil_module.cpp b/src/modules/soil_module.cpp new file mode 100644 index 00000000..ffa2af6d --- /dev/null +++ b/src/modules/soil_module.cpp @@ -0,0 +1,459 @@ +#include "soil_module.hpp" + +REGISTER_MODULE_CPP(soil_module); + +soil_module::soil_module(config_file cfg) : module_base("soil_module", parallel::data, cfg) +{ + depends("swe"); + depends("ET"); + depends("inf"); + depends("runoff"); + depends("surface_temperature"); + // depends("routing_residual"); + + provides("condensation"); + provides("actual_soil_ET"); + provides("soil_excess_to_runoff"); + provides("soil_excess_to_gw"); + provides("runoff_to_depression"); + provides("ground_water_out"); + provides("soil_to_ssr"); + provides("rechr_to_ssr"); + provides("soil_storage"); + provides("soil_rechr_storage"); + provides("depression_storage"); + provides("ground_water_storage"); + provides("detention_storage"); + provides("K_rechr_to_ssr"); + provides("K_lower_to_ssr"); + provides("K_detention_to_runoff"); + provides("K_depression_to_ssr"); + provides("K_depression_to_gw"); + provides("K_ground_water_out"); + provides("K_soil_to_gw"); + provides("thaw_front_depth"); + provides("freeze_front_depth"); + provides("first_front_depth"); + provides("thaw_fraction_rechr"); + provides("thaw_fraction_lower"); +}; + +soil_module::~soil_module() +{ + +}; + +void soil_module::init(mesh& domain) +{ + + SoilDataObj = std::make_unique(); + + for (size_t i = 0; i < domain->size_faces(); i++) + { + auto face = domain->face(i); + auto& d = face->make_module_data(ID); + // I do some evil things here to allow for the submodules to access module_base functions like is_water + // A pointer to face is put in d, likewise a pointer to this instance of this class is also added, see the overridden functions + // get_dt and is_lake below. + // Changed is_lake to a variable from a function, so it stores the result of is_water rather than requiring a copy of face. + //d.my_face = &face; + set_local_module(d); + set_soil_params(face,d); + + // dependency injection of K_estimation + d.K_estimator = std::make_unique(d); + d.soil_layers = std::make_unique(d,*d.K_estimator); + // ET coud (should) have been dependecy injection. So TODO + d.ET = std::make_unique(d); + //if (d.soil_layers) Removed if statments for now because ET and soil modules are coupled. + //if (d.ET) + set_ET_params(face,d); + + initial_soil_conditions(face,d); + + init_param_state_XG(face,d); + + set_soil_outputs(face,d); + + } +}; + + + +void soil_module::run(mesh_elem& face) +{ + + auto& d = face->make_module_data(ID); + + // TODO new newday stuff and is last day new and is crhm test + XG_algorithm XG = get_XG(face,d); + XG.run(); + get_soil_inputs(face,d,XG); + + + // order of operations here is hard coded, but it wouldn't be physically wrong to impose ET before soil + // This is fine because this soil module is run in this order. + + d.soil_layers->run(); + + if (d.swe == 0.0) + d.ET->run(); + else + d.actual_soil_ET = 0.0; + + + + set_soil_outputs(face,d); + + +}; + +bool soil_module::is_new_day() +{ + // TODO This has hard coded elements, Chris suggested something different here: https://godbolt.org/z/3c51T1avT + int td = global_param->posix_time().time_of_day().total_seconds(); + int time_to_midnight = 86400 - td; + if (td >= 0 && td < global_param->dt()) //(time_to_midnight >= global_param->dt()) + { + return true; + } + else + return false; +}; + +void soil_module::get_soil_inputs(mesh_elem& face,soil_module::data& d,XG_algorithm& XG) +{ + d.swe = (*face)["swe"_s]; + d.thaw_front_depth = XG.get_thaw_depth(); //(*face)["thaw_front_depth"_s]; + d.freeze_front_depth = XG.get_freeze_depth(); //(*face)["freeze_front_depth"_s]; + d.freeze_thaw_first_front = XG.get_first_front_depth(); //(*face)["first_front_depth"_s]; + d.actual_ET = (*face)["ET"_s]; + d.infil = (*face)["inf"_s]; + d.runoff = (*face)["runoff"_s]; + d.routing_residual = 0.0; //(*face)["routine_residual"_s]; + d.is_lake = is_water(face); +}; + +void soil_module::set_soil_outputs(mesh_elem& face,soil_module::data& d) +{ + (*face)["condensation"_s] = d.condensation; + (*face)["actual_soil_ET"_s] = d.actual_soil_ET; + (*face)["soil_excess_to_runoff"_s] = d.soil_excess_to_runoff; + (*face)["soil_excess_to_gw"_s] = d.soil_excess_to_gw; + (*face)["runoff_to_depression"_s] = d.runoff_to_depression; + (*face)["ground_water_out"_s] = d.ground_water_out; + (*face)["soil_to_ssr"_s] = d.soil_to_ssr; + (*face)["rechr_to_ssr"_s] = d.rechr_to_ssr; + (*face)["soil_storage"_s] = d.soil_storage; + (*face)["soil_rechr_storage"_s] = d.soil_rechr_storage; + (*face)["depression_storage"_s] = d.depression_storage; + (*face)["ground_water_storage"_s] = d.ground_water_storage; + (*face)["detention_storage"_s] = d.detention_storage; + (*face)["K_rechr_to_ssr"_s] = d.K_rechr_to_ssr; + (*face)["K_lower_to_ssr"_s] = d.K_lower_to_ssr; + (*face)["K_detention_to_runoff"_s] = d.K_detention_to_runoff; + (*face)["K_depression_to_ssr"_s] = d.K_depression_to_ssr; + (*face)["K_depression_to_gw"_s] = d.K_depression_to_gw; + (*face)["K_ground_water_out"_s] = d.K_ground_water_out; + (*face)["K_soil_to_gw"_s] = d.K_soil_to_gw; + + (*face)["thaw_fraction_rechr"_s] = d.thaw_fraction_rechr; + (*face)["thaw_fraction_lower"_s] = d.thaw_fraction_lower; + // XG out + (*face)["thaw_front_depth"_s] = d.thaw_front_depth; + (*face)["freeze_front_depth"_s] = d.freeze_front_depth; + (*face)["first_front_depth"_s] = d.freeze_thaw_first_front; +}; + +void soil_module::set_soil_params(mesh_elem& face, soil_module::data& d) +{ + // TODO actually connect to stuff + if (face->has_soil()) + { + d.soil_storage_max = face->soil_attribute("soil_storage_max"_s); + d.soil_rechr_max = face->soil_attribute("soil_rechr_max"_s); + d.excess_to_ssr = face->soil_attribute("excess_to_ssr"_s); + d.detention_snow_max = face->soil_attribute("detention_snow_max"_s); + d.detention_organic_max = face->soil_attribute("detention_organic_max"_s); + d.depression_max = face->soil_attribute("depression_max"_s); + d.ground_water_max = face->soil_attribute("ground_water_max"_s); + d.local_slope = face->soil_attribute("local_slope"_s)*3.14159265/180; + + d.pore_size_dist = face->soil_attribute("PSD_K_estimator"); + d.pore_size_dist_organic = face->soil_attribute("PSD_K_organic"); + const std::string soil_type = face->soil_attribute("soil_type"_s,"soils"); + d.porosity = SoilDataObj->porosity(soil_type); + d.soil_index = face->soil_attribute("soil_index"); + d.snow_grain_diameter = face->soil_attribute("snow_grain_diameter"); + + d.Ksaturated_rechr = face->soil_attribute("Ksaturated_rechr"); + d.Ksaturated_lower = face->soil_attribute("Ksaturated_lower"); + d.Ksaturated_ground_water = face->soil_attribute("Ksaturated_ground_water"); + d.Ksaturated_organic = face->soil_attribute("Ksaturated_organic"); + // Note: Ksaturated_snow is computed in K_estimate + + d.allow_runoff_from_infiltration = face->soil_attribute("allow_runoff_from_infiltration"_s); + } + else + { + d.soil_storage_max = 0.0; + d.soil_rechr_max = 0.0; + d.excess_to_ssr = 0.0; + d.detention_snow_max = 0.0; + d.detention_organic_max = 0.0; + d.depression_max = 0.0; + d.ground_water_max = 0.0; + d.local_slope = 0.0; + + d.pore_size_dist = 0.0; + d.pore_size_dist_organic = 0.0; + d.soil_index = 0.0; + d.snow_grain_diameter =0.0; + + d.Ksaturated_rechr = 0.0; + d.Ksaturated_lower = 0.0; + d.Ksaturated_ground_water = 0.0; + d.Ksaturated_organic = 0.0; + + } + +}; + + +void soil_module::set_ET_params(mesh_elem& face, soil_module::data& d) +{ + + d.ground_cover_type = face->soil_attribute("soil_groundcover_ET"_s); + std::string type = face->soil_attribute("soil_type"_s,"soils"); + + // I'm about to do something very evil. + // CRHMs inconsistent soil typing is responsible + // likely the result of empirical techniques not + // always classifying soils using the same exact categories + // soil_ET, from CHRM SoiX, only takes sand, loam, or clay + // physics/Soil.h takes 11 types, I use those types to infer the type here + + + // NOTE: This soil_type is not the same as the soil_type used for other parameters like porosity. This is purely used for a computation in the ET calculation in the soil. + int soil_type = 100; + bool sandy = compare_substring(type,"sand"); + bool loamy = compare_substring(type,"loam"); + bool clayy = compare_substring(type,"clay"); + + + if (sandy && !loamy && !clayy) // mainly sand + { + soil_type = 1; + } + else if (sandy && loamy) // loamy sand or sandy loam + { + if (type.substr(0,4) == "loam") //loamy sand + { + soil_type = 1; + } + else //sandy loam + soil_type = 2; + } + else if (sandy && clayy) // sandy clay + soil_type = 3; + else if (clayy && loamy) // loam with clay (clay-y loam) + soil_type = 2; + else if (clayy) // mainly clay + soil_type = 3; + else if (loamy) // mainly loam + soil_type = 2; + + // detaul is 100, this is handled by soil_ET by treating the soil as mainly organic + d.soil_type_rechr = soil_type; + d.soil_type_lower = soil_type; //Same for now + +}; + +int soil_module::compare_substring(std::string& type, std::string sub) +{ + return type.find(sub) != std::string::npos; // find returns npos +}; + + +void soil_module::initial_soil_conditions(mesh_elem& face, soil_module::data& d) +{ + // TODO actually connect to stuff + // requires MESHER or at least data for one station + if (face->has_soil()) + { + d.soil_storage = face->soil_attribute("soil_storage"_s); + d.soil_rechr_storage = face->soil_attribute("soil_rechr_storage"_s); + //d.thaw_fraction_rechr = face->soil_attribute("thaw_fraction_rechr"_s); + //d.thaw_fraction_lower = face->soil_attribute("thaw_fraction_lower"_s); + d.detention_snow_init = face->soil_attribute("detention_snow_init"_s); + d.detention_organic_init = face->soil_attribute("detention_organic_init"_s); + d.depression_storage = face->soil_attribute("depression_storage"_s); + d.ground_water_storage = face->soil_attribute("ground_water_storage"_s); + + } + else + { + d.soil_storage = 0.0; + d.soil_rechr_storage = 0.0; + d.thaw_fraction_rechr = 0.0; + d.thaw_fraction_lower = 0.0; + d.detention_snow_init = 0.0; + d.detention_organic_init = 0.0; + d.depression_storage = 0.0; + d.ground_water_storage = 0.0; + d.pore_size_dist = 1.0; // Set to 1 be default becuase there is a divide by zero with this value. It probably will actually be skipped so not a real worry. + d.local_slope = 0.0; + } + +}; + +//bool soil_module::data::is_lake(soil_ET_DTO& DTO) +//{ +// try +// { +// // TODO resolve this bug +// soil_module::data& d = dynamic_cast(DTO); +// //bool temp = d.local_module->is_water(*d.my_face); +// return false;//d.local_module->is_water(*d.my_face); +// } catch (const std::bad_cast& e) { +// SPDLOG_DEBUG("bad cast"); +// return false; +// } +//}; + +int soil_module::data::get_dt() +{ + if (this->local_module) + return this->local_module->global_param->dt(); + + CHM_THROW_EXCEPTION(module_error,"local_module not set in soil_module::data"); + +}; + +bool soil_module::data::get_new_day() +{ + if (this->first_day) + { + this->first_day = false; + return true; + } + else if (this->local_module) + return this->local_module->is_new_day(); + + CHM_THROW_EXCEPTION(module_error,"local_module not set in soil_module::data"); +}; + +void soil_module::set_local_module(soil_module::data& d) +{ + d.local_module = this; + //d.local_module = new soil_module::data::my_module(*this); +}; + // TODO this is just a copy of is_water and this is a bad practice but currently the is_water function is not accessible by the data class. Fix: create a separate object taht module_base inherits that contains these functions. face_info will also inherit these functions. + // soil_module::data& d = static_cast(DTO); + // bool is = false; + + // if(d.my_face->has_parameter("landcover"_s)) + // { + // int LC = face->parameter("landcover"_s); + // is = global_param->parameters.get("landcover." + std::to_string(LC) + ".is_water",false); + // } + // return is + +//}; +// +void soil_module::init_param_state_XG(mesh_elem& face, soil_module::data& d) +{ + if (face->has_soil()) + { + C.Trigthrhld = cfg.get("Trigthrhld",100.0); + + double depths = cfg.get("depth_per_layer",0.5); + double theta_default = cfg.get("theta_default",0.5); + bool uniform_conductivities = cfg.get("uniform_conductivities",true); + double dry_soil_k, sat_soil_frozen_k, sat_soil_thaw_k; + if (uniform_conductivities) + { + dry_soil_k = cfg.get("dry_soil_k"); + sat_soil_frozen_k = cfg.get("sat_soil_frozen_k"); + sat_soil_thaw_k = cfg.get("sat_soil_thaw_k"); + } + else + { + dry_soil_k = face->soil_attribute("dry_soil_k"_s); + sat_soil_frozen_k = face->soil_attribute("sat_soil_frozen_k"_s); + sat_soil_thaw_k = face->soil_attribute("sat_soil_thaw_k"_s); + } + + + C.num_layers = cfg.get("number_XG_layers",10); + std::vector dry_soil_k_vec(C.num_layers,dry_soil_k); + std::vector sat_soil_frozen_k_vec(C.num_layers,sat_soil_frozen_k); + std::vector sat_soil_thaw_k_vec(C.num_layers,sat_soil_thaw_k); + std::vector depth_vec(C.num_layers,depths); + std::vector por(C.num_layers,d.porosity); + std::vector theta_default_vec(C.num_layers,theta_default); + C.theta_min = cfg.get("moisture_content_min_per_layer",0.001); + double perma_frost_depth = face->soil_attribute("perma_frost_depth"_s); + C.freeze_kw_ki_update = cfg.get("update_k_behind_front_freeze",true); + C.thaw_ki_kw_update = cfg.get("update_k_behind_front_thaw",true); + C.k_update = cfg.get("k_update",1); + C.time_step_per_day = 86400.0/global_param->dt(); + C.calc_conductivity = cfg.get("Johansen_conductivity",false); + + sat_soil_frozen_k_vec.at(0) = 1.55; + sat_soil_thaw_k_vec.at(0) = 0.8; + + d.P = std::make_unique(depth_vec, + C.Trigthrhld, + por, + C.num_layers, + theta_default_vec, + C.theta_min, + dry_soil_k_vec, + sat_soil_frozen_k_vec, + sat_soil_thaw_k_vec, + 0.0, //SWE_k Not actually used but still part of the module for now + perma_frost_depth, + C.freeze_kw_ki_update, + C.thaw_ki_kw_update, + C.k_update, + d.soil_rechr_max, + d.soil_storage_max, + C.time_step_per_day, + C.calc_conductivity); + + d.P->is_crhm_test = true; + d.S = std::make_unique(d.P->N_Soil_layers,*(d.P)); + + d.S->set_layer_moisture_maximums(*(d.P)) + .set_thermal_conductivities(*(d.P),d.soil_storage,d.soil_rechr_storage) + .set_freezethaw_ratios(*(d.P)); + + double Zdf_init = face->soil_attribute("init_freeze_front_depth"_s); + double Zdt_init = face->soil_attribute("init_thaw_front_depth"_s); + + + XG_algorithm XG(0.0,0.0,0.0,*(d.S),*(d.P)); + + + XG.init_freezethaw_degreedays(Zdf_init,Zdt_init,d.P->Zpf_init); + + } +}; + +XG_algorithm soil_module::get_XG(mesh_elem& face,soil_module::data& d) +{ + if (d.first_day) + { + d.S->is_newday = true; + } + else + d.S->is_newday = is_new_day(); + + + XG_algorithm XG((*face)["surface_temperature"_s],d.soil_storage,d.soil_rechr_storage,*(d.S),*(d.P)); + + return XG; +}; + + diff --git a/src/modules/soil_module.hpp b/src/modules/soil_module.hpp new file mode 100644 index 00000000..1045c7f9 --- /dev/null +++ b/src/modules/soil_module.hpp @@ -0,0 +1,146 @@ +// +// Canadian Hydrological Model - The Canadian Hydrological Model (CHM) is a novel +// modular unstructured mesh based approach for hydrological modelling +// Copyright (C) 2018 Christopher Marsh +// +// This file is part of Canadian Hydrological Model. +// +// Canadian Hydrological Model is free software: you can redistribute it and/or +// modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// Canadian Hydrological Model is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with Canadian Hydrological Model. If not, see +// . +// + +#pragma once + +#include "logger.hpp" +#include "triangulation.hpp" +#include "module_base.hpp" +#include "TPSpline.hpp" +#include +#include +#include "Soil.h" +#include "soil_two_layer.hpp" +#include "soil_ET.hpp" +#include "K_estimate.hpp" +#include "XG_algorithm.hpp" + +/** + * \ingroup modules infil soil_module exp + * @{ + * \class soil_module * + * + * Organizes soil process models (separate classes) + * - estimating freeze/thaw fronts using the XG algorithm: XG-freeze_thaw + * - layer per-time-step outflow: k_estimate + * - two layer soil model: soil_two_layer + * - soil Evapotranspiration: soil_ET + * + * Each are found in soil_submodules/ + * + * **Depends:** + * - Snow water equivalent "swe" [mm] + * - thaw front in soil "thaw_front_depth" [mm] + * - freeze front in soil "freeze_front_depth" [mm] + * - potential evapotranspiration "potential_ET" [mm] + * - infiltrated moisture "inf" [mm] + * - runoff/excess from infiltration scheme "runoff" [mm] + * - leftover from routing between triangles "routing_residual" [mm] + * + * **Provides:** TODO match provides with actual .cpp file + * - Infiltration "inf" [\f$mm \cdot dt^{-1}\f$] + * - Total infiltration "total_inf" [mm] + * - Total infiltration excess "total_excess" [mm] + * - Total runoff "runoff" [mm] + * - Total soil storage "soil_storage" + * - Potential infiltration "potential_inf" + * - Opportunity time for infiltration to occur "opportunity_time" + * - Available storage for water of the soil "available_storage" + * + * \rst + * .. note:: TODO add any notes + * Has hardcoded soil parameters that need to be read from the mesh parameters. + * + * \endrst + * + * **References:** TODO add any references + * - Gray, D., Toth, B., Zhao, L., Pomeroy, J., Granger, R. (2001). Estimating areal snowmelt infiltration into frozen soils + * Hydrological Processes 15(16), 3095-3111. https://dx.doi.org/10.1002/hyp.320 + * @} + */ + +class soil_module : public module_base +{ +REGISTER_MODULE_HPP(soil_module); +public: + soil_module(config_file cfg); + + ~soil_module(); + + void run(mesh_elem& face); + void init(mesh& domain); + + class data : public face_info, public two_layer_DTO, public soil_ET_DTO + { + public: + std::unique_ptr soil_layers; + std::unique_ptr ET; + // custom deletor that does nothing to make sure it doesn't try to delete the mesh_elem it points to + // REMOVED same reason as is_lake below, may 2025 + //mesh_elem* my_face;//(nullptr, [](mesh_elem* ptr) {}); + std::unique_ptr K_estimator; + // overridden + // REMOVED is_lake as a function and is now a variable storing the result of is_water (May 2025) + //bool is_lake(soil_ET_DTO& DTO) override; + int get_dt() override; + bool get_new_day() override; + // custom deletor that does nothing to make sure it doesn't try to delete the mesh_elem it points to + soil_module* local_module;//(nullptr, [](soil_module*) {}); + + std::unique_ptr P; + std::unique_ptr S; + + bool first_day = true; + }; + + void set_local_module(soil_module::data& d); + + +private: + + std::unique_ptr SoilDataObj; + + void get_soil_inputs(mesh_elem& face, data& d,XG_algorithm& XG); + void set_soil_outputs(mesh_elem& face, data& d); + void set_soil_params(mesh_elem& face, soil_module::data& d); + void set_ET_params(mesh_elem& face, soil_module::data& d); + void initial_soil_conditions(mesh_elem& face, soil_module::data& d); + + int compare_substring(std::string& type, std::string sub); + + XG_algorithm get_XG(mesh_elem& face,data& d); + void init_param_state_XG(mesh_elem& face,data& d); + bool is_new_day(); + struct XG_shared_const + { + double Trigthrhld; + size_t num_layers; + double theta_min; + size_t freeze_kw_ki_update; + size_t thaw_ki_kw_update; + size_t k_update; + size_t time_step_per_day; + bool calc_conductivity; + }; + XG_shared_const C; +}; diff --git a/src/modules/soil_submodules/I_K_estimate.hpp b/src/modules/soil_submodules/I_K_estimate.hpp new file mode 100644 index 00000000..8469e750 --- /dev/null +++ b/src/modules/soil_submodules/I_K_estimate.hpp @@ -0,0 +1,9 @@ +#pragma once +#include "submodule_base.hpp" + +class I_K_estimate : public submodule_base +{ +public: + virtual ~I_K_estimate() = default; + virtual void run(void) = 0; +}; diff --git a/src/modules/soil_submodules/I_freeze_thaw_depths.hpp b/src/modules/soil_submodules/I_freeze_thaw_depths.hpp new file mode 100644 index 00000000..ec2e7567 --- /dev/null +++ b/src/modules/soil_submodules/I_freeze_thaw_depths.hpp @@ -0,0 +1,22 @@ +#pragma once +#include "submodule_base.hpp" + +class I_freeze_thaw_depths : public submodule_base +{ +public: + virtual ~I_freeze_thaw_depths() = default; + virtual void run(void) = 0; + + + double get_freeze_depth() { return freeze_front_depth;}; + + double get_thaw_depth() { return thaw_front_depth; }; + + double get_first_front_depth() { return first_front_depth; }; + +protected: + double thaw_front_depth = 0.0; + double freeze_front_depth = 0.0; + double first_front_depth = 0.0; + +}; diff --git a/src/modules/soil_submodules/K_estimate.cpp b/src/modules/soil_submodules/K_estimate.cpp new file mode 100644 index 00000000..753624bb --- /dev/null +++ b/src/modules/soil_submodules/K_estimate.cpp @@ -0,0 +1,156 @@ +#include "K_estimate.hpp" + +K_estimate::K_estimate(two_layer_DTO& _DTO) : DTO(_DTO) +{ + Vels = std::make_unique(DTO); +}; + +void K_estimate::run(void) +{ + Vels->init_vels(); + + check_soil_zeros(); + + Vels->calculate(); + + set_K_values(*Vels); + + +}; + +void K_estimate::set_K_values(I_Darcy_Vels& Vels) +{ + double unit_changer = DTO.get_dt() * 1000.0; // m/s * s/step * mm/m = mm/step | m/s -> units of Vels.lateral_rechr (and others) + double unit_changer_lateral = unit_changer / 1000.0; // dW/A from Fang et al, (2013), only applies to lateral flow + DTO.K_rechr_to_ssr = Vels.lateral_rechr * DTO.soil_rechr_max * unit_changer_lateral; + DTO.K_lower_to_ssr = Vels.lateral_lower * (DTO.soil_storage_max - DTO.soil_rechr_max) * unit_changer_lateral; + DTO.K_depression_to_ssr = Vels.lateral_lower * DTO.soil_storage_max * unit_changer_lateral; + DTO.K_depression_to_gw = Vels.vertical_depression * unit_changer; + + DTO.K_soil_to_gw = Vels.vertical_lower * unit_changer; + + DTO.K_ground_water_out = Vels.lateral_ground_water * DTO.ground_water_storage * unit_changer_lateral; // CRHM uses ground_water_storage, rather than ground_water_max in K_estimate. This was a suggestioning by a referee because ground water is saturated flow and so the "max" doesn't really make sense in this context. Remember that this model is conceptual and not a "real" soil model. + DTO.K_detention_to_runoff = Vels.lateral_detention * DTO.detention_max * unit_changer_lateral; +}; + +void K_estimate::check_soil_zeros() +{ + // floating point errors can lead to some things being nearly zero but not quite + // this function simple ensures small values are set to zero. + + if ( DTO.soil_rechr_storage <= 0.0000001 ) + DTO.soil_rechr_storage = 0.0; + + if ( DTO.soil_storage <= 0.0000001 ) + DTO.soil_storage = 0.0; + + if ( DTO.ground_water_storage <= 0.0000001 ) + DTO.ground_water_storage = 0.0; + + if ( DTO.soil_rechr_storage > DTO.soil_storage ) + DTO.soil_rechr_storage = DTO.soil_storage; + +}; + +void I_Darcy_Vels::init_vels(void) +{ + lateral_rechr = 0.0; + + lateral_lower = 0.0; + + vertical_depression = 0.0; + + vertical_lower = 0.0; + + lateral_ground_water = 0.0; + + lateral_detention = 0.0; + +}; + +void Darcy_Vels::calculate() +{ + + if( DTO.soil_storage_max > 0.0 ) + { + if (DTO.swe > 0.0) + { + set_snow(); + } + else + set_clear(); + + } + +}; + +void Darcy_Vels::set_snow() +{ + lateral_rechr = 0.0; + + lateral_lower = get_lateral_lower(); + + vertical_depression = 0.0; + + vertical_lower = get_reused(); + + lateral_ground_water = get_lateral_ground_water(); + + if (DTO.detention_max > 0.0) + { + if (DTO.snow_density > 100) // when snowcover, use Shimizu (1970) to estimate sat. hydraulic conductivity of snow + { + lateral_detention = get_detention_snow(); + } + else + { + lateral_detention = get_detention_organic(); + } + } +}; + +void Darcy_Vels::set_clear() +{ + lateral_rechr = DTO.Ksaturated_rechr * std::pow( DTO.soil_rechr_storage/DTO.soil_rechr_max, exponent) * std::tan(DTO.local_slope); + lateral_lower = get_lateral_lower(); + + vertical_depression = get_reused(); + + vertical_lower = get_reused(); + + lateral_ground_water = get_lateral_ground_water(); + + if (DTO.detention_max > 0.0) + { + lateral_detention = get_detention_organic(); + } +}; + +double Darcy_Vels::get_lateral_lower() +{ + return DTO.Ksaturated_lower * std::pow( (DTO.soil_storage - DTO.soil_rechr_storage) / (DTO.soil_storage_max - DTO.soil_rechr_max), exponent) *std::tan(DTO.local_slope); +}; + +double Darcy_Vels::get_reused() +{ + // This Darcy Vel is used many times, so its called repeated. + return DTO.Ksaturated_lower * std::pow( DTO.soil_storage / DTO.soil_storage_max, exponent); +}; + +double Darcy_Vels::get_lateral_ground_water(void) +{ + return DTO.Ksaturated_ground_water * std::tan(DTO.local_slope); +}; + +double Darcy_Vels::get_detention_snow(void) +{ + double Ksaturated_snow = (0.077*std::pow((DTO.snow_grain_diameter/1000),2.0)*std::exp(-7.8*(DTO.snow_density/1000)))*factor; + + + return Ksaturated_snow * std::pow(DTO.detention_storage/DTO.detention_max,DTO.soil_index) * std::sin(DTO.local_slope); +} + +double Darcy_Vels::get_detention_organic(void) +{ + return DTO.Ksaturated_organic * std::pow(DTO.detention_storage/DTO.detention_max,exponent_organic) * std::tan(DTO.local_slope); +}; diff --git a/src/modules/soil_submodules/K_estimate.hpp b/src/modules/soil_submodules/K_estimate.hpp new file mode 100644 index 00000000..ee7b8bc4 --- /dev/null +++ b/src/modules/soil_submodules/K_estimate.hpp @@ -0,0 +1,73 @@ +#pragma once + +#include "I_K_estimate.hpp" +#include "soil_DTO.hpp" +#include +#include +// using an interface is pointless for now, since both are exposed, but this is more flexible in the future +// and would allow for dependecy injection from the parent module. +class I_Darcy_Vels +{ +public: + virtual ~I_Darcy_Vels() = default; + double lateral_rechr = 0.0; + + double lateral_lower = 0.0; + + double vertical_depression = 0.0; + + double vertical_lower = 0.0; + + double lateral_ground_water = 0.0; + + double lateral_detention = 0.0; + + virtual void calculate() = 0; + + I_Darcy_Vels(two_layer_DTO& _DTO) : DTO(_DTO) {}; + + void init_vels(); +protected: + two_layer_DTO& DTO; +}; + +class Darcy_Vels : public I_Darcy_Vels +{ +public: + explicit Darcy_Vels(two_layer_DTO& _DTO) : I_Darcy_Vels(_DTO), exponent(3.0 + 2.0/DTO.pore_size_dist), exponent_organic(3.0 + 2.0/DTO.pore_size_dist_organic) {}; + ~Darcy_Vels() override {}; + virtual void calculate() override; + +private: + // water density * gravity acceleration / dynamic viscosity of water + // only used exactly in this manner + const double factor = 1000.0 * 9.8 * 0.001787; + const double exponent; + const double exponent_organic; + void set_snow(); + void set_clear(); + double get_lateral_lower(); + double get_reused(); // this is an expression reused many times, stored in a function + double get_lateral_ground_water(); + double get_detention_snow(); + double get_detention_organic(); + +}; + +class K_estimate : public I_K_estimate +{ +public: + explicit K_estimate(two_layer_DTO& _DTO); + ~K_estimate() override {}; + + void run(void) override; + +private: + + two_layer_DTO& DTO; + void check_soil_zeros(); + void set_K_values(I_Darcy_Vels& Vels); + std::unique_ptr Vels; + + +}; diff --git a/src/modules/soil_submodules/XG_algorithm.cpp b/src/modules/soil_submodules/XG_algorithm.cpp new file mode 100644 index 00000000..ee85deaf --- /dev/null +++ b/src/modules/soil_submodules/XG_algorithm.cpp @@ -0,0 +1,802 @@ +#include "XG_algorithm.hpp" +XG_algorithm::XG_algorithm(const double& t, const double& _moist, + const double& _rechr, state& _S, const params& _P) : + surface_temp(t), soil_moist(_moist), soil_rechr(_rechr), + S(_S), P(_P), depths(this->P.getdepths()) +{ + depths.back() = 100.0; +}; + + + +void XG_algorithm::run() +{ + if (S.last_step_new_day) + { + S.reset_degree_day_counter(); + } + + + S.accumulate_degree_days(surface_temp); + + + if (S.is_newday) + { + S.determine_freeze_thaw_idle(); + + S.set_thermal_conductivities(P,soil_moist,soil_rechr) + .set_freezethaw_ratios(P); + + if(S.freezing()) // handle freezing + { + if(S.net_negative_degree_days()) + { + freeze(); // XG-Algorithm - Freezing + + // check for thaw front + if(S.thaw_front_overtaken()) + { + if(S.exist_excess_fronts()) + { + S.merge_freeze_fronts(this); + } + else + { + S.Zdt = 0.0; + S.Bth = 0.0; + S.Zd_front[1] = 0.0; + } + + } + } + S.Zd_front[0] = -S.Zdf; + } // freezing handled + else if(S.thawing()) // Surface thawing lower ground frozen + { + S.Bth += S.B; // Accumulate thawing degree-days + + if(S.Bth <= 0.0) + { + S.Bfr = S.Bth; + S.Zdt = 0.0; + S.Bth = 0.0; + } + else + { + thaw(); // XG-Algorithm - Thawing + + // check for freeze front + if(S.freeze_front_overtaken()) + { + if(S.exist_excess_fronts()) + S.merge_thaw_fronts(this); + else + { + S.Zdf = 0.0; + S.Bfr = 0.0; + S.Zdt = 0.0; + S.Bth = 0.0; + S.TrigState = 0; + S.Zd_front[1] = 0.0; + } + + + } + + } // if check freeze front + + S.Zd_front[0] = S.Zdt; + } // thawing handled + + S.last_step_new_day = true; + } + thaw_front_depth = S.Zdt; + freeze_front_depth = S.Zdf; + first_front_depth = S.Zd_front[0]; +}; + +void XG_algorithm::freeze(void) +{ + size_t layer = 1; + double Za; + double L = 335000; + + S.Zdf = 0.0; + + double ftc; + if (P.k_update == 2) + ftc = Interpolated_ftc(S.Zdf, layer); + else + ftc = S.ftc[layer-1]; + + Za = stefan_equation(S.Bfr,ftc,layer); + + while (Za > depths[layer-1] && layer < P.N_Soil_layers) + { + S.Zdf += depths[layer-1]; + + Za = (Za - depths[layer-1])/S.pf[layer]; + ++layer; + + if(P.k_update > 0 && P.freeze_kw_ki_update && layer > S.Fz_low) + { + S.ftc[S.Fz_low-1] = S.get_ttc(S.Fz_low-1,P); + S.ftc_contents[S.Fz_low-1] = 1; + S.pf[S.Fz_low] = std::sqrt(S.ftc[S.Fz_low-1]*S.layer_h2o[S.Fz_low]/(S.ftc[S.Fz_low]*S.layer_h2o[S.Fz_low-1])); + S.Fz_low = layer; + S.tc_composite[layer-1] = S.ftc[layer-1]; + } + } + S.Zdf += Za; + + S.Zdf = std::min(S.Zdf,P.Zpf_init); + +}; + +void XG_algorithm::thaw(void) +{ + size_t layer = 1; + double Za; + double L = 335000; + + S.Zdt = 0.0; + + double ttc; + if (P.k_update == 2) + ttc = Interpolated_ftc(S.Zdf, layer); + else + ttc = S.ttc[layer-1]; + + Za = stefan_equation(S.Bth,ttc,layer); + + while (Za > depths[layer-1] && layer < P.N_Soil_layers) + { + S.Zdt += depths[layer-1]; + + Za = (Za - depths[layer-1])/S.pt[layer]; + ++layer; + + if(P.k_update > 0 && P.thaw_ki_kw_update && layer > S.Th_low) + { + S.ttc[S.Th_low-1] = S.get_ftc(S.Th_low-1,P); + S.ttc_contents[S.Th_low-1] = 1; + S.pt[S.Th_low] = std::sqrt(S.ttc[S.Th_low-1]*S.layer_h2o[S.Th_low]/(S.ttc[S.Th_low]*S.layer_h2o[S.Th_low-1])); + S.Th_low = layer; + S.tc_composite[layer-1] = S.ttc[layer-1]; + } + } + S.Zdt += Za; + + S.Zdt = std::min(S.Zdt,P.Zpf_init); + +}; + +double XG_algorithm::stefan_equation(double& surface_index, double& thermal_conductivity,size_t& layer) +{ + const double L = 335000; + // TODO Don't know where the 211 and 8640011 comes from. + // Examining the original equation from (Xie and Gough, 2013) + // A factor of 86400 is needed + assert(S.layer_h2o[layer-1] != 0 && "About to trip divide by zero in stefan equation!"); + return std::sqrt(2.0*86400.0*thermal_conductivity * surface_index / + (S.layer_h2o[layer-1]*L)); + +}; + +double XG_algorithm::Interpolated_ttc(double Za, size_t layer) +{ + if(!P.thaw_ki_kw_update) + return (S.ttc[layer-1]); + + double split = (Za - S.Zdt)/depths[layer-1]; + if(split >= 1.0) + split = 1.0; + + double combination = S.ttc[layer-1] - split*(S.ttc[layer-1] - S.ftc[layer-1]); // thawed(18k) to frozen (4k) + + S.tc_composite2[layer-1] = combination; + + return (combination); +}; + +double XG_algorithm::Interpolated_ftc(double Za, size_t layer) { // + + if(!P.freeze_kw_ki_update) + return (S.ftc[layer-1]); + + double split = (Za - S.Zdf)/depths[layer-1]; + if(split >= 1.0) + split = 1.0; + + double combination = S.ftc[layer-1] + split*(S.ttc[layer-1] - S.ftc[layer-1]); // frozen (4k) to thawed(18k) + + S.tc_composite2[layer-1] = combination; + + return (combination); +}; + +#define tolerance 0.000001 +void XG_algorithm::find_thaw_D(double dt) { // XG-Algorithm - Thawing - used by init +// solve for Bth from Zdt using Bisection method + + if(dt == 0) + return; + + auto solution = [&dt](double Zdt) + { return Zdt - dt; }; + + double low = 0.0; + double high = 50000.0; + + do { + double mid = (high + low)/2; + S.Bth = mid; + thaw(); + if (solution(S.Zdt) > 0) + high = mid; + else + low = mid; + + } while (high - low > tolerance); + + //TODO Throw CHM exception here, indicates that Zdt is too large +}; + +void XG_algorithm::find_freeze_D(double df) { // XG-Algorithm - Thawing - used by init +// solve for Bfr from Zdt using Bisection method + if(df == 0) + return; + + auto solution = [&df](double Zdf) + { return Zdf - df; }; + + double low = 0.0; + double high = 50000.0; + + do { + double mid = (high + low)/2; + S.Bfr = mid; + freeze(); + if (solution(S.Zdf) > 0) + high = mid; + else + low = mid; + + } while (high - low > tolerance); + //TODO Throw CHM exception here, indicates Zdf is too large +}; + +void XG_algorithm::init_freezethaw_degreedays(const double& Zdf_init, const double& Zdt_init, const double& Zpf_init) +{ + if (Zdf_init > 0.0 && S.Bfr == 0.0) + find_freeze_D(Zdf_init); + + if (Zdt_init > 0.0 && Zdt_init < Zpf_init && S.Bth == 0.0) + find_thaw_D(Zdt_init); +}; + +void XG_algorithm::state::push_front(double D) { + + // if(nfront >= front_size-3){ // space to allocate plus Zdf/Zdt(2 slots) plus top of stack indicator + // string S = string("'") + Name + " (XG)' too many fronts in hru = " + to_string(hh+1).c_str(); + // CRHMException TExcept(S.c_str(), TExcept::TERMINATE); + // LogError(TExcept); + // throw TExcept; + // } + // TODO, above if statement is needed for some edge cases, but I don't know what front_size even is. I have this on my XG note TODO + + + for(size_t ii = nfront+1; 2 <= ii ; --ii) // move contents up + Zd_front[ii+1] = Zd_front[ii]; + + ++nfront; + Zd_front[2] = D; // add new entry +}; + +double XG_algorithm::state::pop_front(void) { + + double D = std::fabs(Zd_front[2]); // always positive + + for(size_t ii = 2; ii < nfront+1; ++ii) // move contents down + Zd_front[ii] = Zd_front[ii+1]; + + Zd_front[nfront+1] = 0.0; // clear memory + + --nfront; + + return D; +}; + +double XG_algorithm::state::last_front(void){ + + if(!nfront) + return 0.0; + else + return (Zd_front[2]); +}; + + +double XG_algorithm::state::get_ftc(size_t layer,const params& P){ // unfrozen(thawed) soil to be frozen + if(P.calc_conductivity){ + return (P.soil_solid_km_kw[layer] - P.soil_solid_km[layer])*std::pow(this->layer_h2o[layer]/(1000.0*P.por[layer]),2) + P.soil_solid_km[layer]; + } + else + return (1.0 - P.por[layer])*P.soil_solid_km[layer] + this->layer_h2o[layer]/1000.0*kw + (P.por[layer] - this->layer_h2o[layer]/1000.0)*ka; +}; + +double XG_algorithm::state::get_ttc(size_t layer,const params& P){ // frozen soil to be unfrozen(thawed) + if(P.calc_conductivity){ + return P.soil_solid_km[layer]*std::pow(P.soil_solid_km_ki[layer]/P.soil_solid_km[layer], this->layer_h2o[layer]/(1000.0*P.por[layer])); + } + else + return (1.0 - P.por[layer])*P.soil_solid_km[layer] + this->layer_h2o[layer]/1000.0*ki + (P.por[layer] - this->layer_h2o[layer]/1000.0)*ka; +}; + +void XG_algorithm::state::determine_freeze_thaw_idle() +{ + if (TrigAcc > P.Trigthrhld) + TrigAcc = P.Trigthrhld; + else if (TrigAcc < -P.Trigthrhld) + TrigAcc = -P.Trigthrhld; + + // Start Thaw starting from Idle + if (TrigAcc >= P.Trigthrhld / 2.0 && TrigState == 0 && (Zdf > 0.0 || nfront > 0)) + { + TrigState = 1; + Zd_front[1] = -Zdf; + t_trend = 0.0; + } + + // Start Freeze starting from Idle + if (TrigAcc <= -P.Trigthrhld / 2.0 && TrigState == 0) + { + TrigState = -1; + Zd_front[1] = Zdt; + t_trend = 0.0; + } + + // Start Idle starting from Freeze + if (TrigState == -1 && TrigAcc >= P.Trigthrhld/2.0 && t_trend > 0.0) + { + TrigState = 0; + + if (Zdt > 0.0 && Zdf > 0.0) + { + if (Zdt > Zdf) + { + push_front(Zdt); + Zdt = 0.0; + Bth = 0.0; + Zd_front[0] = 0.0; + Zd_front[1] = -Zdf; + } + } + } + + // Start Idle starting from Thaw + if (TrigState == 1 && TrigAcc <= -P.Trigthrhld / 2.0 && t_trend < 0.0) + { + TrigState = 0; + + if (Zdf > 0.0 && Zdt > 0.0) + { + if (Zdf > Zdt) + { + push_front(-Zdf); + Zdf = 0.0; + Bfr = 0.0; + Zd_front[0] = 0.0; + Zd_front[1] = Zdt; + } + } + } +}; + +bool XG_algorithm::state::freezing() +{ + return TrigState < 0; +}; + +bool XG_algorithm::state::thawing() +{ + return TrigState > 0; +}; + +bool XG_algorithm::state::net_negative_degree_days() +{ + Bfr -= B; + return Bfr > 0.0; +}; + +void XG_algorithm::state::reset_degree_day_counter() +{ + last_step_new_day = false; + B = 0.0; +}; + +void XG_algorithm::state::merge_freeze_fronts(XG_algorithm* XG) +{ + double Last = last_front(); + + if(Last < 0.0) // frozen front + { + Zdf = pop_front(); + XG->find_freeze_D(Zdf); + double Last = last_front(); + + if(Last > 0.0) // thaw front + { + Zdt = pop_front(); + XG->find_thaw_D(Zdt); + Zd_front[1] = Zdt; + } + else if(Last < 0.0) // never two frozen fronts together + { + //CRHM has thos throw an error, never two frozen fronts + } + else // no thaw front + { + Zdt = 0.0; + Bth = 0.0; + Zd_front[1] = 0.0; + } + } + else if(Last < 0.0) // never two freeze fronts together + { + // CRHM had another throw here, but I dont think its possible to reach this ever. + } + else // no thaw layer + { + Zdt = 0.0; + Bth = 0.0; + Zd_front[1] = 0.0; + } + +}; + +void XG_algorithm::state::merge_thaw_fronts(XG_algorithm* XG) +{ + double Last = last_front(); + + if(Last > 0.0) // thaw front + { + Zdt = pop_front(); + XG->find_thaw_D(Zdt); + double Last = last_front(); + + if(Last < 0.0) // frozen front + { + Zdf = pop_front(); + XG->find_freeze_D(Zdf); + Zd_front[1] = -Zdf; + } + else if(Last > 0.0) // never two thaw fronts together + { + //CRHM has thos throw an error, never two frozen fronts + } + else // no frozen front + { + Zdf = 0.0; + Bfr = 0.0; + Zd_front[1] = 0.0; + } + } + else if(Last < 0.0) // never two freeze fronts together + { + // CRHM had another throw here, but I dont think its possible to reach this ever. + } + else // no thaw layer + { + Zdf = 0.0; + Bfr = 0.0; + Zd_front[1] = 0.0; + } +}; + +bool XG_algorithm::state::exist_excess_fronts() +{ + return nfront > 0; +}; + +bool XG_algorithm::state::freeze_front_overtaken() +{ + return Zdf > 0.0 && Zdt >= Zdf; +}; + +bool XG_algorithm::state::thaw_front_overtaken() +{ + return Zdt > 0.0 && Zdf >= Zdt; +}; + +XG_algorithm::state& XG_algorithm::state::set_layer_moisture_maximums(const XG_algorithm::params& P) +{ + //XG_algorithm xg(0.0,0.0,0.0,*this,P); + //xg.printdebug(); + double sum = 0.0; + std::vector depths = P.getdepths(); + for (double val : depths) + { + sum += val; + } + if (sum < this->Zdf || sum < this->Zdt) + { + //TODO Add A CHM exception to say that the total soil depth is less than initial Zdt,Zdf + } + double rechrmax = P.soil_rechr_max; + double soilmax = P.soil_moist_max; + //xg.printdebug(); + for (size_t layer = 0; layer < P.N_Soil_layers; ++layer) + { + + //xg.printdebug(); + this->XG_max[layer] = 0; + //xg.printdebug(); + this->XG_max[layer] = P.por[layer] * depths[layer] * 1000.0; + //xg.printdebug(); + + this->theta[layer] = P.theta_default[layer]; + //xg.printdebug(); + if (rechrmax > 0.0) + { + //xg.printdebug(); + if (rechrmax > this->XG_max[layer]) + { + //xg.printdebug(); + this->XG_rechr_d += depths[layer]; + this->rechr_fract[layer] = 1.0; + rechrmax -= this->XG_max[layer]; + } + else + { + //xg.printdebug(); + const double amount = rechrmax / this->XG_max[layer]; + this->rechr_fract[layer] = rechrmax / this->XG_max[layer]; + + this->XG_rechr_d += depths[layer] * amount; + const double amount_remaining = 1.0 - amount; + //xg.printdebug(); + if (soilmax >= this->XG_max[layer]*amount_remaining) + { + this->moist_fract[layer] = amount_remaining; + soilmax -= this->XG_max[layer] * amount_remaining; + this->XG_moist_d += depths[layer]; + } + else + { + this->moist_fract[layer] = (soilmax - rechrmax) / this->XG_max[layer]; + const double used = this->rechr_fract[layer] + this->moist_fract[layer]; + this->default_fract[layer] = 1.0 - used; + this->XG_moist_d += this->XG_rechr_d + depths[layer] * used; + soilmax = 0.0; + } + rechrmax = 0.0; + } + } + else if (soilmax > 0.0) + { + //xg.printdebug(); + if (soilmax >= this->XG_max[layer]) { + this->XG_moist_d += depths[layer]; + this->moist_fract[layer] = 1.0; + soilmax -= this->XG_max[layer]; + } + else + { + const double amount = soilmax / this->XG_max[layer]; + this->XG_moist_d += depths[layer] * amount; + this->moist_fract[layer] = amount; + this->default_fract[layer] = 1.0 - amount; + soilmax = 0.0; + } + } + else + { + this->default_fract[layer] = 1.0; + } + } + + //xg.printdebug(); + if (rechrmax != 0.0 || soilmax != 0.0) + { + // put CHM exception here + } + return *this; +}; + +XG_algorithm::state& XG_algorithm::state::set_thermal_conductivities(const XG_algorithm::params& P,const double& soil_moist, const double& soil_rechr) +{ + Th_low = 1; + Fz_low = 1; + check_XG_moist = 0.0; + XG_moist.resize(P.N_Soil_layers,0.0); + + // Process each soil layer + for (size_t layer = 0; layer < P.N_Soil_layers; ++layer) + { + // Calculate moisture content + if (P.soil_moist_max > 0.0) + { // handle soil_moist_max = 0.0 (slough case) + XG_moist[layer] = + rechr_fract[layer] * XG_max[layer] * (soil_rechr / P.soil_rechr_max) + + moist_fract[layer] * XG_max[layer] * + (soil_moist - soil_rechr) / (P.soil_moist_max - P.soil_rechr_max); + } + else + { + XG_moist[layer] = 0.0; + } + + // Update moisture checks and adjustments + check_XG_moist += XG_moist[layer]; + XG_moist[layer] += + default_fract[layer] * XG_max[layer] * P.theta_default[layer]; + + + // Calculate and validate theta + theta[layer] = XG_moist[layer] / XG_max[layer]; + theta[layer] = (theta[layer] < P.theta_min) ? P.theta_min : theta[layer]; + //if (theta[layer] <= P.theta_min) + //{ + // theta[layer] = P.theta_min; // enforce minimum value + //} + + // Convert to water content (kg/m³) + layer_h2o[layer] = theta[layer] * P.por[layer] * 1000.0; + // Update thermal conductivities + if (P.k_update) + { // dynamic update mode + ftc[layer] = (ftc_contents[layer] == 1) + ? get_ttc(layer,P) + : get_ftc(layer,P); + + ttc[layer] = (ttc_contents[layer] == 1) + ? get_ftc(layer,P) + : get_ttc(layer,P); + + } + else + { + ftc[layer] = get_ftc(layer,P); + ttc[layer] = get_ttc(layer,P); + ftc_contents[layer] = 0; + ttc_contents[layer] = 0; + } + } + return *this; +}; + +XG_algorithm::state& XG_algorithm::state::set_freezethaw_ratios(const XG_algorithm::params& P) +{ + for (size_t layer = 1; layer < P.N_Soil_layers; ++layer) + { + pf[layer] = std::sqrt( + (ftc[layer-1] / layer_h2o[layer-1]) / + (ftc[layer] / layer_h2o[layer]) + ); + + pt[layer] = std::sqrt( + (ttc[layer-1] / layer_h2o[layer-1]) / + (ttc[layer] / layer_h2o[layer]) + ); + }; + + return *this; +}; + +void XG_algorithm::state::accumulate_degree_days(const double& surface_temp) +{ + B += surface_temp / P.time_step_per_day; + TrigAcc += B; + + t_trend -= t_trend / 192; + t_trend += B/192; +}; +//void XG_algorithm::init_freezethaw_front_depths(const double& Zdf_init, const double& Zdt_init, const double& Zpf_init) +//{ +// if (Zdf_init > 0.0 && Bfr == 0.0) +// find_freeze_D(Zdf_init); +// +// if (Zdt_init > 0.0 && Zdt_init < Zpf_init && Bth == 0.0) +// find_thaw_D(Zdt_init); +//}; + +//void XG_algorithm::run() +//{ +// +// if (DTO.is_day_start(DTO)) +// B = 0.0; +// +// B += DTO.t_surface; +// +// if (STO.is_day_start(DTO)) +// { +// // idle to freeze/thaw or freeze/thaw to idle +// +// // Keep T_total within +/- T_threshold +// // TODO might not be needed, might be diagnositic like B +// if (T_total > T_threshold) +// T_total = T_threshold; +// else if (T_total < -T_threshold) +// T_total = -T_threshold; +// +// if (std::abs(T_total) >= T_threshold/2.0 && freeze_thaw_state == 0) +// { +// if (freeze_front_depth < 0.0 || num_fronts) +// freeze_thaw_state = 1; +// else +// freeze_thaw_state = -1; +// +// temp_trend = 0.0; +// } +// +// if (freeze_thaw_state && std::abs(T_total) >= T_threshold/2.0 && std::abs(temp_trend) > 0.0) +// { +// freeze_thaw_state = 0; +// +// if (freeze_front_depth > 0.0 && thaw_front_depth > 0.0) +// { +// if ( is_freezing() ) +// switch_to_idle(thaw_front_depth); +// else +// switch_to_idle(-freeze_front_depth); +// } +// } +// +// // Calculate thermal conductivities +// +// // Issues +// +// thaw_low_layer = 1; +// freeze_low_layer = 1; +// +// double thaw_k_T_rechr = get_k_T(DTO.soil_rechr_storage); +// double freeze_k_T_rechr = ; +// double thaw_k_T_lower = ; +// double freeze_k_T_lower = ; +// +// for (int ii = 0; ii < N_layers; ++ii) +// { +// if (DTO.soil_storage_max > 0.0) +// { +// +// } +// +// } +// } +// +// +// +// +// +// +// +//}; +// +//bool XG_algorithm::is_freezing() +//{ +// if (freeze_thaw_state == -1 && T_total > 0.0 && temp_trend > 0.0) +// return true; +// else if (freeze_thaw_state == 1 && T_total < 0.0 && temp_trend < 0.0) +// return false; +//} +// +//void XG_algorithm::switch_to_idle(double front_depth) +//{ +// push_front(front_depth);// TODO write this function +// // TODO in crhm both parts modify Zd_front_array (something I've decided not to include in XG) come back if necessary +// if (front_depth > 0.0 && thaw_front_depth > freeze_front_depth) +// { +// thaw_front_depth = 0.0; +// thaw_degree_days = 0.0; +// } +// else (front_depth < 0.0 && freeze_front_depth > thaw_front_depth) +// { +// freeze_front_depth = 0.0; +// freeze_degree_days = 0.0; +// } +// +//}; +// +// +// +// diff --git a/src/modules/soil_submodules/XG_algorithm.hpp b/src/modules/soil_submodules/XG_algorithm.hpp new file mode 100644 index 00000000..3a8541be --- /dev/null +++ b/src/modules/soil_submodules/XG_algorithm.hpp @@ -0,0 +1,279 @@ +// This code was written as an exact copy of the CRHM module ClassXG +// Therefore, any bugs in ClassXG must also be present here +// Tests were written by examining the original code to 'guess' at behaviour +#include "I_freeze_thaw_depths.hpp" +#include +#include + +class XG_algorithm : public I_freeze_thaw_depths +{ +public: + class state; + class params; + XG_algorithm(const double& t, const double& _moist, + const double& _rechr, state& _S, const params& _P); + + ~XG_algorithm() {}; + + virtual void run() override; + void init_freezethaw_degreedays(const double& Zdf_init, const double& Zdt_init, const double& Zpf_init); + + static constexpr double ko = 0.21; // W/(m K) organic material + static constexpr double km = 2.50; // W/(m K) mineral + static constexpr double ka = 0.025; // W/(m K) air + static constexpr double ki = 2.24; // W/(m K) ice 2.24 + static constexpr double kw = 0.57; // W/(m K) water 0.57 + +private: + const double& surface_temp; + const double& soil_moist; + const double& soil_rechr; + state& S; + const params& P; + std::vector depths; + + void freeze(void); + void thaw(void); + double stefan_equation(double& surface_index, + double& thermal_conductivity, size_t& layer); + double Interpolated_ttc(double Za, size_t layer); + double Interpolated_ftc(double Za, size_t layer); + void push_front(double D); + void find_thaw_D(double dt); + void find_freeze_D(double df); + +public: + class state + { + public: + // Depth-related + double Zdf; // (m) depth of freezing front + double Zdt; // (m) depth of thawing front + std::vector Zd_front; // (m) depths of all freezing/thawing fronts (thaw: +, freeze: -) + size_t Th_low; // lowest thawed layer + size_t Fz_low; // lowest frozen layer + size_t nfront; // number of freezing/thawing fronts + + // Degree-day counters + double Bfr; // (ºC*d) freeze degree days + double Bth; // (ºC*d) thaw degree days + + // Layer ratios + std::vector pf; // soil layers freezing ratios + std::vector pt; // soil layers thawing ratios + + // Thermal properties + std::vector ttc; // (W/(m*K)) thawing thermal conductivity + std::vector ftc; // (W/(m*K)) freezing thermal conductivity + std::vector tc_composite; // (W/(m*K)) composite ftc/ttc + std::vector tc_composite2; // (W/(m*K)) composite2 ftc/ttc + + // Moisture/water content + std::vector theta; // (m³/m³) layer theta + std::vector layer_h2o; // (kg/m³) layer water content + double XG_moist_d; // (m) layer depth of soil moisture in XG + double XG_rechr_d; // (m) layer depth of soil recharge in XG + std::vector XG_max; // (mm) layer max soil moisture + std::vector XG_moist; // (mm) layer moisture content + double check_XG_moist; // (mm) sum of XG soil_moist + + // Local/temporary + double B; // (ºC*d) interval degree-day sum + double TrigAcc; // (ºC*d) freeze/thaw cycle detection + int TrigState; // 1/0/-1 → thaw/idle/freeze + std::vector ttc_contents; // 0/1 → thaw/freeze + std::vector ftc_contents; // 0/1 → freeze/thaw + double t_trend; // (°C) temperature long-term trend + + bool is_newday = false; + bool last_step_new_day = false; + + // Fractions + std::vector rechr_fract; // fraction of layer (soil_rechr_max) + std::vector moist_fract; // fraction of layer (soil_moist_max) + std::vector default_fract; // fraction of layer (theta_default) + // Constructor that initializes vector sizes + // + + // TODO: refactor P as a unique pointer which is created by the module, then ownership is moved to B + const params& P; + + explicit state(size_t N,const params& _P) : + Zd_front(N,0.0), + pf(N,0.0), + pt(N,0.0), + ttc(N,0.0), + ftc(N,0.0), + tc_composite(N,0.0), + tc_composite2(N,0.0), + theta(N,0.0), + layer_h2o(N,0.0), + XG_max(N,0.0), + XG_moist(N,0.0), + ttc_contents(N,0.0), + ftc_contents(N,0.0), + rechr_fract(N,0.0), + moist_fract(N,0.0), + default_fract(N,0.0), + P(_P) + { + assert(N > 0 && "State requires a positive number of layers"); + // Initialize other members + Zdf = 0.0; + Zdt = 0.0; + Th_low = 1; + Fz_low = 1; + nfront = 0; + Bfr = 0.0; + Bth = 0.0; + XG_moist_d = 0.0; + XG_rechr_d = 0.0; + check_XG_moist = 0.0; + B = 0.0; + TrigAcc = 0.0; + TrigState = 0; + t_trend = 0.0; + } + double get_ftc(size_t layer, const params& P); + double get_ttc(size_t layer, const params& P); + + //void set_XG_max(params& P); + //void set_theta(params& P); + state& set_layer_moisture_maximums(const params& P); + state& set_thermal_conductivities(const params& P, + double const& soil_moist, double const& soil_rechr); + state& set_freezethaw_ratios(const params& P); + + bool freezing(); + bool thawing(); + bool net_negative_degree_days(); + void reset_degree_day_counter(); + void determine_freeze_thaw_idle(); + bool thaw_front_overtaken(); + bool freeze_front_overtaken(); + void merge_freeze_fronts(XG_algorithm* XG); + void merge_thaw_fronts(XG_algorithm* XG); + bool exist_excess_fronts(); + void push_front(double D); + void accumulate_degree_days(const double& surface_temp); + double last_front(void); + double pop_front(void); + + }; + + class params + { + private: + std::vector depths; // (m) soil layer thicknesses + public: + // Core parameters + double Trigthrhld; // (ºC*d) Trigger reference level + const std::vector& getdepths() const + { return depths; }; + std::vector por; // soil porosity + size_t N_Soil_layers; // number of soil layers (≤ nlay) + std::vector theta_default; // (m³/m³) default theta + double theta_min; // (m³/m³) minimum theta + std::vector soil_solid_km; // (W/(m*K)) dry soil conductivity + std::vector soil_solid_km_ki; // (W/(m*K)) saturated frozen conductivity + std::vector soil_solid_km_kw; // (W/(m*K)) saturated unfrozen conductivity + double SWE_k; // (W/(m*K)) snow thermal conductivity UNUSED IN CRHM + //Next 3 are initial conditions + //const double Zdf_init; // (m) initial freezing front depth + //const double Zdt_init; // (m) initial thawing front depth + const double Zpf_init; // (m) initial permafrost depth + bool freeze_kw_ki_update; // update kw→ki behind freeze front + bool thaw_ki_kw_update; // update ki→kw behind thaw front + size_t k_update; // 0=never, 1=post-layer, 2=continuous + double soil_rechr_max; // (mm) max recharge zone capacity + double soil_moist_max; // (mm) max rooting zone capacity + double time_step_per_day; + bool calc_conductivity; + + ~params() {}; + params( + const std::vector d, const double& t, const std::vector p, + const size_t& n, const std::vector td, const double& tm, const std::vector skm, + const std::vector ski, const std::vector skw, const double swk, const double zpf, + const size_t& fku, const size_t& tku, const size_t& ku, const double& srm, const double& smm, + const double& tspd, const bool& cc + ) : + depths(d), + Trigthrhld(t), + por(p), + N_Soil_layers(n), + theta_default(td), + theta_min(tm), + soil_solid_km(skm), + soil_solid_km_ki(ski), + soil_solid_km_kw(skw), + SWE_k(swk), + Zpf_init(zpf), + freeze_kw_ki_update(fku), + thaw_ki_kw_update(tku), + k_update(ku), + soil_rechr_max(srm), + soil_moist_max(smm), + time_step_per_day(tspd), + calc_conductivity(cc) + { + //CHM + }; + }; + // Variation #1 parameters + //const double& n_factor_a; // surface-to-air temp ratio + //const double& n_factor_b; // surface-to-air temp ratio + //const double& n_factor_c; // surface-to-air temp ratio +}; + +//class StateBuilder { +//public: +// explicit StateBuilder(int num_layers) { +// state_ = std::make_unique(); +// initialize_vectors(num_layers); +// } +// +// StateBuilder& set_initial_freezethaw_depths(const double Zdf, const double Zdt) +// { +// state_->Zdf = Zdf; +// state_->Zdt = Zdt; +// +// return *this; +// }; + +//XG_algorithm::state& XG_algorithm::state::check_initial_Zdf_depths(const std::vector& depths) +//{ +// double sum = 0.0; +// for (double val : depths) +// { +// sum += val; +// } +// +// if (sum < this->Zdf || sum < this->Zdt) +// { +// // TODO add exception +// } +// +// return *this; +//}; + +//XG_algorithm::state& XG_algorithm::state::set_XG_max(std::vector por,std::vector depths) +//{ +// +// for (auto&& [m,p,d] : std::views::zip(state_->XG_max,por,depths)) +// m = p * d * 1000.0; +// +// return *this; +// +//}; +// +//XG_algorithm::state& XG_algorithm::state::set_theta(std::vector theta_default) +//{ +// state_->theta = theta_default; +// +// return *this; +//}; + + + + diff --git a/src/modules/soil_submodules/soil_DTO.hpp b/src/modules/soil_submodules/soil_DTO.hpp new file mode 100644 index 00000000..e7fcd7ed --- /dev/null +++ b/src/modules/soil_submodules/soil_DTO.hpp @@ -0,0 +1,96 @@ +#pragma once + +struct shared_DTO +{ + // shared between more than one submodule + double soil_storage_max = 0.0; + double swe = 0.0; //in - snobal TODO should this include input from PBSM3D? + double soil_storage = 0.0; + double soil_rechr_storage = 0.0; + double soil_rechr_max = 0.0; + double depression_storage = 0.0; + double actual_ET = 0.0; +}; + +struct soil_ET_DTO : virtual shared_DTO +{ + + double actual_soil_ET = 0.0; //out + int ground_cover_type = 0; + int soil_type_rechr = 0; + int soil_type_lower = 0; + bool is_lake = false; + //REMOVED and replaced by variable (see above) (MAy 2025) + //virtual bool is_lake(soil_ET_DTO& DTO) = 0; + +}; + +struct two_layer_DTO : virtual shared_DTO +{ + double thaw_front_depth = 0.0; //in - XG + double freeze_front_depth = 0.0; //in - XG + double freeze_thaw_first_front = 0.0; //in - XG + double infil = 0.0; //in - infil_all + double runoff = 0.0; //in - infil_all + double routing_residual = 0.0; //in - NetRoute + + double condensation = 0.0; //out (and used) + double soil_excess_to_runoff = 0.0; //out (and used) + double runoff_to_depression = 0.0; //out + double soil_excess_to_gw = 0.0; //out (and used) + double ground_water_out = 0.0; //out (and used) + double soil_to_ssr = 0.0; //out (and used) + double rechr_to_ssr = 0.0; + double excess = 0.0; + + // This parameter, if false, disables runoff produces from infiltration + // Or, ignores thaw fraction in soil produces by thaw/freeze fronts + bool allow_runoff_from_infiltration; + // All K_ stuff are in + double K_soil_to_gw = 0.0; // All these K's could be swapped to held from previous and could be global/system wide + double K_rechr_to_ssr = 0.0; + double K_lower_to_ssr = 0.0; + double K_detention_to_runoff = 0.0; + double K_depression_to_ssr = 0.0; + double K_depression_to_gw = 0.0; + double K_ground_water_out = 0.0; + double Ksaturated_rechr = 0.0; + double Ksaturated_lower = 0.0; + double Ksaturated_ground_water = 0.0; + double Ksaturated_snow = 0.0; + double Ksaturated_organic = 0.0; + + // held from previous + double thaw_fraction_rechr = 0.0; + double thaw_fraction_lower = 0.0; + double porosity = 0.0; + bool excess_to_ssr = true; + double detention_max = 0.0; + double detention_storage = 0.0; + double detention_snow_max = 0.0; + double detention_organic_max = 0.0; + double detention_snow_init = 0.0; + double detention_organic_init = 0.0; + double depression_max = 0.0; + double depression_to_gw = 0.0; + double ground_water_storage = 0.0; + double ground_water_max = 0.0; + double snow_covered_threshold = 0.0; + + // K_estimate + double pore_size_dist = 0.0; + double pore_size_dist_organic = 0.0; + double soil_index = 0.0; + double local_slope = 0.0; + double snow_grain_diameter = 0.0; // divided by 1000 in K_estimate, likely converting from mm to m, this var is mm + double snow_density = 0.0; + + virtual int get_dt() = 0; + virtual bool get_new_day() = 0; + //virtual bool is_day_start(two_layer_DTO& DTO) = 0; + //virtual void call_model_error(two_layer_DTO& DTO) = 0; +}; + +struct main_DTO : two_layer_DTO, soil_ET_DTO +{ +}; diff --git a/src/modules/soil_submodules/soil_ET.cpp b/src/modules/soil_submodules/soil_ET.cpp new file mode 100644 index 00000000..5666f270 --- /dev/null +++ b/src/modules/soil_submodules/soil_ET.cpp @@ -0,0 +1,153 @@ +#include "soil_ET.hpp" +#include +soil_ET::soil_ET(soil_ET_DTO& _DTO) : DTO(_DTO) +{ + +}; + +soil_ET::~soil_ET() +{ + +}; + +void soil_ET::run() +{ + + DTO.actual_soil_ET = 0.0; + + if (DTO.swe > 0.0) + return; + + std::cout << "in soil ET" << std::endl; + double available_to_evap = DTO.actual_ET; + if (DTO.depression_storage + DTO.soil_storage > 0.0) + { + available_to_evap *= DTO.depression_storage / + (DTO.depression_storage + DTO.soil_storage); + } + else + available_to_evap = 0.0; + + if (DTO.depression_storage > 0.0 && available_to_evap > 0.0) + { + if (DTO.depression_storage > available_to_evap) + DTO.depression_storage = std::max(0.0,DTO.depression_storage - available_to_evap); + else + { + available_to_evap = DTO.depression_storage; + DTO.depression_storage = 0.0; + } + DTO.actual_ET = available_to_evap; + } + else + available_to_evap = 0.0; + + std::cout << "in soil ET" << std::endl; + available_to_evap = DTO.actual_ET - available_to_evap; + + if (available_to_evap > 0.0 && DTO.soil_storage > 0.0 && DTO.ground_cover_type > 0) + { + double percent_available_lower; + double percent_available_rechr; + double ET_lower; + double ET_rechr; + double soil_lower_storage = DTO.soil_storage - DTO.soil_rechr_storage; + double soil_lower_max = DTO.soil_storage_max - DTO.soil_rechr_max; + + if ( soil_lower_max > 0.0 ) // soil_lower > 0.0 + percent_available_lower = soil_lower_storage / soil_lower_max; + else + percent_available_lower = 0.0; + + percent_available_rechr = DTO.soil_rechr_storage / DTO.soil_rechr_max; + + ET_rechr = available_to_evap; + + ET_rechr = set_ET_layer(ET_rechr,percent_available_rechr,DTO.soil_type_rechr); + + if (ET_rechr > available_to_evap) + { + ET_lower = 0.0; + ET_rechr = available_to_evap; + } + else + ET_lower = available_to_evap - ET_rechr; + + if (ET_lower > 0.0) + { + ET_lower = set_ET_layer(ET_lower, percent_available_lower, DTO.soil_type_lower); + } + + double ET = 0.0; + std::cout << "in soil ET" << std::endl; + switch (DTO.ground_cover_type) + { + case 0: // bare soil, no evap + break; + case 1: // recharge layer only, crops + if (ET_rechr > DTO.soil_rechr_storage) + { + ET = DTO.soil_rechr_storage; + DTO.soil_rechr_storage = 0.0; + } + else + { + DTO.soil_rechr_storage -= ET_rechr; + ET = ET_rechr; + } + DTO.soil_storage -= ET_rechr; + break; + case 2: // all soil moisture, grasses & shrubs + if (ET_rechr + ET_lower >= DTO.soil_storage) + { + ET = DTO.soil_storage; + DTO.soil_storage = 0.0; + DTO.soil_rechr_storage = 0.0; + } + else + { + ET = ET_rechr + ET_lower; + DTO.soil_storage -= ET; + + DTO.soil_rechr_storage = std::max(DTO.soil_rechr_storage - ET_rechr, 0.0); + } + break; + } + + DTO.actual_soil_ET += ET; + std::cout << "in soil ET" << std::endl; + if (DTO.is_lake) + { + std::cout << "out" << std::endl; + DTO.actual_soil_ET = DTO.actual_ET; + } + }; +}; + +double soil_ET::set_ET_layer(double et, double& percent, int& type) +{ + switch (type) + { + case 1: // sandy soil + if (percent < 0.25) + et = 0.5 * percent * et; + break; + case 2: // loam soil + if (percent < 0.5) + et = percent * et; + break; + case 3: // clay soil + if (percent <= 0.33) + et = 0.5 * percent * et; + else if (percent < 0.67) + et = percent * et; + break; + default: //organic soil + // do nothing, use default above + break; + } + + return et; +}; + + diff --git a/src/modules/soil_submodules/soil_ET.hpp b/src/modules/soil_submodules/soil_ET.hpp new file mode 100644 index 00000000..67714a88 --- /dev/null +++ b/src/modules/soil_submodules/soil_ET.hpp @@ -0,0 +1,19 @@ +#pragma once + +#include "soil_base.hpp" +#include "soil_DTO.hpp" +class soil_ET : public soil_base +{ +public: + soil_ET(soil_ET_DTO& _DTO); + ~soil_ET(); + + void run() override; + +private: + + soil_ET_DTO& DTO; + + double set_ET_layer(double et, double& percent, int& type); +}; + diff --git a/src/modules/soil_submodules/soil_base.hpp b/src/modules/soil_submodules/soil_base.hpp new file mode 100644 index 00000000..79ba4d00 --- /dev/null +++ b/src/modules/soil_submodules/soil_base.hpp @@ -0,0 +1,11 @@ +#pragma once + +#include "submodule_base.hpp" +#include +class soil_base : public submodule_base +{ +public: + virtual void run() = 0; + +}; + diff --git a/src/modules/soil_submodules/soil_classes.cpp b/src/modules/soil_submodules/soil_classes.cpp new file mode 100644 index 00000000..633f626e --- /dev/null +++ b/src/modules/soil_submodules/soil_classes.cpp @@ -0,0 +1,281 @@ +#include "soil_classes.hpp" +void soil_class_base::_push_excess_down(double& layer_storage, double& layer_max, + double& layer_down) +{ + double excess = layer_storage - layer_max; + + layer_storage = layer_max; + + layer_down += excess; +}; + +void initializer::zero_single_step_vars() +{ + d.condensation = 0.0; + d.soil_excess_to_runoff = 0.0; + d.soil_excess_to_gw = 0.0; + d.ground_water_out = 0.0; + d.soil_to_ssr = 0.0; + d.rechr_to_ssr = 0.0; + d.excess = 0.0; + d.runoff_to_depression = 0.0; + d.routing_residual = 0.0; //For now, no routing so zeroed here, TODO +}; + +void initializer::layer_thaw_fraction() +{ + //rerun once daily + if (!d.get_new_day()) + return; + + + d.thaw_fraction_rechr = 0.0; + d.thaw_fraction_lower = 0.0; + + // reference so intention is more clear + // If we are allowing runoff from infiltration we are also allowing the soil to freeze + bool& allow_soil_to_freeze = d.allow_runoff_from_infiltration; + double rechr_depth = 0.0; + double soil_depth = 0.0; + + // TODO this might not be right, but the above is what CRHM does + // depth = storage / porosity, because porosity = storage / depth; + if (d.soil_storage_max ==0.0) + { + d.thaw_fraction_lower = 0.0; + d.thaw_fraction_rechr = 0.0; + return; + } + + if (d.porosity > 0.0) + { + // TODO, porosity in rechr is the same as lower + rechr_depth = d.soil_rechr_max / d.porosity / 1000.0; + soil_depth = d.soil_storage_max / d.porosity / 1000.0; + } + else + { + d.thaw_fraction_rechr = 0.0; + d.thaw_fraction_lower = 0.0; + return; + }; + + if (!allow_soil_to_freeze || d.freeze_thaw_first_front == 0.0) + { + d.thaw_fraction_rechr = 1.0; + d.thaw_fraction_lower = 1.0; + } + else + { + d.thaw_fraction_rechr = 0.0; + d.thaw_fraction_lower = 0.0; + } + + if (d.soil_storage_max > 0.0 && allow_soil_to_freeze && + d.freeze_thaw_first_front > 0.0) + { + + // TODO Verify this calculation + if (d.thaw_front_depth < rechr_depth) + { + d.thaw_fraction_rechr = d.thaw_front_depth / rechr_depth; + } + else if ( d.thaw_front_depth < soil_depth) + { + d.thaw_fraction_rechr = 1.0; + d.thaw_fraction_lower = (d.thaw_front_depth - rechr_depth) + / (soil_depth - rechr_depth); + } + else + { + d.thaw_fraction_rechr = 1.0; + d.thaw_fraction_lower = 1.0; + } + + } +}; + +void condensator::set() +{ + if (d.actual_ET < 0.0 && d.swe == 0.0) + { + d.condensation = -1.0 * d.actual_ET; + d.actual_ET = 0.0; + } +}; + +void infiltrator::distribute() +{ + if (d.soil_storage_max > 0.0) + { + double soil_lower_storage = d.soil_storage - d.soil_rechr_storage; + + double potential = d.infil + d.condensation; + + double possible = d.thaw_fraction_rechr * (d.soil_rechr_max - d.soil_rechr_storage); + if (possible > potential || !d.allow_runoff_from_infiltration) + possible = potential; + else + d.soil_excess_to_runoff = potential - possible; + d.soil_rechr_storage += possible; + + if (d.soil_rechr_storage > d.soil_rechr_max) + _push_excess_down(d.soil_rechr_storage,d.soil_rechr_max,soil_lower_storage); + + d.soil_storage = soil_lower_storage + d.soil_rechr_storage; + + if (d.soil_storage > d.soil_storage_max) + _push_excess_down(d.soil_storage,d.soil_storage_max,d.soil_excess_to_gw); + + if (d.swe == 0.0) // if there is no snowcover + { + d.rechr_to_ssr = d.soil_rechr_storage / d.soil_rechr_max * d.K_rechr_to_ssr * d.thaw_fraction_rechr; + d.rechr_to_ssr = std::min(d.rechr_to_ssr,d.soil_rechr_storage * d.thaw_fraction_rechr); + + d.soil_rechr_storage = std::max(0.0, d.soil_rechr_storage - d.rechr_to_ssr); + + d.soil_storage -= d.rechr_to_ssr; + d.soil_to_ssr = d.rechr_to_ssr; + + } + + if (d.soil_excess_to_gw > d.K_soil_to_gw * d.thaw_fraction_lower) + { + double excess_to_gw_max = d.K_soil_to_gw * d.thaw_fraction_lower; + _push_excess_down(d.soil_excess_to_gw,excess_to_gw_max,d.excess); + } + + // Line 607 of SoilX crhmcommetns branch, comment says upper layer but code says lower-layer, ask logan about this. + if (d.excess_to_ssr && d.excess > 0.0) + { + double excess_to_ssr_max = d.excess * (1.0 - d.thaw_fraction_lower); + _push_excess_down(d.excess,excess_to_ssr_max,d.soil_to_ssr); + } + + + } + else + { + d.excess = d.infil + d.condensation; + } +}; + +void detention_layer::manage() +{ + double face_area = 1.0; // Later will be pulled from face, but since routine_residual is always zero, ignoring it here. + d.soil_excess_to_runoff += d.runoff + d.excess + d.routing_residual / face_area; // routing_residual comes from the crhm varaible redirected_residual which has units of mm*km^2/int (not sure why), so face_area is there for now for consistency. + + if (d.soil_excess_to_runoff > 0.0) + { + if (d.swe == 0.0) + d.detention_max = d.detention_organic_max; + else + d.detention_max = d.detention_snow_max; + + double detention_space = d.detention_max - d.detention_storage; + + if (detention_space > 0.0) + { + if (d.soil_excess_to_runoff > detention_space) + { + d.soil_excess_to_runoff = std::max(0.0,d.soil_excess_to_runoff - detention_space); + d.detention_storage += detention_space; + } + else + { + d.detention_storage += d.soil_excess_to_runoff; + d.soil_excess_to_runoff = 0.0; + } + } + } + + if (d.detention_storage > 0.0 && d.K_detention_to_runoff > 0.0) + { + double transfer = std::min(d.detention_storage,d.K_detention_to_runoff); + d.soil_excess_to_runoff += transfer; + d.detention_storage -= transfer; + + if (d.detention_storage < 0.0001) // from CRHM, for safety and to drop any Floating-point errors + d.detention_storage = 0.0; + } +}; + +void depression_layer::manage() +{ + if (d.soil_excess_to_runoff > 0.0 && d.depression_max > 0.0) + { + double exponent = -1.0 * std::min(12.0,d.soil_excess_to_runoff / d.depression_max); + + double depression_space = (d.depression_max - d.depression_storage) * (1 - exp(exponent)); + + if (d.soil_storage_max == 0.0) + depression_space = d.depression_max - d.depression_storage; + + if (depression_space > 0.0) + { + if (d.soil_excess_to_runoff > depression_space) + { + d.soil_excess_to_runoff = std::max(0.0,d.soil_excess_to_runoff - depression_space); + d.depression_storage += depression_space; + d.runoff_to_depression += depression_space; + // TODO add total tracker + } + else + { + d.depression_storage += d.soil_excess_to_runoff; + // TODO add total tracker + d.soil_excess_to_runoff = 0.0; + } + + } + } + + if (d.depression_storage > 0.0 && d.K_depression_to_gw > 0.0) + { + double amount_to_move = std::min(d.depression_storage,d.K_depression_to_gw); + d.depression_to_gw += amount_to_move; + d.depression_storage -= amount_to_move; + if (d.depression_storage < 0.0) // floatig point error safety? + d.depression_storage = 0.0; + } +}; + +void groundwater_layer::manage() +{ + d.soil_excess_to_gw += d.depression_to_gw; + d.depression_to_gw = 0.0; + d.ground_water_storage += d.soil_excess_to_gw; + + if (d.ground_water_storage > d.ground_water_max) + { + _push_excess_down(d.ground_water_storage,d.ground_water_max,d.ground_water_out); + } + + if (d.ground_water_max > 0.0) // divide by zero safety + { + double spilled = d.ground_water_storage / d.ground_water_max * d.K_ground_water_out; + d.ground_water_storage -= spilled; + d.ground_water_out += spilled; + } +}; + +void subsurface_runoff::manage() +{ + if (d.depression_storage > 0.0 && d.K_depression_to_ssr > 0.0) + { + double amount_to_move = std::min(d.depression_storage,d.K_depression_to_ssr); + d.soil_to_ssr += amount_to_move; + d.depression_storage -= amount_to_move; + if (d.depression_storage < 0.0) + d.depression_storage = 0.0; + } + + if (d.K_lower_to_ssr > 0.0) + { + double available = d.soil_storage - d.soil_rechr_storage; + double lower_unfrozen = d.K_lower_to_ssr * d.thaw_fraction_lower; + double amount_to_move = std::min(lower_unfrozen,available); + d.soil_storage -= amount_to_move; + d.soil_to_ssr += amount_to_move; + } +}; diff --git a/src/modules/soil_submodules/soil_classes.hpp b/src/modules/soil_submodules/soil_classes.hpp new file mode 100644 index 00000000..a33e4c77 --- /dev/null +++ b/src/modules/soil_submodules/soil_classes.hpp @@ -0,0 +1,86 @@ +#pragma once +#include "soil_DTO.hpp" +#include +#include +class soil_class_base +{ +public: + soil_class_base(two_layer_DTO& DTO) : d(DTO) {}; + + ~soil_class_base() {}; +protected: + two_layer_DTO& d; + + void _push_excess_down(double& layer_storage, double& layer_max, + double& layer_down); + +}; + +class initializer : public soil_class_base +{ +public: + initializer(two_layer_DTO& DTO) : soil_class_base(DTO) {}; + ~initializer() {}; + + void zero_single_step_vars(); + void layer_thaw_fraction(); + +}; + +class condensator : public soil_class_base +{ +public: + condensator(two_layer_DTO& DTO) : soil_class_base(DTO) {}; + ~condensator() {}; + + void set(); + +}; + +class infiltrator : public soil_class_base +{ +public: + infiltrator(two_layer_DTO& DTO) : soil_class_base(DTO) {}; + ~infiltrator() {}; + + void distribute(); + +}; + +class detention_layer : public soil_class_base +{ +public: + detention_layer(two_layer_DTO& DTO) : soil_class_base(DTO) {}; + ~detention_layer() {}; + + void manage(); + +}; + +class depression_layer : public soil_class_base +{ +public: + depression_layer(two_layer_DTO& DTO) : soil_class_base(DTO) {}; + ~depression_layer() {}; + void manage(); + +}; + +class groundwater_layer : public soil_class_base +{ +public: + groundwater_layer(two_layer_DTO& DTO) : soil_class_base(DTO) {}; + ~groundwater_layer() {}; + void manage(); + +}; + +class subsurface_runoff : public soil_class_base +{ +public: + subsurface_runoff(two_layer_DTO& DTO) : soil_class_base(DTO) {}; + ~subsurface_runoff() {}; + + void manage(); + +}; diff --git a/src/modules/soil_submodules/soil_two_layer.cpp b/src/modules/soil_submodules/soil_two_layer.cpp new file mode 100644 index 00000000..3d0a8f7c --- /dev/null +++ b/src/modules/soil_submodules/soil_two_layer.cpp @@ -0,0 +1,331 @@ +#include "soil_two_layer.hpp" + +soil_two_layer::soil_two_layer(two_layer_DTO& DTO,I_K_estimate& estimate) : + init(DTO), + condensation(DTO), + infiltrate(DTO), + detention(DTO), + depression(DTO), + groundwater(DTO), + ssr(DTO), + k_estimator(estimate) +{ +}; +void soil_two_layer::run() +{ + + init.zero_single_step_vars(); + init.layer_thaw_fraction(); + + condensation.set(); + + infiltrate.distribute(); + + detention.manage(); + + depression.manage(); + + groundwater.manage(); + + ssr.manage(); + + k_estimator.run(); +}; + +//void soil_two_layer::initialize_single_step_vars() +//{ +// DTO.condensation = 0.0; +// DTO.soil_excess_to_runoff = 0.0; +// DTO.soil_excess_to_gw = 0.0; +// DTO.ground_water_out = 0.0; +// DTO.soil_to_ssr = 0.0; +// DTO.rechr_to_ssr = 0.0; +// DTO.excess = 0.0; +// +//}; +// +//void soil_two_layer::set_K_values() +//{ +// //K_estimate.calculate(); +//}; +// +//void soil_two_layer::set_layer_thaw_fraction() +//{ +// std::cout << "in set layer" << std::endl; +// DTO.thaw_fraction_rechr = 0.0; +// DTO.thaw_fraction_lower = 0.0; +// +// // reference so intention is more clear +// // If we are allowing runoff from infiltration we are also allowing the soil to freeze +// bool& allow_soil_to_freeze = DTO.allow_runoff_from_infiltration; +// double rechr_depth = 0.0; +// double soil_depth = 0.0; +// +// // TODO this might not be right, but the above is what CRHM does +// // depth = storage / porosity, because porosity = storage / depth; +// if (DTO.porosity > 0.0) +// { +// // TODO, porosity in rechr is the same as lower +// rechr_depth = DTO.soil_rechr_max / DTO.porosity / 1000.0; +// soil_depth = DTO.soil_storage_max / DTO.porosity / 1000.0; +// } +// else +// { +// std::cout << "in set layer to exit" << std::endl; +// DTO.thaw_fraction_rechr = 0.0; +// DTO.thaw_fraction_lower = 0.0; +// return; +// }; +// +// std::cout << "in set layer" << std::endl; +// if (allow_soil_to_freeze || DTO.freeze_thaw_first_front == 0.0) +// { +// DTO.thaw_fraction_rechr = 1.0; +// DTO.thaw_fraction_lower = 1.0; +// } +// else +// { +// DTO.thaw_fraction_rechr = 0.0; +// DTO.thaw_fraction_lower = 0.0; +// } +// std::cout << "in set layer" << std::endl; +// +// if (DTO.soil_storage_max > 0.0 && allow_soil_to_freeze && +// DTO.freeze_thaw_first_front > 0.0) +// { +// +// // TODO Verify this calculation +// if (DTO.thaw_front_depth < rechr_depth) +// { +// DTO.thaw_fraction_rechr = DTO.thaw_front_depth / rechr_depth; +// } +// else if ( DTO.thaw_front_depth < soil_depth) +// { +// DTO.thaw_fraction_rechr = 1.0; +// DTO.thaw_fraction_lower = (DTO.thaw_front_depth - rechr_depth) +// / (soil_depth - rechr_depth); +// } +// else +// { +// DTO.thaw_fraction_rechr = 1.0; +// DTO.thaw_fraction_lower = 1.0; +// } +// +// } +// +//}; +// +//void soil_two_layer::set_condensation() +//{ +// if (DTO.actual_ET < 0.0 && DTO.swe == 0.0) +// { +// DTO.condensation = -1.0 * DTO.actual_ET; +// DTO.actual_ET = 0.0; +// } +//}; +// +//void soil_two_layer::organize_soil_layers() +//{ +// if (DTO.soil_storage_max > 0.0) +// { +// double soil_lower_storage = DTO.soil_storage - DTO.soil_rechr_storage; +// +// double potential = DTO.infil + DTO.condensation; +// +// double possible = DTO.thaw_fraction_rechr * (DTO.soil_rechr_max - DTO.soil_rechr_storage); +// if (possible > potential || !DTO.allow_runoff_from_infiltration) +// possible = potential; +// else +// DTO.soil_excess_to_runoff = potential - possible; +// +// DTO.soil_rechr_storage += possible; +// +// if (DTO.soil_rechr_storage > DTO.soil_rechr_max) +// _push_excess_down(DTO.soil_rechr_storage,DTO.soil_rechr_max,soil_lower_storage); +// +// DTO.soil_storage = soil_lower_storage + DTO.soil_rechr_storage; +// +// if (DTO.soil_storage > DTO.soil_storage_max) +// _push_excess_down(DTO.soil_storage,DTO.soil_storage_max,DTO.soil_excess_to_gw); +// +// if (DTO.swe == 0.0) // if there is no snowcover +// { +// DTO.rechr_to_ssr = DTO.soil_rechr_storage / DTO.soil_rechr_max * DTO.K_rechr_to_ssr * DTO.thaw_fraction_rechr; +// DTO.rechr_to_ssr = std::min(DTO.rechr_to_ssr,DTO.soil_rechr_storage * DTO.thaw_fraction_rechr); +// +// DTO.soil_rechr_storage = std::max(0.0, DTO.soil_rechr_storage - DTO.rechr_to_ssr); +// +// DTO.soil_storage -= DTO.rechr_to_ssr; +// DTO.soil_to_ssr = DTO.rechr_to_ssr; +// +// } +// +// if (DTO.soil_excess_to_gw > DTO.K_soil_to_gw * DTO.thaw_fraction_lower) +// { +// double excess_to_gw_max = DTO.K_soil_to_gw * DTO.thaw_fraction_lower; +// _push_excess_down(DTO.soil_excess_to_gw,excess_to_gw_max,DTO.excess); +// } +// +// // Line 607 of SoilX crhmcommetns branch, comment says upper layer but code says lower-layer, ask logan about this. +// if (DTO.excess_to_ssr && DTO.excess > 0.0) +// { +// double excess_to_ssr_max = DTO.excess * (1.0 - DTO.thaw_fraction_lower); +// _push_excess_down(DTO.excess,excess_to_ssr_max,DTO.soil_to_ssr); +// } +// +// +// } +// else +// { +// DTO.excess = DTO.infil + DTO.condensation; +// } +// +// +// +//}; +// +//void soil_two_layer::manage_detention() +//{ +// double face_area = 1.0; // Later will be pulled from face, but since routine_residual is always zero, ignoring it here. +// DTO.soil_excess_to_runoff += DTO.runoff + DTO.excess + DTO.routing_residual / face_area; // routing_residual comes from the crhm varaible redirected_residual which has units of mm*km^2/int (not sure why), so face_area is there for now for consistency. +// +// if (DTO.soil_excess_to_runoff > 0.0) +// { +// if (DTO.swe == 0.0) +// DTO.detention_max = DTO.detention_snow_max; +// else +// DTO.detention_max = DTO.detention_organic_max; +// +// double detention_space = DTO.detention_max - DTO.detention_storage; +// +// if (detention_space > 0.0) +// { +// if (DTO.soil_excess_to_runoff > detention_space) +// { +// DTO.soil_excess_to_runoff = std::max(0.0,DTO.soil_excess_to_runoff - detention_space); +// DTO.detention_storage += detention_space; +// } +// else +// { +// DTO.detention_storage += DTO.soil_excess_to_runoff; +// DTO.soil_excess_to_runoff = 0.0; +// } +// } +// } +// +// if (DTO.detention_storage > 0.0 && DTO.K_detention_to_runoff > 0.0) +// { +// double transfer = std::min(DTO.detention_storage,DTO.K_detention_to_runoff); +// DTO.soil_excess_to_runoff += transfer; +// DTO.detention_storage -= transfer; +// +// if (DTO.detention_storage < 0.0001) // from CRHM, for safety and to drop any Floating-point errors +// DTO.detention_storage = 0.0; +// } +//}; +// +//void soil_two_layer::manage_depression() +//{ +// if (DTO.soil_excess_to_runoff > 0.0 && DTO.depression_max > 0.0) +// { +// double exponent = -1.0 * std::min(12.0,DTO.soil_excess_to_runoff / DTO.depression_max); +// +// double depression_space = (DTO.depression_max - DTO.depression_storage) * (1 - exp(exponent)); +// +// if (DTO.soil_storage_max == 0.0) +// depression_space = DTO.depression_max - DTO.depression_storage; +// +// if (depression_space > 0.0) +// { +// if (DTO.soil_excess_to_runoff > depression_space) +// { +// DTO.soil_excess_to_runoff = std::max(0.0,DTO.soil_excess_to_runoff - depression_space); +// DTO.depression_storage += depression_space; +// // TODO add total tracker +// } +// else +// { +// DTO.depression_storage += DTO.soil_excess_to_runoff; +// // TODO add total tracker +// DTO.soil_excess_to_runoff = 0.0; +// } +// +// } +// } +// +// if (DTO.depression_storage > 0.0 && DTO.K_depression_to_gw > 0.0) +// { +// double value = transfer_min(DTO.depression_storage,DTO.K_depression_to_gw); +// DTO.depression_to_gw += value; +// DTO.depression_storage -= value; +// if (DTO.depression_storage < 0.0) // floatig point error safety? +// DTO.depression_storage = 0.0; +// } +// +// +//}; +// +// +// +//void soil_two_layer::manage_groundwater() +//{ +// DTO.soil_excess_to_gw += DTO.depression_to_gw; +// DTO.depression_to_gw = 0.0; +// DTO.ground_water_storage += DTO.soil_excess_to_gw; +// +// if (DTO.ground_water_storage > DTO.ground_water_max) +// { +// _push_excess_down(DTO.ground_water_storage,DTO.ground_water_max,DTO.ground_water_out); +// } +// +// if (DTO.ground_water_max > 0.0) // divide by zero safety +// { +// double spilled = DTO.ground_water_storage / DTO.ground_water_max * DTO.K_ground_water_out; +// DTO.ground_water_storage -= spilled; +// DTO.ground_water_out += spilled; +// } +// +// // TODO daily and simulation totals, like in CRHM? +// +//}; +// +//void soil_two_layer::manage_subsurface_runoff() +//{ +// if (DTO.depression_storage > 0.0 && DTO.K_depression_to_ssr > 0.0) +// { +// double value = transfer_min(DTO.depression_storage,DTO.K_depression_to_ssr); +// DTO.soil_to_ssr += value; +// DTO.depression_storage -= value; +// if (DTO.depression_storage < 0.0) +// DTO.depression_storage = 0.0; +// } +// +// if (DTO.K_lower_to_ssr > 0.0) +// { +// double available = DTO.soil_storage - DTO.soil_rechr_storage; +// double lower_unfrozen = DTO.K_lower_to_ssr * DTO.thaw_fraction_lower; +// double value = transfer_min(lower_unfrozen,available); +// DTO.soil_storage -= value; +// DTO.soil_to_ssr += value; +// } +// +//}; +// +// +//void soil_two_layer::_push_excess_down(double& layer_storage, double& layer_max, double& layer_down) +//{ +// double excess = layer_storage - layer_max; +// +// layer_storage = layer_max; +// +// layer_down += excess; +//}; +// +//double soil_two_layer::transfer_min(double& val1, double& val2) +//{ +// return std::min(val1,val2); +//}; +// +// +// +// diff --git a/src/modules/soil_submodules/soil_two_layer.hpp b/src/modules/soil_submodules/soil_two_layer.hpp new file mode 100644 index 00000000..381c528b --- /dev/null +++ b/src/modules/soil_submodules/soil_two_layer.hpp @@ -0,0 +1,41 @@ +#pragma once + +#include "soil_base.hpp" +#include "soil_two_layer.hpp" +#include "soil_DTO.hpp" +#include "I_K_estimate.hpp" +#include "soil_classes.hpp" + +class soil_two_layer : public soil_base +{ +public: + soil_two_layer(two_layer_DTO& DTO,I_K_estimate& estimate); + ~soil_two_layer() {}; + + void run() override; + +private: + initializer init; + condensator condensation; + infiltrator infiltrate; + detention_layer detention; + depression_layer depression; + groundwater_layer groundwater; + subsurface_runoff ssr; + I_K_estimate& k_estimator; + + //void initialize_single_step_vars(); + //void set_K_values(); + //void set_layer_thaw_fraction(); // set the new value in a layer, ensure values are between 0 and 1. + //void set_condensation(); + //void organize_soil_layers(); + //void manage_detention(); + //void manage_depression(); + //void manage_groundwater(); + //void manage_subsurface_runoff(); + + + //void _push_excess_down(double& layer_storage, double& layer_max, double& layer_down); + //double transfer_min(double& val1, double& val2); + +}; diff --git a/src/modules/submodule_base.hpp b/src/modules/submodule_base.hpp new file mode 100644 index 00000000..09ec8b12 --- /dev/null +++ b/src/modules/submodule_base.hpp @@ -0,0 +1,16 @@ +#pragma once + +#include +#include +#include "logger.hpp" +#include + +class submodule_base +{ +public: + virtual ~submodule_base() = default; + + virtual void run() = 0; + +}; + diff --git a/src/modules/submodules/base_step.hpp b/src/modules/submodules/base_step.hpp new file mode 100644 index 00000000..6b40175c --- /dev/null +++ b/src/modules/submodules/base_step.hpp @@ -0,0 +1,24 @@ +#pragma once +#include +#include +#include + +#define THROW_NULL_POINTER_EXCEPTION() \ + throw std::runtime_error( \ + std::string("Null pointer at ") + __File__ + ":" + std::to_string(__Line__) \ + ) + + +template +class base_step +{ +public: + explicit base_step(data& _d) : d(_d) {}; + virtual ~base_step() = default; + + virtual void execute() = 0; +protected: + data& d; +}; + + diff --git a/src/modules/submodules/daily_accumulator.hpp b/src/modules/submodules/daily_accumulator.hpp new file mode 100644 index 00000000..ef013928 --- /dev/null +++ b/src/modules/submodules/daily_accumulator.hpp @@ -0,0 +1,58 @@ +#pragma once +#include "base_step.hpp" + +template +concept accumulator_data = requires(T& t) { + // Must have is_new_day() const member function returning bool + { t.is_new_day() } -> std::same_as; + + // Must have steps_per_day() const member function convertible to double + { t.steps_per_day() } -> std::convertible_to; +}; + +template +class daily_accumulator : public base_step +{ + const double* _target_var = nullptr; + double _accumulator = 0.0; + double _mean_value = 0.0; +public: + // overloaded constructors + // normal construction is the first line + // second constructor if you want to initialize the mean_value (starting in the middle of a day). + explicit daily_accumulator(data& _d) : base_step(_d) {}; + explicit daily_accumulator(data& _d, const double value) : base_step(_d), _mean_value(value) {}; + + void bind_to_var(double& var) + { + if (!_target_var) + { + _target_var = &var; + return; + } + throw std::logic_error("_targer_var already bound"); + }; + + const double& get_last_mean() const + { + return _mean_value; + }; + + void execute() override final + { + if (!_target_var) + { + throw std::logic_error("_target_var should be bound before calling execute()"); + return; + }; + + if (this->d.is_new_day()) + { + _mean_value = _accumulator / this->d.steps_per_day(); + _accumulator = 0.0; + }; + + _accumulator += *_target_var; + + }; +}; diff --git a/src/physics/Atmosphere.h b/src/physics/Atmosphere.h index c0035fcd..19a4e4fa 100644 --- a/src/physics/Atmosphere.h +++ b/src/physics/Atmosphere.h @@ -43,6 +43,9 @@ namespace Atmosphere { double corr_precip_slope(double p, double slope); double saturatedVapourPressure(const double& T); + + const double Cp = 1005; // (J/kg/K) volumetric heat capcity of dry air. + const double kappa = 0.4; // proportionality constant, used at least for the PenmanMonteith ET. } diff --git a/src/physics/Soil.cpp b/src/physics/Soil.cpp new file mode 100644 index 00000000..dbd4ee43 --- /dev/null +++ b/src/physics/Soil.cpp @@ -0,0 +1,236 @@ +#include "Soil.h" +#include + +namespace Soil +{ + double _soils_base::lookup(const mymap& map, const std::string key) const + { + auto it = map.find(key); + if (it != map.end()) + { + return it->second; + } + else + { + + CHM_THROW_EXCEPTION(module_error, "Soil Type does not exist in map"); + } + } + + soils_na::soils_na() + { + _make_hash(); + } + + soils_na::~soils_na() + { + + } + + double soils_na::porosity(std::string soil_type) const + { + return lookup(_porosity,soil_type); + } + + double soils_na::pore_size_dist(std::string soil_type) const + { + return lookup(_pore_size_dist,soil_type);//_pore_size_dist[soil_type]; + } + + double soils_na::wilt_point(std::string soil_type) const + { + return lookup(_wilt_point,soil_type); //_wilt_point[soil_type]; + } + + double soils_na::air_entry_tension(std::string soil_type) const + { + return lookup(_air_entry_tension,soil_type); //_air_entry_tension[soil_type]; + } + + double soils_na::capillary_suction(std::string soil_type) const + { + return lookup(_capillary_suction,soil_type);//_capillary_suction[soil_type]; + } + + double soils_na::saturated_conductivity(std::string soil_type) const + { + return lookup(_saturated_conductivity,soil_type); + } + + double soils_na::ayers_texture(std::string texture, std::string ground_cover) const + { + auto it = _ayers_texture.find(texture); + if (it != _ayers_texture.end()) + { + return lookup(it->second, ground_cover); + } + else + CHM_THROW_EXCEPTION(module_error, "Soil Type does not exist in map"); + + } + + void soils_na::check_map() + { + print_ayers_texture(_ayers_texture); + } + + void soils_na::print_ayers_texture(const auto& texture_map) { + for (const auto& [texture, inner_map] : texture_map) { + for (const auto& [ground_cover, value] : inner_map) { + std::cout << texture << " -> " + << ground_cover << " = " + << value << "\n"; + } + } + } + + void soils_na::_make_hash() + { + // illegal to not include this step for the dense_hash_map + _porosity.set_empty_key(""); + _pore_size_dist.set_empty_key(""); + _wilt_point.set_empty_key(""); + _air_entry_tension.set_empty_key(""); + _capillary_suction.set_empty_key(""); + _saturated_conductivity.set_empty_key(""); + _ayers_texture.set_empty_key(""); + + std::string t1 = "sand"; + std::string t2 = "loam"; + std::string t3 = "clay"; + std::string t4 = "loamsand"; + std::string t5 = "sandloam"; + std::string t6 = "siltloam"; + std::string t7 = "saclloam"; + std::string t8 = "clayloam"; + std::string t9 = "siclloam"; + std::string t10 = "sandclay"; + std::string t11 = "siltclay"; + + _porosity[t1] = 0.395; + _porosity[t2] = 0.451; + _porosity[t3] = 0.482; + _porosity[t4] = 0.410; + _porosity[t5] = 0.435; + _porosity[t6] = 0.485; + _porosity[t7] = 0.420; + _porosity[t8] = 0.476; + _porosity[t9] = 0.477; + _porosity[t10] = 0.426; + _porosity[t11] = 0.492; + + _pore_size_dist[t1] = 4.05; + _pore_size_dist[t2] = 5.39; + _pore_size_dist[t3] = 11.4; + _pore_size_dist[t4] = 4.38; + _pore_size_dist[t5] = 4.9; + _pore_size_dist[t6] = 5.3; + _pore_size_dist[t7] = 7.12; + _pore_size_dist[t8] = 8.52; + _pore_size_dist[t9] = 7.75; + _pore_size_dist[t10] = 10.4; + _pore_size_dist[t11] = 11.4; + + _wilt_point[t1] = 0.040; + _wilt_point[t2] = 0.110; + _wilt_point[t3] = 0.215; + _wilt_point[t4] = 0.060; + _wilt_point[t5] = 0.100; + _wilt_point[t6] = 0.130; + _wilt_point[t7] = 0.140; + _wilt_point[t8] = 0.150; + _wilt_point[t9] = 0.190; + _wilt_point[t10] = 0.200; + _wilt_point[t11] = 0.210; + + _air_entry_tension[t1] = 0.121; + _air_entry_tension[t2] = 0.478; + _air_entry_tension[t3] = 0.405; + _air_entry_tension[t4] = 0.090; + _air_entry_tension[t5] = 0.218; + _air_entry_tension[t6] = 0.786; + _air_entry_tension[t7] = 0.299; + _air_entry_tension[t8] = 0.630; + _air_entry_tension[t9] = 0.356; + _air_entry_tension[t10] = 0.153; + _air_entry_tension[t11] = 0.490; + + _capillary_suction[t1] = 0.4950; + _capillary_suction[t2] = 0.8890; + _capillary_suction[t3] = 0.3163; + _capillary_suction[t4] = 0.6130; + _capillary_suction[t5] = 0.1101; + _capillary_suction[t6] = 0.1668; + _capillary_suction[t7] = 0.2185; + _capillary_suction[t8] = 0.2088; + _capillary_suction[t9] = 0.2733; + _capillary_suction[t10] = 0.2390; + _capillary_suction[t11] = 0.2922; + + _saturated_conductivity[t1] = 117.8; + _saturated_conductivity[t2] = 3.4; + _saturated_conductivity[t3] = 0.3; + _saturated_conductivity[t4] = 29.9; + _saturated_conductivity[t5] = 10.9; + _saturated_conductivity[t6] = 6.5; + _saturated_conductivity[t7] = 1.5; + _saturated_conductivity[t8] = 1.0; + _saturated_conductivity[t9] = 1.0; + _saturated_conductivity[t10] = 0.6; + _saturated_conductivity[t11] = 0.5; + + + std::string c1 = "bare_soil"; + std::string c2 = "row_crop"; + std::string c3 = "poor_pasture"; + std::string c4 = "small_grains"; + std::string c5 = "good_pasture"; + std::string c6 = "forested"; + google::dense_hash_map inner_map; + inner_map.set_empty_key(""); + + + // coarse_over_coarse + inner_map[c1] = 7.6; + inner_map[c2] = 12.7; + inner_map[c3] = 15.2; + inner_map[c4] = 17.8; + inner_map[c5] = 25.4; + inner_map[c6] = 76.2; + + _ayers_texture["coarse_over_coarse"] = inner_map; + + // medium over medium + inner_map[c1] = 2.5; + inner_map[c2] = 5.1; + inner_map[c3] = 7.6; + inner_map[c4] = 10.2; + inner_map[c5] = 12.7; + inner_map[c6] = 15.2; + + _ayers_texture["medium_over_medium"] = inner_map; + + // medium/fine over fine + inner_map[c1] = 1.3; + inner_map[c2] = 1.8; + inner_map[c3] = 2.5; + inner_map[c4] = 3.8; + inner_map[c5] = 5.1; + inner_map[c6] = 6.4; + + _ayers_texture["medium_over_fine"] = inner_map; + _ayers_texture["fine_over_fine"] = inner_map; + + + // soil over shallow bedrock + inner_map[c1] = 0.5; + inner_map[c2] = 0.5; + inner_map[c3] = 0.5; + inner_map[c4] = 0.5; + inner_map[c5] = 0.5; + inner_map[c6] = 0.5; + + _ayers_texture["soil_over_bedrock"] = inner_map; + } + +} diff --git a/src/physics/Soil.h b/src/physics/Soil.h index 391b1384..bdee7b91 100644 --- a/src/physics/Soil.h +++ b/src/physics/Soil.h @@ -26,8 +26,75 @@ #pragma once +#include +#include +#include "exception.hpp" + namespace Soil { /********* Soil ************/ -} + enum GATable {PSI, KSAT, WILT, FCAP, PORG, PORE, AIRENT, PORESZ, AVAIL}; // Used for mapping the soil table, PSI and KSAT are used, the others are unused but may but used in the future or other modules. + // see GreenAmpt module in CRHM wiki for details. + // + + using mymap = google::dense_hash_map; + + class _soils_base + { + public: + virtual double porosity(std::string) const = 0; + virtual double pore_size_dist(std::string) const = 0; + virtual double wilt_point(std::string) const = 0; + virtual double air_entry_tension(std::string) const = 0; + virtual double capillary_suction(std::string) const = 0; + virtual double saturated_conductivity(std::string) const = 0; + + virtual double ayers_texture(std::string texture, std::string ground_cover) const = 0; + double lookup(const mymap& map, const std::string key) const; + virtual ~_soils_base() = default; + + private: + virtual void _make_hash() = 0; + }; + + class soils_na : public _soils_base + { + public: + soils_na(); + ~soils_na(); + + double porosity(std::string soil_type) const override; + double pore_size_dist(std::string soil_type) const override; + double wilt_point(std::string soil_type) const override; + double air_entry_tension(std::string soil_type) const override; + double capillary_suction(std::string soil_type) const override; + double saturated_conductivity(std::string soil_type) const override; + double ayers_texture(std::string texture, std::string ground_cover) const override; + + void check_map(); + private: + + mymap _porosity; + mymap _pore_size_dist; + mymap _wilt_point; + mymap _air_entry_tension; + mymap _capillary_suction; + mymap _saturated_conductivity; + + google::dense_hash_map< std::string, mymap> _ayers_texture; + + void _make_hash() override; + + void print_ayers_texture(const auto& texture_map); + + + }; + + template + const T& get_soil_obj() + { + static const T soil_obj; + return soil_obj; + }; +}; diff --git a/src/tests/googletest b/src/tests/googletest index b155875f..57193061 160000 --- a/src/tests/googletest +++ b/src/tests/googletest @@ -1 +1 @@ -Subproject commit b155875f32dc74e293d96c0de2dfcdfa913804e4 +Subproject commit 571930618fa96eabcd05b573285edbee9fc13bae diff --git a/src/tests/main.cpp b/src/tests/main.cpp index 43ad65d5..49c7d6e9 100644 --- a/src/tests/main.cpp +++ b/src/tests/main.cpp @@ -11,16 +11,16 @@ int main(int argc, char* argv[]) ::testing::InitGoogleTest(&argc, argv); -#ifdef USE_MPI - boost::mpi::environment _mpi_env; - boost::mpi::communicator _comm_world; - - ::testing::TestEventListeners& listeners = - ::testing::UnitTest::GetInstance()->listeners(); - if (_comm_world.rank() != 0) { - delete listeners.Release(listeners.default_result_printer()); - } -#endif +//#ifdef USE_MPI +// boost::mpi::environment _mpi_env; +// boost::mpi::communicator _comm_world; +// +// ::testing::TestEventListeners& listeners = +// ::testing::UnitTest::GetInstance()->listeners(); +// if (_comm_world.rank() != 0) { +// delete listeners.Release(listeners.default_result_printer()); +// } +//#endif // MPI_Init(&argc, &argv); result = RUN_ALL_TESTS(); diff --git a/src/tests/submoduletests/CSVreader.hpp b/src/tests/submoduletests/CSVreader.hpp new file mode 100644 index 00000000..87e808fe --- /dev/null +++ b/src/tests/submoduletests/CSVreader.hpp @@ -0,0 +1,163 @@ +#include +#include +#include +#include +#include +#include +#include + +class CSVReader { +private: + std::string filePath; + std::vector headers; + std::vector> data; + bool isLoaded = false; + + + + // Helper function to convert string to type T + template + T convert(const std::string& value) const { + std::istringstream iss(value); + T result; + if (!(iss >> result)) { + throw std::runtime_error("Failed to convert value: " + value); + } + return result; + } + + // Specialization for std::string (no conversion needed) + template<> + std::string convert(const std::string& value) const { + return value; + } + + void loadData() + { + std::ifstream file(filePath); + if (!file.is_open()) { + throw std::runtime_error("Error opening file at: " + filePath); + } + + std::string line; + // Read headers + if (std::getline(file, line)) { + std::stringstream ss(line); + std::string value; + while (std::getline(ss, value, ',')) { + headers.push_back(value); + } + } + + // Read data rows + while (std::getline(file, line)) { + std::vector row; + std::stringstream ss(line); + std::string value; + while (std::getline(ss, value, ',')) { + row.push_back(value); + } + data.push_back(row); + } + + isLoaded = true; + } + +public: + + // Constructor with hardcoded Mac path + CSVReader() : filePath("/Users/hin601/Documents/TestBuild/CRHM_Compare/PythonScripts/unittestscripts/Cleaned_data.csv") { + if (!std::filesystem::exists(filePath)) { + throw std::runtime_error("File not found: " + filePath); + } + + loadData(); + } + + + + // Function to get value from CSV + template + T getValue(const std::string& columnName, int rowNumber) + { + if (!isLoaded) + { + throw std::runtime_error("Data not loaded"); + } + + // Find the target column + auto it = std::find(headers.begin(), headers.end(), columnName); + if (it == headers.end()) + { + throw std::runtime_error("Column '" + columnName + "' not found."); + } + int targetColumn = std::distance(headers.begin(), it); + + // Check row bounds + if (rowNumber < 0 || rowNumber >= static_cast(data.size())) + { + throw std::runtime_error("Row " + std::to_string(rowNumber) + " is out of bounds"); + } + + // Check column bounds + if (targetColumn >= static_cast(data[rowNumber].size())) + { + throw std::runtime_error("Column index out of bounds for row " + + std::to_string(rowNumber)); + } + + return convert(data[rowNumber][targetColumn]); + } + + //std::ifstream file(filePath); + //std::string line, value; + //std::vector columnNames; + //int targetColumn = -1; + //int currentRow = 0; + + //// Check if file opened successfully + //if (!file.is_open()) { + // throw std::runtime_error("Error opening file at: " + filePath); + //} + + //// Read header row to get column names + //if (std::getline(file, line)) { + // std::stringstream ss(line); + // while (std::getline(ss, value, ',')) { + // columnNames.push_back(value); + // } + + // // Find the target column + // for (size_t i = 0; i < columnNames.size(); ++i) { + // if (columnNames[i] == columnName) { + // targetColumn = static_cast(i); + // break; + // } + // } + + // if (targetColumn == -1) { + // throw std::runtime_error("Column '" + columnName + "' not found."); + // } + //} + + //// Read data rows until we reach the target row + //while (std::getline(file, line)) { + // if (currentRow == rowNumber) { + // std::vector rowData; + // std::stringstream ss(line); + // while (std::getline(ss, value, ',')) { + // rowData.push_back(value); + // } + + // if (targetColumn < rowData.size()) { + // return convert(rowData[targetColumn]); + // } else { + // throw std::runtime_error("Column index out of bounds for row " + std::to_string(rowNumber)); + // } + // } + // currentRow++; + //} + // + //throw std::runtime_error("Row " + std::to_string(rowNumber) + " not found in the CSV file."); + //} +}; diff --git a/src/tests/submoduletests/test_CSVreader.cpp b/src/tests/submoduletests/test_CSVreader.cpp new file mode 100644 index 00000000..df6f7303 --- /dev/null +++ b/src/tests/submoduletests/test_CSVreader.cpp @@ -0,0 +1,20 @@ +#include "gtest/gtest.h" +#include "CSVreader.hpp" + +TEST(CSVreaderTeat, CheckFunction) +{ + CSVReader reader; + double value = reader.getValue("net_rain",1); + + ASSERT_DOUBLE_EQ(value,0.0); + +}; + +TEST(CSVreaderTeat, CheckLaterValue) +{ + CSVReader reader; + double value = reader.getValue("net_rain",136768); + + ASSERT_DOUBLE_EQ(value,0.213297); + +}; diff --git a/src/tests/submoduletests/test_XG.cpp b/src/tests/submoduletests/test_XG.cpp new file mode 100644 index 00000000..b82e7078 --- /dev/null +++ b/src/tests/submoduletests/test_XG.cpp @@ -0,0 +1,815 @@ +#include "XG_algorithm.hpp" +#include +#include +#include "gtest/gtest.h" +/* + * XGStateTest: Wrapper class for tests + * XGStateTest is effectively a mock of Infil_All module but done indirectly. Due to the complexity of the module classes, it was easier to write this. + * The member variables with the _ prefix are inputs that are supplied to the constructor of Crack. + * Default values are given and used for most tests. + * Other member variables are parameters that are also supplied to Crack unless it has the const specifier, then it is just useful for these tests. + * Member functions are just tools to enable the tests. + * Initialization of XGStateTest assumes that the frozen period has just begun. + * + */ + +class XGStateTest : public testing::Test +{ +protected: + XGStateTest() + { + init_vectors(); + }; + + double surface_temperature = 0.0; + std::unique_ptr set_default_P(); + std::unique_ptr set_default_S(const int& N); + std::unique_ptr P; + std::unique_ptr S; + void distribute_moisture(); + + // params + int num_layers = 4; + std::vector depths; + std::vector porosity; + std::vector theta_default; + std::vector dry_soil_k; + std::vector saturated_frozen_soil_k; + std::vector saturated_unfrozen_soil_k; + double Trigthrhld = 100.0; + double theta_min = 0.001; + double SWE_k = 0.35; + double Zpf_init = 2.0; + bool freeze_kw_ki_update = true; + bool thaw_ki_kw_update = true; + int k_update = 1; + double soil_rechr_max = 350.0; + double soil_moist_max = 625.0; + int time_step_per_day = 24; + bool calc_conductivity = false; + + void init_vectors() + { + depths.resize(num_layers,0.5); + porosity.resize(num_layers,0.5); + theta_default.resize(num_layers,0.5); + dry_soil_k.resize(num_layers,2.5); + saturated_frozen_soil_k.resize(num_layers,1.98); + saturated_frozen_soil_k[0] = 1.55; + saturated_unfrozen_soil_k.resize(num_layers,1.67); + saturated_unfrozen_soil_k[0] = 0.8; + }; +}; + +// Use as first XGTest +//TEST_F(XGStateTest, NoRunDefaultOutput) { +// +// XG_algorithm xg(surface_temperature,S,P); +// +// double thaw = xg.get_thaw_front_depth(); +// +// double freeze = xg.get_freeze_front_depth(); +// +// double first = xg.get_first_front_depth(); +// +// ASSERT_EQ(thaw,0.0); +// ASSERT_EQ(freeze,0.0); +// ASSERT_EQ(first,0.0); +//}; +// + +std::unique_ptr XGStateTest::set_default_S(const int& N) +{ + std::unique_ptr S = + std::make_unique(N,*P); + //S.set_layer_moisture_maximums(P) + // .set_thermal_conductivities(P) + // .set_freezethaw_ratios(P) + // .set_initial_freezethaw_depths(Zdf,Zdt); + + + //S.nfront = 0; + //S.Fz_low = 1; + //S.Th_low = 1; + //S.t_trend = 0.0; + //S.Zd_front.assign(P.n,0.0); + //S.XG_moist_d = 0.0; + //S.XG_rechr_d = 0.0; + //S.tc_composite.assign(P.n,0.0); + //S.tc_composite2.assign(P.n,0.0); + //S.XG_max.assign(P.n); + // + //std::transform(P.por.begin(),P.por.end(), + // P.depths.begin(), + // S.XG_max.begin(), + // [](double x,double y) {return x* y;} + // ); + // + //S.theta = P.theta_default; + // + return S; + +}; + +std::unique_ptr XGStateTest::set_default_P(void) +{ + std::unique_ptr P = + std::make_unique(depths, + Trigthrhld, + porosity, + num_layers, + theta_default, + theta_min, + dry_soil_k, + saturated_frozen_soil_k, + saturated_unfrozen_soil_k, + SWE_k, + Zpf_init, + freeze_kw_ki_update, + thaw_ki_kw_update, + k_update, + soil_rechr_max, + soil_moist_max, + time_step_per_day, + calc_conductivity); + + return P; +}; + +TEST_F(XGStateTest,TestParams) +{ + P = set_default_P(); + + EXPECT_EQ(P->por.size(),num_layers); + EXPECT_EQ(P->getdepths().size(),num_layers); + EXPECT_EQ(P->theta_default.size(),num_layers); +} + +TEST_F(XGStateTest,MoistureMaxTest) +{ + P = set_default_P(); + S = set_default_S(P->N_Soil_layers); + S->set_layer_moisture_maximums(*P); + for (int i=0; i < P->N_Soil_layers; ++i) + { + EXPECT_EQ(S->XG_max[i],0.5*0.5*1000.0); + EXPECT_EQ(S->theta[i],theta_default[i]); + + if (i == 0) + { + EXPECT_EQ(S->rechr_fract[i],1.0); + EXPECT_EQ(S->moist_fract[i],0.0); + EXPECT_EQ(S->default_fract[i],0.0); + } + else if (i == 1) + { + EXPECT_EQ(S->rechr_fract[i], + (soil_rechr_max - S->XG_max[0])/S->XG_max[1]); + EXPECT_EQ(S->moist_fract[i], + 1.0 - S->rechr_fract[i]); + EXPECT_EQ(S->default_fract[i],0.0); + } + else if (i == 2) + { + EXPECT_EQ(S->rechr_fract[i],0.0); + EXPECT_EQ(S->moist_fract[i],1.0); + EXPECT_EQ(S->default_fract[i],0.0); + } + else if (i == 3) + { + EXPECT_EQ(S->rechr_fract[i],0.0); + EXPECT_EQ(S->moist_fract[i], + (soil_moist_max - S->XG_max[1]*(1 - S->rechr_fract[1]) - S->XG_max[2])/S->XG_max[3]); + EXPECT_EQ(S->default_fract[i], + 1 - S->moist_fract[i]); + } + else if (i == 4) + { + EXPECT_EQ(S->rechr_fract[i],0.0); + EXPECT_EQ(S->moist_fract[i],0.0); + EXPECT_EQ(S->default_fract[i],1.0); + } + + + } + + EXPECT_EQ(S->XG_rechr_d, + 0 + depths[0] + depths[1]*(soil_rechr_max - S->XG_max[0])/S->XG_max[1]); + EXPECT_EQ(S->XG_moist_d, + 0 + depths[1] + depths[2] + depths[3]*(soil_moist_max - S->XG_max[1]*(1 - (soil_rechr_max - S->XG_max[0])/S->XG_max[1]) - S->XG_max[2])/S->XG_max[3]); + + + +}; + +TEST_F(XGStateTest,MoistureMaxTest_Case2) +{ + num_layers = 1; + soil_rechr_max = 100; + soil_moist_max = 250; + P = set_default_P(); + S = set_default_S(P->N_Soil_layers); + S->set_layer_moisture_maximums(*P); + + EXPECT_EQ(S->XG_rechr_d, + P->getdepths()[0]*soil_rechr_max/S->XG_max[0]); + EXPECT_EQ(S->XG_moist_d, + P->getdepths()[0]*soil_moist_max/S->XG_max[0]); + EXPECT_EQ(S->rechr_fract[0],soil_rechr_max/S->XG_max[0]); + EXPECT_EQ(S->moist_fract[0],(soil_moist_max - soil_rechr_max)/S->XG_max[0]); + EXPECT_EQ(S->default_fract[0],1 - soil_moist_max/S->XG_max[0]); +}; + +TEST_F(XGStateTest,GetTtcAndGetFtcTest) +{ + //calc_conductivity off + P = set_default_P(); + S = set_default_S(P->N_Soil_layers); + int i = 2; + S->layer_h2o[i] = 23.0; + + double k_f = S->get_ftc(i,*P); + + double k_t = S->get_ttc(i,*P); + + double k_f_val = (1.0 - P->por[i])*P->soil_solid_km[i] + + S->layer_h2o[i]/1000.0*XG_algorithm::kw + (P->por[i] - S->layer_h2o[i]/1000.0)*XG_algorithm::ka; + double k_t_val = (1.0 - P->por[i])*P->soil_solid_km[i] + + S->layer_h2o[i]/1000.0*XG_algorithm::ki + (P->por[i] - S->layer_h2o[i]/1000.0)*XG_algorithm::ka; + + EXPECT_EQ(k_f,k_f_val) << "off"; + EXPECT_EQ(k_t,k_t_val) << "off"; + + // calc_conductivity on + calc_conductivity = true; + P = set_default_P(); + + k_f = S->get_ftc(i,*P); + k_t = S->get_ttc(i,*P); + + k_f_val = (P->soil_solid_km_kw[i] - P->soil_solid_km[i])*std::pow(S->layer_h2o[i]/(1000.0*P->por[i]),2) + P->soil_solid_km[i]; + k_t_val = P->soil_solid_km[i]*std::pow(P->soil_solid_km_ki[i]/P->soil_solid_km[i],S->layer_h2o[i]/(1000.0*P->por[i])); + + EXPECT_EQ(k_f,k_f_val) << "on"; + EXPECT_EQ(k_t,k_t_val) << "on"; + +}; + +TEST_F(XGStateTest,SetThermalConductivityTest) +{ + num_layers = 2; + P = set_default_P(); + S = set_default_S(P->N_Soil_layers); + S->ftc_contents[0] = 1; S->ttc_contents[0] = 1; + S->ftc_contents[1] = 0; S->ttc_contents[1] = 0; + const double soil_rechr = 100.0; + const double soil_moist = 300.0; + S->set_layer_moisture_maximums(*P) + .set_thermal_conductivities(*P,soil_moist,soil_rechr); + + double XG_moist_compare = (S->rechr_fract[0]*soil_rechr/soil_rechr_max + S->moist_fract[0]*(soil_moist - soil_rechr)/(soil_moist_max - soil_rechr_max) + S->default_fract[0]*P->theta_default[0])*S->XG_max[0]; + EXPECT_DOUBLE_EQ(S->XG_moist[0],XG_moist_compare) << "Loop 0, XG moist"; + + double check_XG_moist_compare = (S->rechr_fract[0]*soil_rechr/soil_rechr_max + S->moist_fract[0]*(soil_moist - soil_rechr)/(soil_moist_max - soil_rechr_max)) * S->XG_max[0]; + + EXPECT_DOUBLE_EQ(S->theta[0], + S->XG_moist[0] / S->XG_max[0]) << "Loop 0, theta"; + + EXPECT_DOUBLE_EQ(S->layer_h2o[0],1000*S->theta[0]*P->por[0]) << "Loop 0, h2o"; + + EXPECT_DOUBLE_EQ(S->ftc[0],S->get_ttc(0,*P)) << "Loop " << 0 << " ftc"; + + EXPECT_DOUBLE_EQ(S->ttc[0],S->get_ftc(0,*P)) << "Loop 0, ttc"; + + EXPECT_DOUBLE_EQ(S->pf[0],0.0); + EXPECT_DOUBLE_EQ(S->pt[0],0.0); + + XG_moist_compare = (S->rechr_fract[1]*soil_rechr/soil_rechr_max + S->moist_fract[1]*(soil_moist - soil_rechr)/(soil_moist_max - soil_rechr_max) + S->default_fract[1]*P->theta_default[1])*S->XG_max[1]; + EXPECT_DOUBLE_EQ(S->XG_moist[1],XG_moist_compare) << "Loop 1, XG moist"; + + check_XG_moist_compare += (S->rechr_fract[1]*soil_rechr/soil_rechr_max + S->moist_fract[1]*(soil_moist - soil_rechr)/(soil_moist_max - soil_rechr_max))*S->XG_max[1]; + + EXPECT_DOUBLE_EQ(S->check_XG_moist,check_XG_moist_compare) << "check_XG_moist"; + + EXPECT_DOUBLE_EQ(S->theta[1], + S->XG_moist[1] / S->XG_max[1]) << "Loop 1, theta"; + + EXPECT_DOUBLE_EQ(S->layer_h2o[1],1000*S->theta[1]*P->por[1]) << "Loop 1, h2o"; + + EXPECT_DOUBLE_EQ(S->ftc[1],S->get_ftc(1,*P)) << "Loop " << 1 << " ftc"; + + EXPECT_DOUBLE_EQ(S->ttc[1],S->get_ttc(1,*P)) << "Loop 1, ttc"; + + + +}; + +TEST_F(XGStateTest,SetThermalConductivityTest2) +{ + k_update = 0; + num_layers = 1; + P = set_default_P(); + S = set_default_S(P->N_Soil_layers); + S->ftc_contents[0] = 1; S->ttc_contents[0] = 1; + const double soil_rechr = 100.0; + const double soil_moist = 300.0; + S->set_layer_moisture_maximums(*P) + .set_thermal_conductivities(*P,soil_moist,soil_rechr); + + double XG_moist_compare = (S->rechr_fract[0]*soil_rechr/soil_rechr_max + S->moist_fract[0]*(soil_moist - soil_rechr)/(soil_moist_max - soil_rechr_max) + S->default_fract[0]*P->theta_default[0])*S->XG_max[0]; + EXPECT_EQ(S->XG_moist[0],XG_moist_compare) << "Loop 0, XG moist"; + + double check_XG_moist_compare = (S->rechr_fract[0]*soil_rechr/soil_rechr_max + S->moist_fract[0]*(soil_moist - soil_rechr)/(soil_moist_max - soil_rechr_max)) * S->XG_max[0]; + EXPECT_EQ(S->check_XG_moist,check_XG_moist_compare) << "Loop 0, check_XG_moist"; + + EXPECT_EQ(S->theta[0], + S->XG_moist[0] / S->XG_max[0]) << "Loop 0, theta"; + + EXPECT_EQ(S->layer_h2o[0],1000*S->theta[0]*P->por[0]) << "Loop 0, h2o"; + + EXPECT_EQ(S->ftc[0],S->get_ftc(0,*P)) << "Loop " << 0 << " ftc"; + + EXPECT_EQ(S->ttc[0],S->get_ttc(0,*P)) << "Loop 0, ttc"; +}; + +TEST_F(XGStateTest,SetThermaConductivityTest3) +{ + soil_moist_max = 0.0; + num_layers = 1; + P = set_default_P(); + S = set_default_S(P->N_Soil_layers); + S->ftc_contents[0] = 1; S->ttc_contents[0] = 1; + const double soil_rechr = 100.0; + const double soil_moist = 300.0; + S->set_layer_moisture_maximums(*P) + .set_thermal_conductivities(*P,soil_moist,soil_rechr); + + EXPECT_EQ(S->XG_moist[0],S->default_fract[0]*P->theta_default[0]*S->XG_max[0]); + +}; + +TEST_F(XGStateTest,SetFreezeThawRatiosTest) +{ + num_layers = 7; + init_vectors(); + P = set_default_P(); + S = set_default_S(P->N_Soil_layers); + for (int i = 0; i < P->N_Soil_layers; ++i) + { + S->ftc_contents[i] = 1; S->ttc_contents[i] = 1; + } + const double soil_rechr = 100.0; + const double soil_moist = 300.0; + + S->set_layer_moisture_maximums(*P) + .set_thermal_conductivities(*P,soil_moist,soil_rechr) + .set_freezethaw_ratios(*P); + + for (int i = 0; i < P->N_Soil_layers; ++i) + { + if (i == 0) + { + EXPECT_DOUBLE_EQ(S->pf[i], 0.0); + EXPECT_DOUBLE_EQ(S->pt[i], 0.0); + } + else + { + EXPECT_DOUBLE_EQ(S->pf[i],std::sqrt(S->ftc[i-1]*S->layer_h2o[i]/(S->ftc[i]*S->layer_h2o[i-1]))) + << "Loop " << i; + EXPECT_DOUBLE_EQ(S->pt[i],std::sqrt(S->ttc[i-1]*S->layer_h2o[i]/(S->ttc[i]*S->layer_h2o[i-1]))) + << "Loop " << i; + } + } + +}; + +TEST_F(XGStateTest,PushFrontFunction) +{ + P = set_default_P(); + S = set_default_S(P->N_Soil_layers); + + S->nfront=0; + double input = 15.222; + S->push_front(input); + + std::string message = "nfront 0"; + EXPECT_EQ(S->nfront,1) << message; + EXPECT_EQ(S->Zd_front[2],input) << message; + + S->nfront=0; + std::vector Zd{10.0,3.0}; + S->Zd_front[0] = Zd[0]; + S->Zd_front[1] = Zd[1]; + S->push_front(input); + + message = "nonzero ZdFront, nfront 0"; + EXPECT_EQ(S->nfront,1) << message; + EXPECT_EQ(S->Zd_front[0], Zd[0]) << message; + EXPECT_EQ(S->Zd_front[1], Zd[1]) << message; + EXPECT_EQ(S->Zd_front[2],input) << message; + for (int i = 3; i < S->Zd_front.size(); ++i) + EXPECT_EQ(S->Zd_front[i],0.0) << message; + + S->nfront = 2; + Zd.push_back(1.3); + Zd.push_back(4.9); + + for (int i = 0; i < S->Zd_front.size(); ++i) + S->Zd_front[i] = Zd[i]; + + S->push_front(input); + + message = "Nonzero Zd_front, nfront 2"; + for (int i = 0; i < S->Zd_front.size(); ++i) + { + if (i == 2) + EXPECT_EQ(S->Zd_front[i],input) << message << " Loop: " << i; + else if (i == 3 || i == 4) + EXPECT_EQ(S->Zd_front[i],Zd[i-1]) << message << " Loop: " << i; + else + EXPECT_EQ(S->Zd_front[i],Zd[i]) << message << " Loop: " << i; + } + EXPECT_EQ(S->nfront,3) << message; +}; + +TEST_F(XGStateTest,AccumulateDegreeDays) +{ + P = set_default_P(); + S = set_default_S(P->N_Soil_layers); + + std::vector t{13.5,2.3,-1.0,3.9,-4.6}; + double my_B = 0.0, my_TrigAcc = 0.0, myt_trend = 0.0; + + for (int i = 0; i < t.size(); ++i) + { + S->accumulate_degree_days(t[i]); + my_B += t[i]/24; + my_TrigAcc += my_B; + myt_trend -= myt_trend/192; + myt_trend += my_B/192; + EXPECT_EQ(S->B,my_B); + EXPECT_EQ(S->t_trend,myt_trend); + EXPECT_EQ(S->TrigAcc,my_TrigAcc); + } +}; + +TEST_F(XGStateTest,DetermineFreezeThawIdle) +{ + num_layers = 6; + init_vectors(); + P = set_default_P(); + S = set_default_S(P->N_Soil_layers); + + S->TrigAcc = P->Trigthrhld+100.0; + S->TrigState = 0; + + S->determine_freeze_thaw_idle(); + + EXPECT_EQ(S->TrigAcc,P->Trigthrhld); + + S->TrigAcc = -P->Trigthrhld - 100.0; + S->t_trend = 999.9; + + S->determine_freeze_thaw_idle(); + + EXPECT_EQ(S->TrigAcc,-P->Trigthrhld); + EXPECT_EQ(S->TrigState,-1); + EXPECT_EQ(S->t_trend,0.0); + + S->TrigAcc = P->Trigthrhld; + S->t_trend = 100.0; + double Zdt = 1.2; + double Zdf = 0.75; + S->Zdf = Zdf; + S->Zdt = Zdt; + std::vector Zd_init{-Zdf,Zdt,1.4,-1.5,0.0,0.0}; + S->Zd_front = Zd_init; + S->nfront = 2; + S->determine_freeze_thaw_idle(); + + EXPECT_EQ(S->TrigState,0); + + EXPECT_EQ(S->Zd_front[0],0.0); + EXPECT_EQ(S->Zd_front[1],-Zdf); + EXPECT_EQ(S->Zd_front[2],Zdt); + EXPECT_EQ(S->Zd_front[3],Zd_init[2]); + EXPECT_EQ(S->Zd_front[4],Zd_init[3]); + EXPECT_EQ(S->Zdt,0.0); + + S->TrigAcc = -P->Trigthrhld; + S->TrigState = 1; + S->t_trend = -100.0; + + Zdt = 0.75; + Zdf = 1.2; + S->Zdf = Zdf; + S->Zdt = Zdt; + Zd_init = {Zdt,-Zdf,-1.5,1.6,-1.8,0.0}; + S->Zd_front = Zd_init; + S->nfront = 3; + S->determine_freeze_thaw_idle(); + + EXPECT_EQ(S->TrigState,0); + + EXPECT_EQ(S->Zd_front[0],0.0); + EXPECT_EQ(S->Zd_front[1],Zdt); + EXPECT_EQ(S->Zd_front[2],-Zdf); + EXPECT_EQ(S->Zd_front[3],Zd_init[2]); + EXPECT_EQ(S->Zd_front[4],Zd_init[3]); + EXPECT_EQ(S->Zdf,0.0); +}; + + + + + + + +// TODO Mock state and param and move this to another file +class XGTest : public XGStateTest +{ +protected: + XGTest() {}; + const double soil_moist = 300; + const double soil_rechr = 100; + + void do_my_set_up() + { + P = set_default_P(); + S = set_default_S(P->N_Soil_layers); + + S->set_layer_moisture_maximums(*P) + .set_thermal_conductivities(*P,soil_moist,soil_rechr) + .set_freezethaw_ratios(*P); + } + + struct CRHM + { + double TrigAcc; + double TrigState; + double t_trend; + double Zdf; + double Zdt; + double B; + double Bth; + double Bfr; + double hru_tsf; + + CRHM(const int& i,CSVReader& reader) + { + TrigAcc = reader.getValue("TrigAcc",i); + TrigState = reader.getValue("TrigState",i); + t_trend = reader.getValue("t_trend",i); + Zdf = reader.getValue("Zdf",i); + Zdt = reader.getValue("Zdt",i); + B = reader.getValue("B",i); + Bth = reader.getValue("Bth",i); + Bfr = reader.getValue("Bfr",i); + hru_tsf = reader.getValue("hru_tsf",i); + }; + }; +}; + + + +TEST_F(XGTest,BasicOvertakenThawFrontTest) +{ + std::vector TrigAcc{-Trigthrhld,Trigthrhld, + Trigthrhld,-Trigthrhld,-Trigthrhld,-Trigthrhld}; + std::vector t_trend(TrigAcc.size(),1.0); + t_trend[0] *= -1.0; + t_trend[3] *= -1.0; + t_trend[4] *= -1.0; + t_trend[5] *= -1.0; + + std::vector Zd_front; + + std::vector t_surface{-20.0,15.0,15.0,-10.0,-10.0,-8.5}; // Set this, essential for computing z_f and z_t, use it to control what kind of phase we are in (e.g., do fronts merge after a frozen/thaw front is overtaken by the surface thaw/freeze front. + + do_my_set_up(); + + for (int ii = 0; ii < TrigAcc.size(); ++ii) + { + + S->is_newday = true; + S->TrigAcc = TrigAcc.at(ii); + S->t_trend = t_trend.at(ii); + XG_algorithm XG(t_surface.at(ii),soil_moist,soil_rechr,*S,*P); + XG.run(); + if (ii == 0 || ii == 1) + { + // Freeze at 0, Idle at 1 + // + // Array at this stage: + // [Z[0]] < 0 + // [0] + // [0] + // [0] + // [0] + + EXPECT_TRUE(S->Zd_front.at(0) < 0.0) << "Loop: " << ii; + for (int jj = 1; jj < S->Zd_front.size(); ++jj) + { + EXPECT_TRUE(S->Zd_front.at(jj) == 0.0) << "Loop: " << ii; + } + if (ii == 0) + Zd_front.push_back(S->Zd_front.at(0)); + } + else if (ii == 2) + { + // Thaw at 2 + // + // Array at this stage: + // [Z[1]] > 0 + // [Z[0]] < 0 + // [0] + // [0] + // [0] + + EXPECT_TRUE(S->Zd_front.at(0) > 0.0) << "Loop: " << ii; + EXPECT_EQ(S->Zd_front.at(1),Zd_front.at(0)) << "Loop: " << ii; + for (int jj = 2; jj < S->Zd_front.size(); ++jj) + { + EXPECT_TRUE(S->Zd_front.at(jj) == 0.0) << "Loop: " << ii; + } + Zd_front.push_back(S->Zd_front.at(0)); + } + else if (ii == 3) + { + // Idle at 3 + // + // Array at this stage: + // [0] + // [Z[1]] > 0 + // [Z[0]] < 0 + // [0] + // [0] + + EXPECT_EQ(S->Zd_front.at(0),0.0) << "Loop: " << ii; + EXPECT_EQ(S->Zd_front.at(1),Zd_front.at(1)) << "Loop: " << ii; + EXPECT_EQ(S->Zd_front.at(2),Zd_front.at(0)) << "Loop: " << ii; + for (int jj = 3; jj < S->Zd_front.size(); ++jj) + { + EXPECT_TRUE(S->Zd_front.at(jj) == 0.0) << "Loop: " << ii; + } + } + else if (ii == 4) + { + // Freeze at 4 + // + // Array at this stage: + // [Z[2]] < 0 + // [Z[1]] > 0 + // [Z[0]] < 0 + // [0] + // [0] + EXPECT_TRUE(S->Zdf < S->Zdt) << "Loop: " << ii; + EXPECT_TRUE(S->Zd_front.at(0) < 0.0) << "Loop: " << ii; + EXPECT_EQ(S->Zd_front.at(1),Zd_front.at(1)) << "Loop: " << ii; + EXPECT_EQ(S->Zd_front.at(2),Zd_front.at(0)) << "Loop: " << ii; + for (int jj = 3; jj < S->Zd_front.size(); ++jj) + { + EXPECT_TRUE(S->Zd_front.at(jj) == 0.0) << "Loop: " << ii; + } + Zd_front.push_back(S->Zd_front.at(0)); + } + else if (ii == 5) + { + // Freeze at 5 + // + // Array at this stage: + // [Z[3]] < 0 + // [0] + // [0] + // [0] + // [0] + EXPECT_EQ(S->Zdf,-S->Zd_front[0]) << "Loop: " << ii; + EXPECT_TRUE(S->Zd_front[0] < 0) << "Loop: " << ii; + EXPECT_EQ(S->Zdt,0.0) << "Loop: " << ii; + + for (int jj = 1; jj < S->Zd_front.size(); ++jj) + EXPECT_TRUE(S->Zd_front.at(jj) == 0.0) << "Loop: " << ii; + } + } + //EXPECT_EQ(S->Zd_front[0],z_1); + //for (int ii = 1; ii < S->Zd_front.size(); ++ii) + //{ + // EXPECT_TRUE(S->Zd_front[ii] == 0.0); + //} + + //S->TrigAcc = P->Trigthrhld; + //S->t_trend = 1.0; + + //EXPECT_TRUE(S->Zd_front[0] > 0.0); + //std::cout << S->TrigState <Zd_front[1],z_1); + //for (int ii = 1; ii < S->Zd_front.size(); ++ii) + //{ + // EXPECT_TRUE(S->Zd_front[ii] == 0.0); + //} + + //double z_2 = S->Zd_front[0]; + +}; + +#define diff 0.000001 +TEST_F(XGTest,FullZdFrontOrganizeTest) +{ + num_layers=6; + do_my_set_up(); + std::vector Z{-0.1,0.234,-0.758,0.9,-1.3,0.}; + S->Zd_front = Z; + S->nfront=3; + S->TrigState = -1; + S->t_trend = -1; + S->Zdt = S->Zd_front[1]; + S->Zdf = -S->Zd_front[0]; + double t_surface = -300.0; + S->TrigAcc = -Trigthrhld; + S->is_newday = true; + + XG_algorithm XG(t_surface,soil_moist,soil_rechr,*S,*P); + XG.run(); + + for (int i = 0; i < Z.size()-2; ++i) + { + EXPECT_NEAR(S->Zd_front.at(i),Z.at(i+2),diff) << "Loop: " << i; + } + for (int i = Z.size()-2; i < Z.size(); ++i) + EXPECT_EQ(S->Zd_front.at(i),0.0) << "Loop: " << i; + + EXPECT_NEAR(S->Zdf,std::fabs(Z[2]),diff); + EXPECT_NEAR(S->Zdt,Z[3],diff); + +}; + +#define diff4 0.0001 +#define diff3 0.001 +#define diff5 0.00001 +//TEST_F(XGTest,LongTimeTest) +//{ +////#ifdef NDEBUG +//// std::cout << "NDEBUG is defined (asserts are disabled)\n"; +////#else +//// std::cout << "NDEBUG is NOT defined (asserts work)\n"; +////#endif +// soil_rechr_max = 250.0; +// soil_moist_max = 750.0; +// num_layers = 10; +// +// init_vectors(); +// +// double soil_storage = 375.0; +// double soil_rechr_storage = 125.0; +// +// P = set_default_P(); +// P->is_crhm_test = true; +// S = set_default_S(P->N_Soil_layers); +// S->set_layer_moisture_maximums(*P) +// .set_thermal_conductivities(*P,soil_storage,soil_rechr_storage) +// .set_freezethaw_ratios(*P); +// { +// XG_algorithm XG(0.0,0.0,0.0,*S,*P); +// +// double Zdf_init = 0.0; +// double Zdt_init = 0.0; +// XG.init_freezethaw_degreedays(Zdf_init,Zdt_init,P->Zpf_init); +// } +// +// int start = 0; +// int end = 140000; +// +// for (int i = start; i < end; ++i) +// { +// CRHM crhm(i,reader); +// +// XG_algorithm XG(crhm.hru_tsf,soil_storage,soil_rechr_storage, +// *S,*P); +// +// S->is_newday = i % 24 == 23; +// //S->last_step_new_day = i % 24 == 0; +// XG.run(); +// +// //std::cout << " " << std::endl; +// //std::cout << "Loop: " << i << std::endl; +// //std::cout << "Zdf: " << XG.get_freeze_depth() << std::endl; +// //std::cout << "CRHM Zdf: " << crhm.Zdf << std::endl; +// //std::cout << "Zdt: " << XG.get_thaw_depth() << std::endl; +// //std::cout << "CRHM Zdt: " << crhm.Zdt << std::endl; +// //std::cout << "TrigAcc: " << S->TrigAcc << std::endl; +// //std::cout << "CRHM TrigAcc: " << crhm.TrigAcc << std::endl; +// //std::cout << "TrigState: " << S->TrigState << std::endl; +// //std::cout << "CRHM TrigState: " << crhm.TrigState << std::endl; +// //std::cout << "t_trend: " << S->t_trend << std::endl; +// //std::cout << "CRHM t_trend: " << crhm.t_trend << std::endl; +// //std::cout << "B: " << S->B << std::endl; +// //std::cout << "CRHM B: " << crhm.B << std::endl; +// EXPECT_NEAR(XG.get_thaw_depth(),crhm.Zdt,diff3) << "Loop: " << i; +// EXPECT_NEAR(XG.get_freeze_depth(),crhm.Zdf,diff3) << "Loop: " << i; +// EXPECT_NEAR(S->B,crhm.B,diff4) << "Loop: " << i; +// EXPECT_NEAR(S->TrigAcc,crhm.TrigAcc,diff3) << "Loop : " << i; +// EXPECT_EQ(S->TrigState,crhm.TrigState) << "Loop :" << i; +// soil_storage = reader.getValue("soil_moist",i); +// soil_rechr_storage = reader.getValue("soil_rechr",i); +// }; +//}; diff --git a/src/tests/submoduletests/test_ayers.cpp b/src/tests/submoduletests/test_ayers.cpp new file mode 100644 index 00000000..7fecffa4 --- /dev/null +++ b/src/tests/submoduletests/test_ayers.cpp @@ -0,0 +1,112 @@ +#include "Ayers.hpp" +#include "gtest/gtest.h" +#include "Soil.h" +/* + * CrackTest: Wrapper class for tests + * CrackTest is effectively a mock of Infil_All module but done indirectly. Due to the complexity of the module classes, it was easier to write this. + * The member variables with the _ prefix are inputs that are supplied to the constructor of Crack. + * Default values are given and used for most tests. + * Other member variables are parameters that are also supplied to Crack unless it has the const specifier, then it is just useful for these tests. + * Member functions are just tools to enable the tests. + * Initialization of CrackTest assumes that the frozen period has just begun. + * + */ +class AyersTest : public testing::Test +{ +protected: + + AyersTest() + { + }; + + typedef const Soil::soils_na S; + + typedef Ayers MyAyers; + + double _snowmelt = 1.0; + double _rainfall = 0.0; + std::string _ground_cover = "bare_soil"; + std::string _texture = "coarse_over_coarse"; + + S& soils = Soil::get_soil_obj(); + + MyAyers DoAyers(double& snowmelt, double& rainfall) + { + MyAyers ayers(rainfall,snowmelt,_texture,_ground_cover, soils); + + ayers.run(); + return ayers; + }; + + void DoAssert(MyAyers& ayers,const double inf,const double runoff,const double snowinf) + { + ASSERT_EQ(ayers.get_inf(),inf); + ASSERT_EQ(ayers.get_runoff(),runoff); + ASSERT_EQ(ayers.get_snow_inf(),snowinf); + }; +}; + +TEST_F(AyersTest, ZeroInputs) +{ + _rainfall = 0.0; + _snowmelt = 0.0; + + MyAyers ayers = DoAyers(_snowmelt, _rainfall); + + DoAssert(ayers,0.0,0.0,0.0); +}; + +TEST_F(AyersTest, NonZeroRainfallSmall) +{ + _rainfall = 1e-6; + _snowmelt = 0.0; + + MyAyers ayers = DoAyers(_snowmelt,_rainfall); + + DoAssert(ayers,_rainfall,0.0,0.0); + +}; + +TEST_F(AyersTest, NonZeroBigRainfall) +{ + _rainfall = 1e4; + _snowmelt = 0.0; + + MyAyers ayers = DoAyers(_snowmelt,_rainfall); + + DoAssert(ayers,7.6,_rainfall - 7.6,0.0); + +}; + +TEST_F(AyersTest, NonZeroSnowMelt) +{ + _rainfall = 0.0; + _snowmelt = 1.3e2; + + MyAyers ayers1 = DoAyers(_snowmelt,_rainfall); + + DoAssert(ayers1,_snowmelt,0.0,_snowmelt); + + _snowmelt = 3.14159; + + MyAyers ayers2 = DoAyers(_snowmelt,_rainfall); + + DoAssert(ayers2,_snowmelt,0.0,_snowmelt); + +}; + +TEST_F(AyersTest, MeltAndRain) +{ + _rainfall = 99; + _snowmelt = 13; + + _ground_cover = "small_grains"; + _texture = "medium_over_medium"; + + MyAyers ayers = DoAyers(_snowmelt,_rainfall); + double maxinfil = 10.2; + DoAssert(ayers,maxinfil+_snowmelt,_rainfall - maxinfil,_snowmelt); + +}; + + diff --git a/src/tests/submoduletests/test_crack.cpp b/src/tests/submoduletests/test_crack.cpp new file mode 100644 index 00000000..5b544593 --- /dev/null +++ b/src/tests/submoduletests/test_crack.cpp @@ -0,0 +1,477 @@ +#include "Crack.hpp" +#include "gtest/gtest.h" +#include "CSVreader.hpp" +/* + * CrackTest: Wrapper class for tests + * CrackTest is effectively a mock of Infil_All module but done indirectly. Due to the complexity of the module classes, it was easier to write this. + * The member variables with the _ prefix are inputs that are supplied to the constructor of Crack. + * Default values are given and used for most tests. + * Other member variables are parameters that are also supplied to Crack unless it has the const specifier, then it is just useful for these tests. + * Member functions are just tools to enable the tests. + * Initialization of CrackTest assumes that the frozen period has just begun. + * + */ +class CrackTest : public testing::Test +{ +protected: + + CrackTest() + { + status.init(); + status.begin_freeze(); + set_steps_per_day(seconds_per_hour); + }; + + Crack::info status; + + double _snowmelt = 1.0; + double _rainfall = 0.0; + double _swe = 100.0; + double _soil_storage_at_freeze = 50.0; + double _airtemp = 10.0; + bool _newday = false; + + double major = 5; + double min_swe_to_freeze = 25; + unsigned int infDays = 6; + bool AllowPriorInf = true; + double lenstemp = -10.0; + static constexpr double seconds_per_hour = 3600.0; + double steps_per_day; + + const double diff = 1e-5; + + void set_steps_per_day(double seconds_per_step) + { + steps_per_day = 86400 / seconds_per_step; + }; + + void set_newday(bool newday) + { + _newday = newday; + }; + + // Crack is designed to be declared, initialized and run on every timestep. This function does a single step, while taking a snowmelt as an input. + Crack run_a_step(double snowmelt) + { + Crack model(major,min_swe_to_freeze,infDays,AllowPriorInf,lenstemp,steps_per_day,status); + + _swe -= snowmelt; + model.init_inputs(snowmelt,_rainfall,_swe,_soil_storage_at_freeze,_airtemp,_newday); + + model.run(); + return model; + }; + + // similar to the above but allows the `soil_storage_at_freeze` to be set. + Crack restricted_or_unlimited(double melt, double storage) + { + const double day_melt = melt; + set_newday(true); + _soil_storage_at_freeze = storage; + EXPECT_TRUE(storage == 100.0 || storage == 0.0); + status.daily_melt_total = day_melt; + + return run_a_step(_snowmelt); + }; + + // Expected infiltration for comparison. + double expected_inf(double melt) + { + return std::min(melt * status.index,status.index / infDays * status.init_SWE); + }; + +}; + + +TEST_F(CrackTest, AccumulatesInputsUntilNewDay) { + + for (int i = 0; i < steps_per_day; i++) + { + Crack model = run_a_step(_snowmelt); + + // inf should be 0.0 because the snow-covered period + // just stated and the Crack model applies melt from + // day i as infiltration on day i+1 + EXPECT_EQ(model.get_inf(),0.0); + // checking that snowmelt is accumulating + EXPECT_EQ(status.daily_melt_total, _snowmelt * (i+1)); + } + + set_newday(true); + Crack model = run_a_step(_snowmelt); + // Checks that it restarts properly on a new day + EXPECT_EQ(status.daily_melt_total, _snowmelt); + + +}; +TEST_F(CrackTest, ComputeInfOnDailyTotal) { + set_newday(true); + double yesterday_melt = _snowmelt * steps_per_day / 2.0; + status.daily_melt_total = yesterday_melt; + Crack model = run_a_step(_snowmelt); + + // Check that the initial SWE is recorded correctly + EXPECT_EQ(status.init_SWE,_swe); + + // Computed index by hand (prior to dividing by SWE) + const double index = 36.5924; + + // Check index computation + EXPECT_NEAR(status.index,index/_swe,diff); + + // Check that the major melt count increments + EXPECT_EQ(status.major_melt_count,1); + + // check that maximum per day is properly set + EXPECT_NEAR(status.max_major_per_melt,index / infDays, diff); + + // Computation of inf is conditional. Most likely only one ius true at a time. + EXPECT_NEAR(model.get_snow_inf(),expected_inf(yesterday_melt),diff); +}; + +TEST_F(CrackTest, ConstantInfCheck) +{ +const double start_melt = _snowmelt * steps_per_day; + status.daily_melt_total = start_melt; + + double yesterday; + set_newday(true); + for (int i = 0; i < steps_per_day; i++) + { + + Crack model = run_a_step(_snowmelt); + // checking that infiltration from hour i is equal to hour i-1 (skipping the first hour). + if (i > 0) + EXPECT_EQ(model.get_snow_inf(),yesterday) << "Current hour: " << i; + set_newday(false); + yesterday = model.get_snow_inf(); + } +}; + +TEST_F(CrackTest, MultiDayInfiltration) +{ + const double start_melt = 6.0; + const std::vector melt = {start_melt, 2.0 * start_melt, 3.0 * start_melt}; + for (int j = 0; j < melt.size(); j++) + { + status.daily_melt_total = melt[j]; + set_newday(true); + double inf; + for (int i = 0; i < steps_per_day; i++) + { + Crack model = run_a_step(_snowmelt); + set_newday(false); + + inf = model.get_snow_inf(); + + //Test infiltration on every hour over several major days + EXPECT_NEAR(inf,expected_inf(melt[j]),diff); + //inf and snow_inf should be equal + EXPECT_EQ(inf,model.get_inf()); + } + } +}; + +TEST_F(CrackTest, MajorMeltCounterTest) +{ + const double start_melt = 6.0; + std::vector melt; + const int num_major = 12; + for (int j = 0; j < num_major; j ++) + melt.push_back(start_melt * (j + 1)); + + int num_skip = 7; + melt[num_skip] = 1.0; + for (int j = 0; j < melt.size(); j++) + { + status.daily_melt_total = melt[j]; + set_newday(true); + for (int i = 0; i < steps_per_day; i++) + { + Crack model = run_a_step(_snowmelt); + set_newday(false); + + } + // One day has a non-major, expect major_melt_count to not increment + if (j == num_skip) + EXPECT_EQ(status.major_melt_count,j); + else if (j + 2 < infDays) //Check if incremented + EXPECT_EQ(status.major_melt_count,j+1); + + // checks that infiltration stops after the limit is found + if (j > infDays + 1) + EXPECT_EQ(status.current_inf,0.0) << "Day: " << j+1; + } +}; + +TEST_F(CrackTest, RestrictedTestNewDay) +{ + double day_melt = 48.0; + set_newday(true); + double storage = 100.0; + Crack model = restricted_or_unlimited(day_melt,storage); + + // Check that restricted moves all melt to runoff. + EXPECT_DOUBLE_EQ(model.get_inf(),0.0); + EXPECT_DOUBLE_EQ(model.get_snow_inf(),0.0); + EXPECT_DOUBLE_EQ(model.get_runoff(),day_melt); + EXPECT_DOUBLE_EQ(model.get_melt_runoff(),day_melt); + +}; + +TEST_F(CrackTest,RestrictedTestManyDays) +{ + const double start_melt = 6.0; + _soil_storage_at_freeze = 100.0; + const std::vector melt = {start_melt, 2.0 * start_melt, 3.0 * start_melt}; + for (int j = 0; j < melt.size(); j++) + { + status.daily_melt_total = melt[j]; + set_newday(true); + double inf; + for (int i = 0; i < steps_per_day; i++) + { + Crack model = run_a_step(_snowmelt); + set_newday(false); + + inf = model.get_snow_inf(); + // Same as last, buyt many days. + EXPECT_EQ(inf,0.0); + EXPECT_EQ(model.get_snow_inf(),0.0); + EXPECT_EQ(model.get_runoff(),melt[j]); + EXPECT_EQ(model.get_melt_runoff(),melt[j]); + + } + } + +}; + +TEST_F(CrackTest,UnlimitedTestManyDays) +{ + const double start_melt = 6.0; + _soil_storage_at_freeze = 0.0; + const std::vector melt = {start_melt, 2.0 * start_melt, 3.0 * start_melt}; + for (int j = 0; j < melt.size(); j++) + { + status.daily_melt_total = melt[j]; + set_newday(true); + double inf; + for (int i = 0; i < steps_per_day; i++) + { + Crack model = run_a_step(_snowmelt); + set_newday(false); + + inf = model.get_inf(); + //Same as last but now everything infiltrates because it is unlimited conditions + EXPECT_EQ(inf,melt[j]); + EXPECT_EQ(model.get_snow_inf(),melt[j]); + EXPECT_EQ(model.get_runoff(),0.0); + EXPECT_EQ(model.get_melt_runoff(),0.0); + + } + } + +}; + +TEST_F(CrackTest,IceLensTest) +{ + set_newday(true); + status.daily_melt_total = _snowmelt * steps_per_day; + std::vector T; + + for (int i = 0; i < steps_per_day; i++) + T.push_back(-20.0); + + for (int i = 0; i < steps_per_day;i++) + { + _airtemp = T[i]; + Crack model = run_a_step(_snowmelt); + set_newday(false); + if (i > 0) + EXPECT_EQ(status.tmax,T[0]); + } + set_newday(true); + Crack model = run_a_step(_snowmelt*3.0); + + // Check that the ice-lens is set by increasing major_melt_count safely above InfDays. + EXPECT_EQ(status.major_melt_count,infDays+4); + +}; + +TEST_F(CrackTest,RainOnSnowTest) +{ + status.init(); + status.begin_freeze(); + double rain_on_snow,sum,oldsum; + sum = 0; + + for (int ii = 0; ii < 48; ++ii) + { + _rainfall = 2.0*ii; + sum += _rainfall; + _newday = ii % 24 == 0; + if (_newday) + { + oldsum = sum; + sum = 0.0; + } + Crack model = run_a_step(0.0); + + rain_on_snow = model.get_rain_on_snow(); + + EXPECT_EQ(sum,status.daily_rain_total) << "Step: " << ii; + if (_newday) + EXPECT_EQ(oldsum,rain_on_snow) << "Step: " << ii; + } +}; + +class CrackImplTest : public testing::Test +{ +protected: + + CrackImplTest() + { + status.init(); + }; + CSVReader reader; + Crack::info status; + static constexpr double seconds_per_hour = 3600.0; + double steps_per_day = 24; + double rainfall; + double snowfall; + + double major = 5.0; + double min_swe_to_freeze = 25.0; + unsigned int infDays = 6; + bool AllowPriorInf = true; + double lenstemp = -10.0; + + struct CRHM + { + double infil; + double snowinfil; + double melt_runoff; + double runoff; + double rain_on_snow; + + CRHM(const int& i,CSVReader& reader) + { + infil = reader.getValue("infil",i); + snowinfil = reader.getValue("snowinfil",i) / 24; + melt_runoff = reader.getValue("meltrunoff",i) / 24; + runoff = reader.getValue("runoff",i); + rain_on_snow = reader.getValue("RainOnSnow",i); + }; + }; + + template + void print(std::string text,T val) + { + //std::cout << text << val << std::endl; + }; + + +}; + +TEST_F(CrackImplTest,FullImplementTest) +{ + status.init(); + int start = 0; + int end = 140000; + + double total_rain_on_snow = 0.0; + for (int i = start; i < end; ++i) + { + //std::cout << " " <("net_rain",i); + double snowmelt = reader.getValue("snowmeltD",i) / 24; + double swe = reader.getValue("SWE",i); + double soil_storage_at_freeze = 50; + double airtemp = reader.getValue("hru_t",i); + bool crackon = reader.getValue("crackon",i); + std::string datetime = reader.getValue("datetime",i); + //std::cout << datetime << std::endl; + CRHM crhm(i,reader); + print("rainfall: ", rainfall); + print("snowmelt: ", snowmelt); + print("swe: ", swe); + + double runoff = 0.0; + double melt_runoff = 0.0; + double inf = 0.0; + double snowinf = 0.0; + double rain_on_snow = 0.0; + bool is_day_over = i % 24 == 23; + print("Day over tracker: ", i % 24); + + if ((swe > min_swe_to_freeze && !status.frozen && is_day_over) || (crackon && i == start) ) + { + status.begin_freeze(); + status.end_freeze_tomorrow = false; + } + + print("snowinf: ", snowinf); + print("inf: ", inf); + print("melt_runoff: ", melt_runoff); + print("runoff: ", runoff); + print("frozen: ",status.frozen); + print("major melt count: ",status.major_melt_count); + print("Total input: ", snowmelt+rainfall); + + if (status.frozen) + { + Crack crack(major, min_swe_to_freeze, infDays, + AllowPriorInf, lenstemp,steps_per_day,status); + + crack.init_inputs(snowmelt, rainfall, swe, soil_storage_at_freeze, + airtemp, is_day_over); + status.daily_melt_total = snowmelt * steps_per_day; + crack.is_CRHM_compare_test = true; + crack.run(); + + runoff = crack.get_runoff() / steps_per_day; + melt_runoff = crack.get_melt_runoff() / steps_per_day; + inf = crack.get_inf() / steps_per_day; + snowinf = crack.get_snow_inf() / steps_per_day; + rain_on_snow = crack.get_rain_on_snow(); + + total_rain_on_snow += rain_on_snow; + + if (is_day_over && swe <= 0.0 && status.major_melt_count > 0) + status.end_freeze(); + //if (swe <= 0.0 && status.major_melt_count > 0) + // status.end_freeze_tomorrow = true; + + } + else + { + crhm.infil = 0.0; + crhm.snowinfil = 0.0; + crhm.melt_runoff = 0.0; + crhm.runoff = 0.0; + //crhm.rain_on_snow = 0.0; + }; + print("Melt total: ",status.daily_melt_total); + print("Rain total: ",status.daily_rain_total); + print("index: ", status.index); + print("Max major per melt: ", status.max_major_per_melt); + print("init_SWE", status.init_SWE); + print("CRHM: rain on snow:",crhm.rain_on_snow); + + print("snowinf: ", snowinf); + print("inf: ", inf); + print("melt_runoff: ", melt_runoff); + print("runoff: ", runoff); + double diff = 1e-5; + EXPECT_NEAR(crhm.infil,inf - snowinf,diff) << "Step: " << i; + EXPECT_NEAR(crhm.snowinfil,snowinf,diff) << "Step: " << i; + EXPECT_NEAR(crhm.melt_runoff,melt_runoff,diff) << "Step: " << i; + EXPECT_NEAR(crhm.runoff,runoff - melt_runoff,diff) << "Step: " << i; + EXPECT_EQ(status.frozen,crackon) << "Step: " << i; + EXPECT_NEAR(crhm.rain_on_snow,total_rain_on_snow,diff*10) << "Step: " << i; + + + } +}; diff --git a/src/tests/submoduletests/test_daily_accumulator.cpp b/src/tests/submoduletests/test_daily_accumulator.cpp new file mode 100644 index 00000000..1dc990c1 --- /dev/null +++ b/src/tests/submoduletests/test_daily_accumulator.cpp @@ -0,0 +1,146 @@ +#include +#include "daily_accumulator.hpp" + +class data +{ +public: + bool is_new_day(); + int steps_per_day(); + + void set_steps(int& step) + { _steps_per_day = step; }; + + int get_steps() + { return _steps_per_day; }; + + void step_forward() + { ++_count; }; + + void reset() + { _count = 0; }; + +private: + int _steps_per_day; + int _count = 0; +}; + +int data::steps_per_day() +{ + return get_steps(); +}; + +bool data::is_new_day() +{ + return _count % _steps_per_day == 0 && _count > 0; +}; + +class DailyAccumulatorTest : public testing::Test +{ +protected: + template + double sum_to_step(int i,int steps_per_day,T arr) + { + int start = nearest_day_start(i,d.get_steps()); + double sum = 0; + for (int j = start; j < i; ++j) + { + sum += arr.at(j); + } + return sum; + }; + + int nearest_day_start(int i, int step_per_day) + { + if (step_per_day == 0) {return 0;} + + int day_start_index = (i / step_per_day) * step_per_day; + + return (day_start_index >= step_per_day) ? day_start_index : 0; + }; + + data d; +}; + +TEST_F(DailyAccumulatorTest,HelperFunctionSumToStepTest) +{ + EXPECT_EQ(1,1); + const int steps = 5; + + std::vector T{-1.0,0.0,2.0,3.0,6.2}; + + // sum to step test + double sum = sum_to_step(0,steps,T); + + ASSERT_EQ(sum,0.0); + + sum = sum_to_step(2,steps,T); + + ASSERT_EQ(sum,-1.0); + + sum = sum_to_step(5,steps,T); + + ASSERT_EQ(sum,10.2); +}; + +TEST_F(DailyAccumulatorTest,HelperFunctionNearestDayStartTest) +{ + int start = nearest_day_start(3,12); + + ASSERT_EQ(start,0); + + start = nearest_day_start(99,400); + + ASSERT_EQ(start,0); + + start = nearest_day_start(29,12); + + ASSERT_EQ(start,24); + + start = nearest_day_start(35,4); + + ASSERT_EQ(start,32); + + start = nearest_day_start(35+4,4); + + ASSERT_EQ(start,36); +}; + +TEST_F(DailyAccumulatorTest,Test) +{ + double temperature = 0.0; + data d; + const int steps = 4; + int steps_copy = steps; + d.set_steps(steps_copy); + + daily_accumulator accumulate_temperature(d); + accumulate_temperature.bind_to_var(temperature); + + std::array T{-1.0,-3.2,4.0,1.2,-2.9,-10.1,-5.1,5.0}; + + for (int i = 0; i < 2 * steps; ++i) + { + temperature = T[i]; + accumulate_temperature.execute(); + d.step_forward(); + std::cout << "test: " << i % steps << "," << i << std::endl; + if (i < steps) + { + ASSERT_EQ(accumulate_temperature.get_last_mean(),0.0) << "Loop: " << i; + } + else if (i % steps == 0) + { + std::cout << "har" <