diff --git a/docs/submodule_instructions.rst b/docs/submodule_instructions.rst new file mode 100644 index 00000000..6208d008 --- /dev/null +++ b/docs/submodule_instructions.rst @@ -0,0 +1,447 @@ +How To: Submodule Framework +============================ + +The CHM modules, for good reason, allow the user to design their module however they'd like as long as they couple the calculations to the CHM coupler that coordinates the modules and the data passing between them. Modules can be subdivided into two categories: + +1. External models coupled to CHM. +2. Models coded in CHM directly. + +Being able to accommodate both is one of CHMs biggest strengths (the same can be said for most modular models). Modules of the first type are only possible when the models are reasonably portable and adaptable to any framework. It can be readily seen by reading these modules that a lot of work is required in order to make it mesh (pun intended) with CHM. And ultimately that is not surprising. It was not designed with CHM in mind. However, given that we have full control over CHM modules, it makes sense to design a consistent framework for writing them. And that framework should go beyond implementing the virtual functions ``void module_base::run(mesh_elem&)`` and ``void module_base::init(mesh&)``, and using a class derived from ``face_info`` to store per-triangle data. + +If no framework - or *encouraged* - design exists, then modules will either become so completely different that learning each is equivalent to learning a new model each time, or they will be coded in such a manner that testability, extensibility, flexibility, and readability fall to the back burner. + +In this document, I'll present basic instructions and guidelines for designing submodules and then how to write (or modify) a module that uses them. + +Designing Submodules +-------------------- + +Avoiding the dirty details of the why, here is the how. + +Step 1 +^^^^^^ + +Open a new file in ``src/modules/submodules`` with a descriptive name. For this we will use ``submodule_example.hpp``. In this file, immediately create a namespace with the same name as the file and nest inside of it a class called ``Model`` derived from ``base_step,Data>`` where ``Data`` is a template parameter. More on this later... + +The ``Model`` class should have one public function: ``execute_impl(Data& d)``. It will look something like this: + +.. code-block:: cpp + + #include "base_step.hpp" + + namespace submodule_example + { + template + class Model : public base_step,Data> + { + public: + void execute_impl(Data& d); + }; + }; + +Step 2 +^^^^^^ + +Take stock of your calculation and write out a list of inputs it needs and outputs it provides (similar to CHM). Then write a concept to enforce that ``Data`` has getters and setters for both of these. + +.. code-block:: cpp + + #include "base_step.hpp" + #include + + namespace submodule_example + { + template + concept submodule_data = requires(T& t) + { + { t.input1(); } -> std::floating_point; + { t.input2(); } -> std::floating_point; + + { t.output1(std::declval()); } -> std::same_as; + { t.output2(std::declval()); } -> std::same_as; + }; + + template + class Model : public base_step,Data> + { + public: + void execute_impl(Data& d); + }; + }; + +I've left out implementing ``execute_impl(Data&)`` on purpose for now. + +Step 3 (optional) +^^^^^^^^^^^^^^^^^ + +Immediately write tests. Navigate to ``src/tests/submoduletests/`` and create a file with a good name like ``test_submodule_example.cpp``. The body will look something like this: + +.. code-block:: cpp + + #include + #include "submodule_example.hpp" + + class data_for_test + { + // Implement example data class to satisfy the concept above + }; + + class SubmoduleExampleTest : public ::testing::Test + { + protected: // Must be protected!!! + submodule_example submodule; + }; + + TEST_F(SubmoduleExampleTest,TestZerosReturnsZeros) + { + // TEST_F is a macro to create a class derived from SubmoduleExampleTest. So the `submodule` object already exists. + // Operate on that are use macros like EXPECT_EQ(expected,produced); to check if the submodule is working as expected. + // Test BEHAVIOURS. + } + +If you are having trouble designing a test as the result of designing increasingly contrived situations to test obscure edge cases that will actually compile, or you feel yourself wanting to test private methods, this is a `code smell `_. I would recommend then further encapsulating these calculations into another class, defined within the ``submodule_example`` namespace, with a public, testable interface. Then test that. If you had an ``std::vector`` private member, you wouldn't test the ``std::vector`` through the public interface, you'd assume (rightly) it was tested elsewhere. + +Why write tests first +^^^^^^^^^^^^^^^^^^^^^ + +It's called Test Driven Development. It exists to speed up development by spending time early on tests and later benefiting from the pre-planning. A common expression I've read is "If you have a bug, that means there is a missing test". + +Step 4 +^^^^^^ + +Return to ``submodule_example.hpp`` and iterate until it passes your tests. Tada! Now you know its working. Don't forget to add ``test_submodule_example.cpp`` to your ``CMakeLists.txt`` and enable the testing flag. I prefer to compile once, then open ``CMakeCache.txt`` in my build directory and change ``BUILD_TESTS`` to ``ON``, then recompile. It saves having to deal with changing the ``CMakeLists.txt`` every time while trying to commit changes to the actual source later. + +**You're done!** You've successfully written a submodule. The next step is to include it in a CHM module. + +Designing Modules with Submodules +--------------------------------- + +For simplicity, we will assume that we are starting from scratch to build this module, but these instructions can easily be applied to an existing module, either to change a calculation to a submodule or add a submodule as a separate calculation. Write the header as normal with a nice an descriptive name, here we will use ``module_example.hpp``. The way we do this is make use of the class ``data_base`` for automatic caching and some other tools. + +Step 1: Write the Cache +^^^^^^^^^^^^^^^^^^^^^^^ + +First, you need a struct (classes are OK) to act as the cache. The short of it is that you require a member for outputs and a member for inputs. Therefore our Cache is + +.. code-block:: cpp + + struct Cache : public cache_base + { + double input1 = default_value(); + double input2 = default_value(); + + double output1 = 0.0; + double output2 = 0.0; + }; + +``default_value()`` comes from the ``cache_base`` parent class. The details don't matter too much here, but it initializes the input as a sentinel value so that each member can separately be lazy-initialized. + +Step 2: Write the data class +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The ``data`` class is the class that holds the per-triangle data for a module between time steps. This is the class that will be passed as a template parameter to the submodule. Therefore, it must satisfy the concept we defined in ``submodule_example.hpp``. In order to use the cache, it also must be derived from ``data_base``. The Cache as a template parameter sets the type of the cache used behind the scenes. + +Lazy initialization of the input variables on the cache is only possible if ``data`` has a reference to the ``mesh_elem`` object it corresponds to, the ``config_file`` object, and the ``global`` shared pointer. So ``data_base`` has a constructor to pass these references. therefore, the data class is as follows: + +.. code-block:: cpp + + class data : data_base + { + public: + double input1() const; + double input2() const; + + double output1(const double out_val) const; + double output2(const double out_val) const; + + using data_base::data_base; + // Alternatively, you can write out the constructor as follows: + // data(mesh_elem& face_in, std::shared_ptr param, config_file& cfg) + // : data_base(face_in,param,cfg); + }; + +Step 3: Write the Module header +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Put it all together. I recommend placing the ``data`` and ``Cache`` classes inside the body of the module. Like so: + +.. code-block:: cpp + + // Essential CHM includes here + #include "submodule_example.hpp" + + class module_example : public module_base + { + public: + void run(mesh_elem&) override; + void init(mesh&) override; + + struct Cache : public cache_base + { + double input1 = default_value(); + double input2 = default_value(); + + double output1 = 0.0; + double output2 = 0.0; + }; + + class data : data_base + { + public: + double input1() const; + double input2() const; + + double output1(const double out_val) const; + double output2(const double out_val) const; + + using data_base::data_base; + }; + private: + submodule_example submodule_obj; + } + +Step 4: Implement the ``data`` member functions +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Switching to ``module_example.cpp``, we implement the auto-caching using the protected member function from the ``data_base`` class for the inputs and another protected member for outputs. These functions automatically lazy initialize the cache. I will admit here that this is possibly the ugliest part of the implementation and could use improvements. Nonetheless, here is how you do it. You call the function ``update_value`` which expects two arguments, each argument must be a callable. I intend for these to be lambda functions but it may work with function pointers, but I have not tested that. The first lambda returns a reference to the cache member to be accessed, and the second lambda returns a copy of the dereferenced ``face`` object. ``update_value`` calls the second lambda if the cache has not be initialized yet or if the input variable is equal to ``default_value()`` that we initialized the inputs to in the ``module_example::Cache``. + +.. code-block:: cpp + + double module_example::data::input1() + { + update_value( [this]() -> auto& { return cache_->input1; }, + [this]() { return (*face)["input1"_s];} + ); + + return cache_->input1; + } + + double module_example::data::input2() + { + update_value( [this]() -> auto& { return cache_->input2; }, + [this]() { return (*face)["input2"_s];} + ); + + return cache_->input2; + } + +Note that in all four lambdas, ``this`` is captured. ``This`` allows the use of ``face`` and ``cache_`` objects. I means that both lambdas have a copy of the ``this`` pointer (not the whole class) so the overhead is minimal and only stored during the call to ``update_value`` and freed immeidately after the function call. + +``update_value`` is a templated function and therefore the compiler generates a unique function for each unique pair of lambdas. In addition, the arguments are r-value references so these temporary lambdas are passes correctly. The first lambda **must** return a reference. Luckily, concepts hidden in ``data_base`` ensure this is the case, so it won't compile if you forget. This allows the first lambda to be assigned. + +Now for the output functions, its very similar and uses a function called ``set_output`` that likewise accepts two arguments but in this case, only the first argument is a lambda reference to the cache member, and the second is simply forwarding the function argument. + +.. code-block:: cpp + + double module_example::data::output1(const double out_val) + { + set_output( [this]() -> auto& { return cache_->output1; }, out_val); + }; + + double module_example::data::output2(const double out_val) + { + set_output( [this]() -> auto& { return cache_->output2; }, out_val); + }; + +Other kinds of inputs +^^^^^^^^^^^^^^^^^^^^^ + +Not all inputs are created equal. Some are variables, accessed by a dereference from ``face``, but others might be per-triangle or domain-wide parameters. You can approach this however you'd like. But remember that ``cfg`` and ``global`` and ``face`` are all available. So... make use of them. + +For domain-wide parameters, I recommend adding a static variable to the data class, and assigning it in the ``module_example::init`` function, like so: + +.. code-block:: cpp + + // .hpp + + data : public data_base + { + public: + // other stuff + static double param1; + }; + + // .cpp + + void module_example::init(mesh& domain) + { + data::param1 = cfg.get("param1"); + }; + + static double module_example::data::param1() // Doesn't modify the state and only uses static members, so this can be static. + { + return param1; + }; + + // Another option + double module_example::data::param2() + { + //thread safe + //initialized once and then never again + //techincally has a slight overhead to check if its been set yet. + // Doesn't pollute the public interface of the data class. + static double param2 = cfg.get("param2"); + return param2; + } + +For triangle specific parameters, I would proceed exactly the same way as ``param1`` without the static keyword. + +.. code-block:: cpp + + // .hpp + + data : public data_base + { + public: + // other stuff + double param1; + }; + + // .cpp + + void module_example::init(mesh& domain) + { + // loop over every face + // Get the local instance of data via make_module_data, skipped here because I've modified it. + param1 = face->veg_parameter("param1"_s); + }; + + double module_example::data::param1() + { + return param1; + }; + + // Another option + // Static variagble in a functin cannot work because + // it is per CLASS not per instance and you need a unique value per triangle + // However, you could do the following: + + module_example::data::data(mesh_elem& face_in,std::shared_ptr param,config_file& cfg) + : data_base(face_in,param,cfg) + { + param1 = face->veg_parameter("param1"_s); + param2 = cfg.get("param2"_s); // here param2 is a static member of data. + }; + +In fact, that last example of setting them in the constructor is a bit genius and I actually haven't done it yet myself... welp guess I'll start doing that. + +The final type is something accessed via the ``global`` instance ``param``. These will often have to be called dynamically because things change. If it is a time step, set it in the constructor for ``data``, but if it is the time of day then one must instead set it at each call. If one wants to avoid calling ``global`` more than once per time step, feel free to add it as a cached member of ``Cache``. + +.. code-block:: cpp + + double module_example::data::hour() + { + // or whatever other function + return param->hour(); + } + +Step 5: Implement the module +^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Most proceeds as normal. Constructor sets the depends/provides, destructor is empty unless you were evil and used ``new``, be sure to ``delete`` those. + +If we use the method to set parameters in the constructor, the ``.cpp`` file will look like this: + +.. code-block:: cpp + + #include "module_example.hpp" + REGISTER_MODULE_CPP(module_example); + + module_example::module_example(config_file cfg) : module_base("module_example", parallel::data, cfg) + { + depends("input1"); + depends("input2"); + + provides("output1"); + provides("output2"); + } + + //empty deconstructor + + void module_example::init(mesh& domain) + { + // Recommend setting these outside the loop and not in the constructor. + // You create it so many times, it is unwise to set it each time, but + // overhead might be minimal + data::static_parameter = cfg.get("static_parameter"); + for (size_t i = 0; i < domain->size_local_faces(); ++i) + { + auto face = domain->face(i); + face->make_module_data(ID,face,global_param,cfg); + // Everything is set in the constructor... so no need do actually store a reference + // there may be some situations where a parameter determines how something is set + // In that case, store a reference to data and make whatever assignments you'd like. + } + } + + void module_example::run(mesh_elem& face) + { + auto& d = face->get_module_data(ID); + + submodule_example(d); + + d.set_module_outputs(); + } + + void module_example::data::set_module_outputs() + { + auto& c = get_cache(); + if (c) + { + (*face)["output1"_s] = (c->output1 == default_value()) ? c.output1 : 0.0; + (*face)["output2"_s] = (c->output2 == default_value()) ? c.output2 : 0.0; + } + else + { + (*face)["output1"_s] = 0.0; + (*face)["output2"_s] = 0.0; + } + } + +And now your module using submodules is complete! Notice that a function called ``make_module_data`` is called in the module virtual function ``init`` and its constructor now accepts more than just the module ``ID`` as an argument. I've rewritten just the ``make_module_data`` function to have a template of variadic parameters to pass to the constructor of the ``data`` object. + +Final Notes +----------- + +There are a few things that could be improved: + +1. The ``update_value`` and ``set_output`` protected ``data_base`` member functions. They are a bit ugly. + + One idea is to set these lambdas in the constructor and store them as std::functions. Not ideal really but doable. + +2. The style of setting input cache members is not ideal. ``double input1 = default_value();`` is redundant and ugly. + + One idea here is to create an ``input`` class that takes a parameter to a primitive type and auto sets it to ``default_value()``. Like so: + + .. code-block:: cpp + + class Cache : public cache_base + { + input input1; + input input2; + + output output1; // initialized to 0.0 + output output2; + } + + Perhaps ``input`` and ``output`` classes are where we could store lambdas for note 1 above? + +3. When you have many state variables, I recommend defining a struct like: + + .. code-block:: cpp + + namespace submodule_example + { + struct State + { + //state members here + }; + + // modify concept to expect a get_state function as a reference. + // Store a State object privately on data. + } + +Then, if during the testing you decide to break-up your submodule into several subclasses you can design them to operate with the ``State`` object, rather than on the ``data`` object. This nifty trick means you can test these classes (defined in the ``submodule_example`` namespace) without having to worry about the concept to enforce on the template parameter of ``submodule_example``. This also means that these classes can be compiled in a ``submodule_example.cpp`` file. + +4. Consider a situation with N modules. Now you'll have a very busy ``data`` class. Depending on what functions are enforced by their respective concepts, you could have multiple functions for air temperature: ``air_temp()``, ``air_temperature()``, ``temp()``, and so on. For these, implement one and have the others simply pass through to the implemented version. + + Likewise, only store a single cache member for all the air temperatures! Unit conversions can be done in the respective calls, if one module expects Celsius but the other Kelvin! diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ecd6505e..db665119 100755 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -62,12 +62,21 @@ set(MODULE_SRCS modules/PBSM3D.cpp modules/snobal.cpp modules/fsm.cpp + modules/katabatic_routing_glacier.cpp modules/testing/mpi_example.cpp CACHE INTERNAL "" FORCE ) +set(SUBMODULE_SRCS + # modules/submodules/water_flux.cpp + modules/submodules/Glacier.cpp + modules/submodules/katabatic_melt_energy.cpp + modules/submodules/melt_routing_glacier.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 +86,7 @@ set(CHM_SRCS metdata.cpp physics/Atmosphere.cpp + physics/PhysConst.cpp mesh/triangulation.cpp mesh/ugrid_writer.cpp @@ -110,6 +120,7 @@ set (HEADER_FILES modules/interp_met modules/snobal modules/testing + modules/submodules math libmaw interpolation @@ -234,6 +245,7 @@ add_executable( ${CHM_SRCS} ${FILTER_SRCS} ${MODULE_SRCS} + ${SUBMODULE_SRCS} ) @@ -412,17 +424,23 @@ if (BUILD_TESTS) message(STATUS "Tests enabled. Run with make check") 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 # test_daily.cpp - tests/test_triangulation.cpp + # tests/test_triangulation.cpp + tests/submoduletests/test_data_base.cpp + #tests/submoduletests/test_water_flux.cpp + tests/submoduletests/test_glacier.cpp + tests/submoduletests/test_katabatic_melt_energy.cpp + tests/submoduletests/test_melt_routing.cpp + tests/submoduletests/test_bisection.cpp tests/main.cpp ) @@ -435,6 +453,7 @@ if (BUILD_TESTS) ${CHM_SRCS} ${FILTER_SRCS} ${MODULE_SRCS} + ${SUBMODULE_SRCS} ${TEST_SRCS} ) set_target_properties(runUnitTests diff --git a/src/core.cpp b/src/core.cpp index b0390cdd..b2f80af1 100644 --- a/src/core.cpp +++ b/src/core.cpp @@ -1155,10 +1155,17 @@ void core::config_output(pt::ptree &value) OGRCoordinateTransformation::DestroyCT(coordTrans); } - if(!_mesh->is_geographic()) - out.face = _mesh->locate_face(out.x, out.y); - else - out.face = _mesh->locate_face(out.longitude, out.latitude); + + + // if(!_mesh->is_geographic()) + // { + // out.face = _mesh->locate_face(out.x, out.y); + // } + // else + // { + // out.face = _mesh->locate_face(out.longitude, out.latitude); + // } + out.face = _mesh->face(0); if(out.face != nullptr) { diff --git a/src/math/Bisection.hpp b/src/math/Bisection.hpp new file mode 100644 index 00000000..8f560841 --- /dev/null +++ b/src/math/Bisection.hpp @@ -0,0 +1,51 @@ +#pragma once + +#include +#include +#include +namespace math +{ + enum class Root { Found, DidNotConverge, NoRootGuaranteed }; + struct Result + { + Root root; + double value; + size_t iterations; + }; + + inline Result bisection( + std::function f, + double low_guess, + double high_guess, + double tol = 1e-10, + size_t max_iter = 1000 + ) { + double fa = f(low_guess), fb = f(high_guess); + + if (fa * fb > 0) + return Result { + Root::NoRootGuaranteed, + std::numeric_limits::quiet_NaN(), + 0}; + + for (size_t i = 0; i < max_iter; ++i) { + double c = low_guess + (high_guess - low_guess) / 2; + double fc = f(c); + + if (std::abs(fc) < tol || (high_guess - low_guess) / 2 < tol) + return Result{Root::Found,c,i} ; + + if (fa * fc < 0) { + high_guess = c; + fb = fc; + } else { + low_guess = c; + fa = fc; + } + } + + return Result{Root::DidNotConverge, + std::numeric_limits::quiet_NaN(), + max_iter}; // Did not converge + }; +}; diff --git a/src/mesh/triangulation.hpp b/src/mesh/triangulation.hpp index 517df2e5..2f5346b5 100644 --- a/src/mesh/triangulation.hpp +++ b/src/mesh/triangulation.hpp @@ -436,9 +436,10 @@ class face : public Fb // void set_module_data(const std::string &module, face_info *fi); - template - T& make_module_data(const std::string &module); - + // Overload for module data constructor + template + T& make_module_data(const std::string &module, Args... args); + std::string _debug_name; //for debugging to find the elem that we want int _debug_ID; //also for debugging. ID == the position in the output order, starting at 0 size_t cell_global_id; @@ -1788,22 +1789,22 @@ timeseries::iterator face::now() } +// Overloaded for construtor that requires arguments template < class Gt, class Vb> -template -T& face::make_module_data(const std::string &module) +template +T& face::make_module_data(const std::string &module, Args... args) { //we don't already have this, make a new one. if(!_module_face_data[module]) { // T* data = new T; - _module_face_data[module] = std::make_unique(); + _module_face_data[module] = std::make_unique(args...); } return get_module_data(module); } - template < class Gt, class Fb> template < typename T> T& face::get_module_data(const std::string &module) diff --git a/src/modules/PBSM3D.cpp b/src/modules/PBSM3D.cpp index 90b4ef0e..bd545a18 100644 --- a/src/modules/PBSM3D.cpp +++ b/src/modules/PBSM3D.cpp @@ -423,7 +423,7 @@ void PBSM3D::run(mesh& domain) auto tol = [](double a, double b) -> bool { return fabs(a - b) < 1e-8; }; // ice density - double rho_p = PhysConst::rho_ice; + double rho_p = PhysConst::rho_ice(); #pragma omp for for (size_t i = 0; i < domain->size_local_faces(); i++) { @@ -700,7 +700,7 @@ void PBSM3D::run(mesh& domain) // c_4 = 0.5; // g = 9.81; - return u2 * PhysConst::kappa / log(2.0 / (0.6131702345e-2 * ustar * ustar + .5 * lambda)) - + return u2 * PhysConst::kappa() / log(2.0 / (0.6131702345e-2 * ustar * ustar + .5 * lambda)) - ustar; }; try @@ -717,7 +717,7 @@ void PBSM3D::run(mesh& domain) else { // follow PBSM (Pom & Li 2000; Alpine3D) and don't calculate the feedback of z0 on u* - ustar = u2 * PhysConst::kappa / log(2.0 / 0.0002); + ustar = u2 * PhysConst::kappa() / log(2.0 / 0.0002); } if (ustar >= u_star_saltation_threshold) @@ -746,7 +746,7 @@ void PBSM3D::run(mesh& domain) { // we still need a u* for spatial K estimation later d.z0 = Snow::Z0_SNOW; - ustar = std::max(0.01, PhysConst::kappa * uref / log(Atmosphere::Z_U_R / d.z0)); + ustar = std::max(0.01, PhysConst::kappa() * uref / log(Atmosphere::Z_U_R / d.z0)); } // sanity checks @@ -1153,7 +1153,7 @@ void PBSM3D::run(mesh& domain) } } // Li and Pomeroy 2000 - double l = PhysConst::kappa * (cz + d.z0) * l__max / (PhysConst::kappa * (cz + d.z0) + l__max); + double l = PhysConst::kappa() * (cz + d.z0) * l__max / (PhysConst::kappa() * (cz + d.z0) + l__max); if (debug_output) (*face)["l"_s] = l; diff --git a/src/modules/Simple_Canopy.cpp b/src/modules/Simple_Canopy.cpp index 75dd5509..617ef176 100644 --- a/src/modules/Simple_Canopy.cpp +++ b/src/modules/Simple_Canopy.cpp @@ -172,20 +172,20 @@ void Simple_Canopy::run(mesh_elem &face) // Canopy temperature is first approximated by the air temperature. double T1 = ta + mio::Cst::t_water_freezing_pt; // Canopy temperature (C to K) - double rho = air_pressure*1000/(PhysConst::Rgas*T1); // density of Air (pressure kPa to Pa = *1000) + double rho = air_pressure*1000/(PhysConst::RgasDry()*T1); // density of Air (pressure kPa to Pa = *1000) double U1 = U_R; // Wind speed (m/s) at height Z_vw [m] (top of canopy) // Aerodynamic resistance of canopy - ra = (log(Zref/Z0snow)*log(Zwind/Z0snow))/pow(PhysConst::kappa,2)/U1; // (s/m) + ra = (log(Zref/Z0snow)*log(Zwind/Z0snow))/pow(PhysConst::kappa(),2)/U1; // (s/m) - double deltaX = 0.622*PhysConst::Ls*Qs(air_pressure, T1)/(PhysConst::Rgas*(pow(T1,2))); // Must be (kg K-1) + double deltaX = 0.622*PhysConst::Ls()*Qs(air_pressure, T1)/(PhysConst::RgasDry()*(pow(T1,2))); // Must be (kg K-1) double q = (rh/100)*Qs(air_pressure, T1); // specific humidity (kg/kg) // snow surface temperature of snow in canopy - Ts = T1 + (Snow::emiss*(ilwr - PhysConst::sbc*pow(T1, 4.0)) + PhysConst::Ls*(q - Qs(air_pressure, T1))*rho/ra)/ - (4.0*Snow::emiss*PhysConst::sbc*pow(T1, 3.0) + (PhysConst::Cp + PhysConst::Ls*deltaX)*rho/ra); + Ts = T1 + (Snow::emiss*(ilwr - PhysConst::sbc()*pow(T1, 4.0)) + PhysConst::Ls()*(q - Qs(air_pressure, T1))*rho/ra)/ + (4.0*Snow::emiss*PhysConst::sbc()*pow(T1, 3.0) + (PhysConst::Cp() + PhysConst::Ls()*deltaX)*rho/ra); Ts -= mio::Cst::t_water_freezing_pt; // K to C @@ -227,7 +227,7 @@ void Simple_Canopy::run(mesh_elem &face) Kstar_H = iswr * (1.0 - Alpha_c - Tauc * (1.0 - Albedo)); // what is Kstar_H??? // Incident long-wave at surface, "(W/m^2)" - Qlisn = ilwr * Vf_ + (1.0 - Vf_) * Vegetation::emiss_c * PhysConst::sbc * pow(T1, 4.0) + B_canopy * Kstar_H; + Qlisn = ilwr * Vf_ + (1.0 - Vf_) * Vegetation::emiss_c * PhysConst::sbc() * pow(T1, 4.0) + B_canopy * Kstar_H; // Incident short-wave at surface, "(W/m^2)" Qsisn = iswr * Tauc; @@ -289,7 +289,7 @@ void Simple_Canopy::run(mesh_elem &face) Kd = iswr * (1.0 - Alpha_c - Tau_b_gap * (1.0 - Albedo)); Qlisn = Vgap * ilwr + (1.0 - Vgap) * ((ilwr * Tau_b_gap + - (1.0 - Tau_b_gap) * Vegetation::emiss_c * PhysConst::sbc * + (1.0 - Tau_b_gap) * Vegetation::emiss_c * PhysConst::sbc() * pow(T1, 4.0f)) + B_canopy * Kd); Qsisn = cosxs * Qdfo * Tau_b_gap + Vgap * (iswr - Qdfo) + (1.0 - Vgap) * Tau_d * (iswr - Qdfo); @@ -373,7 +373,7 @@ void Simple_Canopy::run(mesh_elem &face) double Es = 611.15 * exp(22.452 * ta / (ta + 273.0)); // {sat pressure} - double SvDens = Es * PhysConst::M / (PhysConst::R * (ta + 273.0)); // {sat density} + double SvDens = Es * PhysConst::M() / (PhysConst::R() * (ta + 273.0)); // {sat density} Lamb = 6.3e-4 * (ta + 273.0) + 0.0673; // thermal conductivity of atmosphere Nr = 2.0 * Snow::Radius * uVent / Atmosphere::KinVisc; // Reynolds number @@ -381,19 +381,19 @@ void Simple_Canopy::run(mesh_elem &face) SStar = M_PI * pow(Snow::Radius, 2) * (1.0 - Snow::AlbedoIce) * iswr; // SW to snow particle !!!! changed A1 = Lamb * (ta + 273) * Nu; - B1 = PhysConst::Ls * PhysConst::M / (PhysConst::R * (ta + 273.0)) - 1.0; + B1 = PhysConst::Ls() * PhysConst::M() / (PhysConst::R() * (ta + 273.0)) - 1.0; J = B1 / A1; Sigma2 = rh / 100 - 1; D = 2.06e-5 * pow((ta + 273.0) / 273.0, -1.75); // diffusivity of water vapour C1 = 1.0 / (D * SvDens * Nu); Alpha = 5.0; - Mpm = 4.0 / 3.0 * M_PI * PhysConst::rho_ice * pow(Snow::Radius, 3) * + Mpm = 4.0 / 3.0 * M_PI * PhysConst::rho_ice() * pow(Snow::Radius, 3) * (1.0 + 3.0 / Alpha + 2.0 / pow(Alpha, 2)); // sublimation rate of single 'ideal' ice sphere: - double Vs = (2.0 * M_PI * Snow::Radius * Sigma2 - SStar * J) / (PhysConst::Ls * J + C1) / Mpm; + double Vs = (2.0 * M_PI * Snow::Radius * Sigma2 - SStar * J) / (PhysConst::Ls() * J + C1) / Mpm; // snow exposure coefficient (Ce): @@ -409,7 +409,7 @@ void Simple_Canopy::run(mesh_elem &face) // calculate 'ice-bulb' temperature of intercepted snow: - double IceBulbT = ta - (Vi * PhysConst::Ls / 1e6 / PhysConst::Ci); + double IceBulbT = ta - (Vi * PhysConst::Ls() / 1e6 / PhysConst::Ci()); // determine whether canopy snow is unloaded: @@ -431,8 +431,8 @@ void Simple_Canopy::run(mesh_elem &face) // limit sublimation to canopy snow available and take sublimated snow away from canopy snow at timestep start //Subl_Cpy = -data.Snow_load*Vi*Hs*Global::Interval*24*3600/Hs; // make W/m2 (original in CRHM) - Subl_Cpy = -data.Snow_load * Vi * PhysConst::Ls * global_param->dt() / - PhysConst::Ls; // make W/m2 TODO: check Interval is same as dt() (in seconds + Subl_Cpy = -data.Snow_load * Vi * PhysConst::Ls() * global_param->dt() / + PhysConst::Ls(); // make W/m2 TODO: check Interval is same as dt() (in seconds // TODO: Hs/HS = 1 !!! (in CRHM, kept here for conistency...) if (Subl_Cpy > data.Snow_load) { @@ -719,4 +719,4 @@ void Simple_Canopy::load_checkpoint(mesh& domain, netcdf& chkpt) } -} \ No newline at end of file +} diff --git a/src/modules/katabatic_routing_glacier.cpp b/src/modules/katabatic_routing_glacier.cpp new file mode 100644 index 00000000..35ea4215 --- /dev/null +++ b/src/modules/katabatic_routing_glacier.cpp @@ -0,0 +1,446 @@ +#include "katabatic_routing_glacier.hpp" +#include "Atmosphere.h" +#include "Glacier.hpp" +#include "PhysConst.h" +#include "melt_routing_glacier.hpp" + +// data and view constructors +REGISTER_MODULE_CPP(katabatic_routing_glacier); + +katabatic_routing_glacier::katabatic_routing_glacier(config_file cfg) + : module_base("katabatic_routing_glacier", parallel::data, cfg) +{ + depends("Pa"); + depends("t"); + depends("rh"); + depends("vapour_pressure_surface"); + // vapour pressure surface?? + depends("t_lapse_rate"); + depends("snowmelt_int"); + depends("swe"); + // TODO uncomment, only commented for CRHM coupling + //depends("iswr_net"); + //depends("iswr_subcanopy"); + //depends("ilwr_subcanopy"); + // TODO remove, only for CRHM coupling + // also should check if this exists in CHM + depends("T_rain"); + + // All submodule outputs + provides("glacier_water_equivalent"); + provides("total_depth"); + provides("firnmelt"); + provides("icemelt"); + provides("latent_heat"); + provides("sensible_heat"); + provides("snowmelt_delayed"); + provides("firnmelt_delayed"); + provides("icemelt_delayed"); + provides("total_delayed"); + provides("firn"); + provides("ice"); + +}; +katabatic_routing_glacier::~katabatic_routing_glacier() {}; +// module constructor +// module init +// module run + +void katabatic_routing_glacier::init(mesh& domain) +{ + for (size_t i = 0; i < domain->size_local_faces(); ++i) + { + auto face = domain->face(i); + auto& p_glacier = glacier.get_params(); + auto& p_routing = routing.get_params(); + auto& d = face->make_module_data(ID,face,global_param,&cfg,&p_glacier,&p_routing); + d.ice_emissivity = cfg.get("ice_emissivity"); + d.firn_emissivity = cfg.get("firn_emissivity"); + }; + + { + using namespace Glacier; + auto& p = glacier.get_params(); + + auto densify_version = cfg.get("densify_version","Linear"); + if (densify_version == "Linear") + p.densify_version = DensifyVersion::Linear; + else if (densify_version == "HerronLangway") + p.densify_version = DensifyVersion::HerronLangway; + else + { + p.densify_version = DensifyVersion::HerronLangway; + SPDLOG_DEBUG("Unknown densify version specified. Using HerronLangway as the default"); + } + + switch (p.densify_version) + { + case DensifyVersion::Linear: + p.small_increment = cfg.get("small_increment",25.0); + p.big_increment = cfg.get("big_increment",50.0); + break; + case DensifyVersion::HerronLangway: + break; + }; + + p.critical_density = cfg.get("critical_density",550.0); + p.firn_to_ice_density = cfg.get("firn_to_ice_density",830.0); + + p.thermal_factor = cfg.get("thermal_factor",0.95); + p.seconds_per_step = global_param->dt(); + } + { + auto& p = katabatic.get_params(); + p.prandtl = cfg.get("Pr",5.0); + p.k = cfg.get("k",4e-4); + p.k2 = cfg.get("k2",1.0); + p.seconds_per_step = global_param->dt(); + } + { + auto& p = routing.get_params(); + p.chain_length = cfg.get("routing_lag",10u); + p.seconds_per_step = global_param->dt(); + } +}; + +double katabatic_routing_glacier::rain_sun_energy(const mesh_elem& face) +{ + using namespace PhysConst; + constexpr auto M_PER_MM = 1 / 1000.0; + auto Qsun = (*face)["Qnsn_Var"_s]; + auto Qrain = Cw() * water_reference_density() * (*face)["rainfall_int"] * M_PER_MM + * ((*face)["T_rain"_s] - 0.0) / global_param->dt(); + + return Qsun + Qrain; +}; + +double katabatic_routing_glacier::rain_sun_energy(const mesh_elem& face,data& d) +{ + //auto iswr = (*face)["iswr_subcanopy"_s]; + //auto ilwr = (*face)["ilwr_subcanopy"_s]; + //const auto& firn = d.glacier_state.firn; + //const auto& ice = d.glacier_state.ice; + //auto emissivity = 0.0; + //if (d.swe().value > 0.0) + // return 0.0; + //else if (firn.water_equivalent().value) + // emissivity = d.firn_emissivity; + //else if (ice.water_equivalent().value) + // emissivity = d.ice_emissivity; + //else + // return 0.0; + // + //auto olwr = PhysConst::sbc() * emissivity + // * std::pow(d.air_temperature().value,4.0); + //auto oswr = (*face)["glacier_albedo"_s] * iswr; + //auto Qsun = (ilwr - olwr) + (iswr - oswr); + //auto Qrain = (*face)["Qrain"_s]; + auto Qsun = 0.0; + auto Qrain = 0.0; + return Qsun + Qrain; +}; + +void katabatic_routing_glacier::run(mesh_elem& face) +{ + auto& d = face->get_module_data(ID); + + { + // run every step, make sure the outputs are accumulating + // somewhere + katabatic_view data_k(d,face); + katabatic.execute(data_k); + const auto& cache = d.get_cache(); + // TODO this version of rain_run_energy is for CRHM coupling, do not include + d.total_energy += rain_sun_energy(face) + cache->latent_heat + + cache->sensible_heat; + // only run once per day + if (is_new_day()) + { + glacier_view data_g(d,face); + glacier.execute(data_g); + + // TODO reset latent/sensible heat here? + // could get_cache + // accumulate total_energy on public d member + // get d.melt_energy() return this value + // reset total_energy here + d.total_energy = 0.0; + } + + // run every step + routing_view data_r(d,face); + routing.execute(data_r); + //Implicitly output to face via deconstructors + //of *_view objects + } + + // *_view objects out of scope, safe to reset + d.reset_cache(); + +}; + +bool katabatic_routing_glacier::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; + +}; + +katabatic_routing_glacier::katabatic_view::~katabatic_view() +{ + auto& cache = d.get_cache(); + + (*face)["latent_heat"_s] = cache->latent_heat; + (*face)["sensible_heat"_s] = cache->sensible_heat; +}; + +katabatic_routing_glacier::routing_view::~routing_view() +{ + auto& cache = d.get_cache(); + + (*face)["snowmelt_delayed"_s] = cache->snowmelt_delayed; + (*face)["firnmelt_delayed"_s] = cache->firnmelt_delayed; + (*face)["icemelt_delayed"_s] = cache->icemelt_delayed; + (*face)["total_delayed"_s] = cache->snowmelt_delayed + cache->firnmelt_delayed + cache->icemelt_delayed; +}; + +katabatic_routing_glacier::glacier_view::~glacier_view() +{ + auto& cache = d.get_cache(); + (*face)["glacier_water_equivalent"] = d.glacier_state.total_water_equiv().value; + (*face)["total_depth"] = d.glacier_state.total_depth().value; + (*face)["firnmelt"] = cache->firnmelt; + (*face)["icemelt"] = cache->icemelt; + + (*face)["firn"] = d.glacier_state.firn.water_equivalent().value; + (*face)["ice"] = d.glacier_state.ice.water_equivalent().value; +}; + +static Glacier::Ice construct_ice_1(const Glacier::Params* p) +{ + using namespace Glacier; + auto init = Units::Milimetres{1e6}; + Ice ice(p,init); + + return ice; +}; + +static Glacier::LayeredFirn construct_firn_1(const Glacier::Params* p) +{ + using namespace Glacier; + struct simple_layer + { + double h; + double rho; + }; + + std::vector simple_layers = {{ + {1.0,450.0}, + {1.0,550.0}, + {1.0,650.0} + }}; + + std::vector layers; + + for (auto l : simple_layers) + { + layers.emplace_back(Layer::Height{l.h},Layer::Density{l.rho}); + }; + + LayeredFirn firn(p,layers); + return firn; +}; +katabatic_routing_glacier::data::data(const mesh_elem& face, std::shared_ptr global_param, const config_file* cfg, + const Glacier::Params* p_g, const GlacierRouting::Params* p_r) + : data_base(face,global_param,cfg), routing_state(p_r), glacier_state(construct_firn_1(p_g),construct_ice_1(p_g)) {}; + +Units::Kelvin katabatic_routing_glacier::data::glacier_temperature() +{ + if (!cfg) + { + std::string err = std::format("config file is null in katabatic_routing_glacier::data::glacier_temperature()"); + CHM_THROW_EXCEPTION(module_error,err); + }; + static const double T = 273.15;//cfg_.get("glacier_temp"_s,273.15); + + return Units::Kelvin{T}; +}; +Units::Kelvin katabatic_routing_glacier::data::air_temperature() +{ + update_value( [this]() -> auto& { return cache_->air_temperature; }, + [this]() { return (*face)["t"_s]; } ); + + return Units::Kelvin{cache_->air_temperature + 273.15}; +}; +Units::Pa katabatic_routing_glacier::data::air_pressure() +{ + update_value( [this]() -> auto& { return cache_->air_pressure; }, + [this]() { return (*face)["Pa"_s]; } ); + + return Units::Pa{cache_->air_pressure}; +}; +Units::Pa katabatic_routing_glacier::data::vapour_pressure() +{ + update_value( [this]() -> auto& { return cache_->vapour_pressure; }, + [this]() { + auto rh = (*face)["rh"_s]; + auto output = rh * Atmosphere::saturatedVapourPressure(air_temperature().value); + return output; } + ); + + return Units::Pa{cache_->vapour_pressure}; +}; +Units::Pa katabatic_routing_glacier::data::vapour_pressure_surface() +{ + update_value( [this]() -> auto& { return cache_->vapour_pressure_surface; }, + [this]() { return (*face)["vapour_pressure_surface"_s]; } ); + + return Units::Pa{cache_->vapour_pressure_surface}; +}; +Units::LapseRateSI katabatic_routing_glacier::data::lapse_rate() +{ + update_value( [this]() -> auto& { return cache_->lapse_rate; }, + [this]() { return (*face)["t_lapse_rate"_s]; } ); + + return Units::LapseRateSI{cache_->lapse_rate}; +}; +void katabatic_routing_glacier::data::latent_heat(const double v) +{ + set_output( [this]() -> auto& { return cache_->latent_heat; }, + v); +}; +void katabatic_routing_glacier::data::sensible_heat(const double v) +{ + set_output( [this]() -> auto& { return cache_->sensible_heat; }, + v); +}; + +const Units::Milimetres katabatic_routing_glacier::data::snowmelt() +{ + update_value( [this]() -> auto& { return cache_->snowmelt; }, + [this]() { return (*face)["snowmelt_int"_s]; } ); + + return Units::Milimetres{cache_->snowmelt}; +}; +const Units::Milimetres katabatic_routing_glacier::data::firnmelt() +{ + if (!cache_ || std::isnan(cache_->firnmelt)) + return Units::Milimetres{0.0}; + + return Units::Milimetres{cache_->firnmelt}; +}; +const Units::Milimetres katabatic_routing_glacier::data::icemelt() +{ + if (!cache_ || std::isnan(cache_->icemelt)) + return Units::Milimetres{0.0}; + + return Units::Milimetres{cache_->icemelt}; +}; +void katabatic_routing_glacier::data::snowmelt_delayed(double v) +{ + set_output( [this]() -> auto& { return cache_->snowmelt_delayed; }, + v); +}; +void katabatic_routing_glacier::data::firnmelt_delayed(double v) +{ + set_output( [this]() -> auto& { return cache_->firnmelt_delayed; }, + v); + +}; +void katabatic_routing_glacier::data::icemelt_delayed(double v) +{ + set_output( [this]() -> auto& { return cache_->icemelt_delayed; }, + v); +}; +void katabatic_routing_glacier::data::total_delayed(double v) +{ + set_output( [this]() -> auto& { return cache_->total_delayed; }, + v); +}; + +const Units::Milimetres katabatic_routing_glacier::data::swe() +{ + update_value( [this]() -> auto& { return cache_->swe; }, + [this]() { return (*face)["swe"_s]; } ); + + return Units::Milimetres{cache_->swe}; +}; +const Units::Watts_per_m2 katabatic_routing_glacier::data::melt_energy() +{ + auto iswr = (*face)["iswr_subcanopy"_s]; + auto ilwr = (*face)["ilwr_subcanopy"_s]; + const auto& firn = glacier_state.firn; + const auto& ice = glacier_state.ice; + auto emissivity = 0.0; + if (swe().value > 0.0) + return Units::Watts_per_m2{0.0}; + else if (firn.water_equivalent().value) + emissivity = firn_emissivity; + else if (ice.water_equivalent().value) + emissivity = ice_emissivity; + else + return Units::Watts_per_m2{0.0}; + + + auto olwr = PhysConst::sbc() * emissivity + * std::pow(air_temperature().value,4.0); + auto oswr = (*face)["glacier_albedo"_s] * iswr; + auto Qsun = (ilwr - olwr) + (iswr - oswr); + auto Qrain = (*face)["Qrain"_s]; + // Means that katabatic_melt_energy didn't run + if (!cache_) + return Units::Watts_per_m2{Qsun + Qrain}; + + auto Q_sensible = 0.0; + auto Q_latent = 0.0; + + if ( !std::isnan(cache_->latent_heat) ) + Q_latent = cache_->latent_heat; + + if ( !std::isnan(cache_->sensible_heat) ) + Q_sensible = cache_->sensible_heat; + + return Units::Watts_per_m2{Qsun + Qrain + Q_sensible + Q_latent}; +}; +bool katabatic_routing_glacier::data::update_now() +{ + static const size_t day = cfg->get("change_day"_s,300u); + + auto p = global_param->posix_time(); + auto d = p.date(); + + auto jd = d.day_of_year(); + std::cout << jd << std::endl; + + return day == jd; +}; + +void katabatic_routing_glacier::data::glacier_water_equivalent(const double out) +{ + set_output( [this]() -> auto& { return cache_->glacier_water_equivalent; }, + out); +}; +void katabatic_routing_glacier::data::total_depth(const double out) +{ + set_output( [this]() -> auto& { return cache_->total_depth; }, + out); +}; +void katabatic_routing_glacier::data::firn_melt(const double out) +{ + set_output( [this]() -> auto& { return cache_->firnmelt; }, + out); +}; +void katabatic_routing_glacier::data::ice_melt(const double out) +{ + set_output( [this]() -> auto& { return cache_->icemelt; }, + out); +}; + + diff --git a/src/modules/katabatic_routing_glacier.hpp b/src/modules/katabatic_routing_glacier.hpp new file mode 100644 index 00000000..af16a2f8 --- /dev/null +++ b/src/modules/katabatic_routing_glacier.hpp @@ -0,0 +1,214 @@ +// +// 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 "triangulation.hpp" +#include "module_base.hpp" +#include "Glacier.hpp" +#include "katabatic_melt_energy.hpp" +#include "melt_routing_glacier.hpp" +#include "data_base.hpp" +#include "PhysConst.h" + +/** + * \ingroup modules infil soils exp + * @{ + * \class Gray_inf + * + * + * 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 katabatic_routing_glacier : public module_base +{ +REGISTER_MODULE_HPP(katabatic_routing_glacier) +public: + katabatic_routing_glacier(config_file cfg); + + ~katabatic_routing_glacier(); + + void run(mesh_elem &face) override; + void init(mesh& domain) override; + + struct Cache : public cache_base + { + Cache() {}; + // All submodule inputs + double swe = cache_base::default_value(); + double melt_energy = cache_base::default_value(); + double glacier_temperature = cache_base::default_value(); + double update_now = cache_base::default_value(); + double air_temperature = cache_base::default_value(); + double air_pressure = cache_base::default_value(); + double vapour_pressure = cache_base::default_value(); + double vapour_pressure_surface = cache_base::default_value(); + double lapse_rate = cache_base::default_value(); + double snowmelt = cache_base::default_value(); + + // All submodule outputs + double glacier_water_equivalent = 0.0; + double total_depth = 0.0; + double firnmelt = 0.0; + double icemelt = 0.0; + double latent_heat = 0.0; + double sensible_heat = 0.0; + double snowmelt_delayed = 0.0; + double firnmelt_delayed = 0.0; + double icemelt_delayed = 0.0; + double total_delayed = 0.0; + }; + + class data : public face_info, public data_base + { + public: + data(const mesh_elem& face, const std::shared_ptr param, const config_file* cfg, + const Glacier::Params* g_p, const GlacierRouting::Params* r_p); + + GlacierRouting::State routing_state; + Glacier::State glacier_state; + + Units::Kelvin glacier_temperature(); + Units::Kelvin air_temperature(); + Units::Pa air_pressure(); + Units::Pa vapour_pressure(); + Units::Pa vapour_pressure_surface(); + Units::LapseRateSI lapse_rate(); + void latent_heat(const double v); + void sensible_heat(const double v); + + const Units::Milimetres snowmelt(); + const Units::Milimetres firnmelt(); + const Units::Milimetres icemelt(); + void snowmelt_delayed(double v); + void firnmelt_delayed(double v); + void icemelt_delayed(double v); + void total_delayed(double v); + + const Units::Milimetres swe(); + const Units::Watts_per_m2 melt_energy(); + bool update_now(); + + void glacier_water_equivalent(const double out); + void total_depth(const double out); + void firn_melt(const double out); + void ice_melt(const double out); + + double firn_emissivity = 0.0; + double ice_emissivity = 0.0; + double total_energy = 0.0; + }; + + // Adapter that satisfies KatabaticData concept + class katabatic_view + { + data& d; + mesh_elem& face; + public: + explicit katabatic_view(data& d,mesh_elem& face) + : d(d), face(face){}; + ~katabatic_view(); + Units::Kelvin glacier_temperature() { return d.glacier_temperature(); } + Units::Kelvin air_temperature() { return d.air_temperature(); } + Units::Pa air_pressure() { return d.air_pressure(); } + Units::Pa vapour_pressure() { return d.vapour_pressure(); } + Units::Pa vapour_pressure_surface() { return d.vapour_pressure_surface(); } + Units::LapseRateSI lapse_rate() { return d.lapse_rate(); } + void latent_heat(const double v) { d.latent_heat(v); } + void sensible_heat(const double v) { d.sensible_heat(v); } + }; + + + // Adapter that satisfies GlacierRoutingData concept + class routing_view { + data& d; + mesh_elem& face; + public: + explicit routing_view(data& d,mesh_elem& face) : d(d),face(face) {} + ~routing_view(); + GlacierRouting::State& get_state() { return d.routing_state; } // no name collision! + const Units::Milimetres snowmelt() { return d.snowmelt(); } + const Units::Milimetres firnmelt() { return d.firnmelt(); } + const Units::Milimetres icemelt() { return d.icemelt(); } + void snowmelt_delayed(const double v) { d.snowmelt_delayed(v); } + void firnmelt_delayed(const double v) { d.firnmelt_delayed(v); } + void icemelt_delayed(const double v) { d.icemelt_delayed(v); } + void total_delayed(const double v) { d.total_delayed(v); } + }; + + // Similarly for Glacier::Model's concept + class glacier_view { + data& d; + mesh_elem& face; + public: + explicit glacier_view(data& d,mesh_elem& face) + : d(d), face(face) {} + ~glacier_view(); + Glacier::State& get_state() { return d.glacier_state; } + const Units::Milimetres swe() { return d.swe(); } + const Units::Watts_per_m2 melt_energy() { return Units::Watts_per_m2{d.total_energy}; } + const Units::Celsius glacier_temperature() { return d.glacier_temperature(); } + bool update_now() { return d.update_now(); } + + void glacier_water_equivalent(const double out) { d.glacier_water_equivalent(out); } + void total_depth(const double out) { d.total_depth(out); } + void firn_melt(const double out) { d.firn_melt(out); } + void ice_melt(const double out) { d.ice_melt(out); } + }; + +private: + katabatic_melt_energy::Model katabatic; + GlacierRouting::Model routing; + Glacier::Model glacier; + + bool is_new_day(); + double rain_sun_energy(const mesh_elem& face); + double rain_sun_energy(const mesh_elem& face,data& d); +}; diff --git a/src/modules/point_mode.cpp b/src/modules/point_mode.cpp index c740cb1a..c7acce8d 100644 --- a/src/modules/point_mode.cpp +++ b/src/modules/point_mode.cpp @@ -22,24 +22,117 @@ // #include "point_mode.hpp" +#include "Atmosphere.h" REGISTER_MODULE_CPP(point_mode); point_mode::point_mode(config_file cfg) : module_base("point_mode", parallel::data, cfg) { - t = cfg.get("provide.t",true); - rh = cfg.get("provide.rh",true); - U_R = cfg.get("provide.U_R",true); - U_2m_above_srf = cfg.get("provide.U_2m_above_srf",true); - p = cfg.get("provide.p",true); - ilwr = cfg.get("provide.ilwr",true); - iswr = cfg.get("provide.iswr",true); - vw_dir = cfg.get("provide.vw_dir",true); - iswr_diffuse = cfg.get("provide.iswr_diffuse",false); - iswr_direct = cfg.get("provide.iswr_direct",false); - T_g = cfg.get("provide.T_g",false); + t = cfg.get("provide.t",false); + rh = cfg.get("provide.rh",false); + U_R = cfg.get("provide.U_R",false); + U_2m_above_srf = cfg.get("provide.U_2m_above_srf",false); + p = cfg.get("provide.p",false); + ilwr = cfg.get("provide.ilwr",false); + iswr = cfg.get("provide.iswr",false); + vw_dir = cfg.get("provide.vw_dir",false); + iswr_diffuse = cfg.get("provide.iswr_diffuse",false); + iswr_direct = cfg.get("provide.iswr_direct",false); + T_g = cfg.get("provide.T_g",false); + svf = cfg.get("provide.svf",false); + + unit_test_new_modules = cfg.get("provide.unit_test_Donovan",false); + unit_test_glacier_module = cfg.get("provide.unit_test_glacier_module",false); + if (unit_test_new_modules && unit_test_glacier_module) + CHM_THROW_EXCEPTION(module_error,"Point_mode: conflicting booleans are both true"); + + if(unit_test_new_modules) + { + // Infil_All + depends_from_met("SWE"); + provides("swe"); + + depends_from_met("snowmeltD"); + provides("snowmelt_int"); + + depends_from_met("net_rain"); + provides("rainfall_int"); + + depends_from_met("fallstat_V"); + provides("soil_storage_at_freeze"); + + // t is included in another section + + // Evapotranspiration_All + depends_from_met("Rn"); + provides("netall"); + + depends_from_met("Pa"); + provides("P_atm"); + + depends_from_met("hru_ea"); + provides("ea"); + + // rh, t, U_2m_above_srf included already here + // soil_storage comes from the soil module + + // soil_module + depends_from_met("rho"); + provides("snow_density"); + // All depends/provides come from other modules + // swe, ET, inf, runoff + // + depends_from_met("hru_t"); + provides("t"); + + depends_from_met("hru_u"); + provides("U_2m_above_srf"); + + depends_from_met("hru_tsf"); + provides("surface_temperature"); + } + else if(unit_test_glacier_module) + { + depends_from_met("SWE"); + provides("swe"); + + // TODO Convert to per-interval before setting to face + depends_from_met("snowmeltD"); + provides("snowmelt_int"); + + // TODO Adjust katabatic module to accept thius and convert to Qrain + depends_from_met("net_rain"); + provides("rainfall_int"); + + depends_from_met("hru_t"); + provides("t"); + + depends_from_met("Pa"); + provides("Pa"); + + // TODO CRHM uses this as a system wide parameter, provide it here but dont look at met + // Look at (*face)-> instead + provides("t_lapse_rate"); + + // TODO look at module, adjust functions to accept rh, that means undoing the conversion here + // Careful with units... + depends_from_met("hru_ea"); + provides("rh"); + + // Constant in crhm, provdided here + // TODO, adjust module to compute this itself + provides("vapour_pressure_surface"); + + depends_from_met("T_rain"); + provides("T_rain"); + + + depends_from_met("Qnsn_Var"); + provides("Qnsn_Var"); + + } if(t) { depends_from_met("t"); @@ -110,11 +203,93 @@ point_mode::~point_mode() } +void point_mode::init(mesh &domain) +{ + size_t i=0; + auto face = domain->face(i); + elevation = face->parameter("elevation"_s); +}; void point_mode::run(mesh_elem &face) { // at this point, if the user has provided more than 1 station, they've been stripped out. // we can safetly take the 1st (and only) station. + if(unit_test_new_modules) + { + double swe = (*face->nearest_station())["SWE"_s]; + (*face)["swe"_s] = swe; + + double snowmelt_int = (*face->nearest_station())["snowmeltD"_s]; + (*face)["snowmelt_int"_s] = snowmelt_int/24; + + double rainfall_int = (*face->nearest_station())["net_rain"_s]; + (*face)["rainfall_int"_s] = rainfall_int; + + double soil_storage_at_freeze = (*face->nearest_station())["fallstat_V"_s]; + (*face)["soil_storage_at_freeze"_s] = soil_storage_at_freeze; + + double netall = (*face->nearest_station())["Rn"_s]; + (*face)["netall"_s] = netall; + + double P_atm = 101.3*pow( (293.0 - 0.0065*elevation)/293.0,5.26); + (*face)["P_atm"_s] = P_atm; + + double ea = (*face->nearest_station())["hru_ea"_s]; + (*face)["ea"_s] = ea; + + double rho = (*face->nearest_station())["rho"_s]; + (*face)["snow_density"_s] = rho; + + double t = (*face->nearest_station())["hru_t"_s]; + (*face)["t"_s] = t; + + double u = (*face->nearest_station())["hru_u"_s]; + (*face)["U_2m_above_srf"_s] = u; + + double temp = (*face->nearest_station())["hru_tsf"_s]; + (*face)["surface_temperature"_s] = temp; + } + else if (unit_test_glacier_module) + { + double SWE = (*face->nearest_station())["SWE"_s]; + (*face)["swe"_s] = SWE; + + // TODO Convert to per-interval before setting to face + double snowmeltD = (*face->nearest_station())["snowmeltD"_s]; + (*face)["snowmelt_int"_s] = snowmeltD * global_param->dt() / (86400.0); + + // TODO Adjust katabatic module to accept thius and convert to Qrain + double net_rain = (*face->nearest_station())["net_rain"_s]; + (*face)["rainfall_int"_s] = net_rain; + + // TODO use equation above to compute pressure, not from provides + double Pa = (*face->nearest_station())["Pa"_s]; + (*face)["Pa"_s] = Pa * 1000.0; //in CRHM Pa is kPa + + // EXPECTS KELVIN + double t = (*face->nearest_station())["hru_t"_s]; + (*face)["t"_s] = t; + + // TODO CRHM uses this as a system wide parameter, provide it here but dont look at met + // Look at (*face)-> instead + + (*face)["t_lapse_rate"_s] = cfg.get("t_lapse_rate"); + + // TODO look at module, adjust functions to accept rh, that means undoing the conversion here + // Careful with units... + double hru_ea = (*face->nearest_station())["hru_ea"_s]; + double es = Atmosphere::saturatedVapourPressure(t + 273.15); + (*face)["rh"_s] = hru_ea / ( es ) * 1000.0; + + double vapour_pressure_surface = 611.3; + (*face)["vapour_pressure_surface"_s] = vapour_pressure_surface; + + double T_rain = (*face->nearest_station())["T_rain"_s]; + (*face)["T_rain"_s] = T_rain; + + double Q = (*face->nearest_station())["Qnsn_Var"_s]; + (*face)["Qnsn_Var"_s] = Q; + } if(t) { double st =(*face->nearest_station())["t"_s]; @@ -180,7 +355,8 @@ void point_mode::run(mesh_elem &face) double T_g =(*face->nearest_station())["T_g"_s]; (*face)["T_g"_s]= T_g; } - - face->parameter("svf") = cfg.get("override.svf", face->parameter("svf"_s)); - + if(svf) + { + face->parameter("svf") = cfg.get("override.svf", face->parameter("svf"_s)); + } } diff --git a/src/modules/point_mode.hpp b/src/modules/point_mode.hpp index 5442aa5f..36f582a2 100644 --- a/src/modules/point_mode.hpp +++ b/src/modules/point_mode.hpp @@ -113,7 +113,7 @@ REGISTER_MODULE_HPP(point_mode); ~point_mode(); virtual void run(mesh_elem &face); - + virtual void init(mesh &domain); bool t ; bool rh ; @@ -126,6 +126,9 @@ REGISTER_MODULE_HPP(point_mode); bool iswr_direct ; bool U_2m_above_srf ; bool T_g; - - + bool svf; + bool unit_test_new_modules; + bool unit_test_glacier_module; + double elevation = 0.0; + double soil_storage_at_freeze = 0.0; }; diff --git a/src/modules/submodules/Glacier.cpp b/src/modules/submodules/Glacier.cpp new file mode 100644 index 00000000..bf4c7b46 --- /dev/null +++ b/src/modules/submodules/Glacier.cpp @@ -0,0 +1,428 @@ +#include "Glacier.hpp" +#include "PhysConst.h" +#include "math/Bisection.hpp" +#include +#include +#include + +namespace Glacier +{ + constexpr auto MM_PER_M = 1000.0; + constexpr auto MegaGram_per_KiloGram = 1 / 1000.0; + namespace Densification + { + static constexpr auto ice_density = PhysConst::rho_ice(); + + static const double k_1(const double T) + { + // Empirical constants + constexpr auto a = 11.0 * MegaGram_per_KiloGram; + constexpr auto b = 10160.0; + return a * std::exp(-b/(PhysConst::R() * T)); + + }; + static const double Z(const double h, const double rho, const double T) + { + auto k = k_1(T); + return std::exp(ice_density * k * h) * rho / ( ice_density - rho); + }; + static const double Z(const double h, const double rho, const double T, const double A, + const double rho_c) + { + constexpr double a = 575.0, b = 21400.0; + auto k_2 = a * std::exp(-b/(PhysConst::R() * T)) * MegaGram_per_KiloGram; + + auto critical_depth = 1/(ice_density * k_1(T)) * + std::log( rho_c * ( ice_density - rho ) / (rho * (ice_density - rho_c))); + + return std::exp(ice_density * k_2 * (h - critical_depth) / std::sqrt(A)) * rho_c / (ice_density - rho_c); + }; + static const double HerronLangwayFormula(const double Z) + { + return ice_density * Z / (1 + Z); + }; + const double HerronLangway(const double depth, const double density, + const double glacier_temperature,const double accumulation_rate,const double critical_density) + { + if (density >= ice_density) [[unlikely]] + { + std::string err = std::format("Herron and Langway Densification formula cannot handle a density at or above {} kg/m^3. Received {} kg/m^3",ice_density,density); + throw std::runtime_error(err); + } + + double z; + + if (density <= critical_density) + z = Z(depth,density,glacier_temperature); + else + z = Z(depth,density,glacier_temperature,accumulation_rate,critical_density); + + return HerronLangwayFormula(z); + }; + + const double Linear(const double rho, const double delta_rho) + { + return rho + delta_rho; + }; + + }; // Densification namespace + + static constexpr double swe_to_firn(const double h) { + return 450.0 - 204.7 / h * (1 - std::exp(-h / 0.673)); + }; + + Layer::Layer(const Units::Milimetres water_eq) + { + using namespace math; + auto test_function = [=,WE = water_eq](const double h) { + auto height = std::max(h,0.6); + auto rho = swe_to_firn(height); + return rho * height * MM_PER_M / PhysConst::water_reference_density() - WE.value; + }; + constexpr double low_guess = 0.6; + constexpr double high_guess = 8.0; + constexpr double tolerance = 1e-15; + + auto result = bisection(test_function, + low_guess,high_guess, + tolerance); + + switch (result.root) { + case Root::Found: + _height.value = result.value; + _density.value = swe_to_firn(_height.value); + break; + case Root::DidNotConverge: + { + std::string err = std::format("Layer Bisection did not converge, with guess bounds ({},{}) and water equivalent {} mm",low_guess,high_guess,water_eq.value); + throw std::runtime_error(err); + break; + } + case Root::NoRootGuaranteed: + { + /* + * swe_to_firn is only well defined for snow depths above 0.6 m + * Therefore, for WE < rho(0.6) * 0.6, it means that the SWE value is + * unsuited for this equation. Instead, take the minimal density value + * and recompute height + */ + if ( swe_to_firn(0.6) * 0.6 - water_eq.value > 0.0 ) + { + _density.value = swe_to_firn(0.6); + update_height(water_eq.value); + } + else + { + std::string err = std::format("Layer Bisection might never converge, with guess bounds ({},{}) and water equivalent {} mm",low_guess,high_guess,water_eq.value); + throw std::runtime_error(err); + } + break; + } + } + + }; + + Layer::Layer(const Height h, const Density rho) + { + _height = h; + _density = rho; + }; + + const Layer::Height Layer::height() const { return _height;}; + const Layer::Density Layer::density() const {return _density;}; + + void Layer::densifyLinear(const Density small, const Density big, const Density critical) + { + Layer::Density increment{0.0}; + auto WE = water_equivalent(); + + if (_density.value <= critical.value) + increment += small; + else + increment += big; + + _density.value = Densification::Linear(_density.value,increment.value); + update_height(WE.value); + }; + + void Layer::densifyHerronLangway(const Units::Kelvin T, double depth, + const double accumulation_rate, const Density critical_density) + { + depth += _height.value; + auto WE = water_equivalent(); + + auto new_value = Densification::HerronLangway(depth,_density.value,T.value,accumulation_rate,critical_density.value); + if (new_value > _density.value) + { + _density.value = new_value; + update_height(WE.value); + } + }; + + void Layer::update_height(const double old) + { + + _height.value = old * PhysConst::water_reference_density() / (_density.value * MM_PER_M); + }; + + void Layer::remove(const Units::Milimetres amount) + { + auto current = water_equivalent(); + if (amount >= current) + throw std::invalid_argument("Layer::remove amount must not be in excess of or equal to current"); + + // Not the same as update_height, this is a flat height reduction without changing density + current -= amount; + _height.value = current.value * PhysConst::water_reference_density() / + ( _density.value * MM_PER_M ); + }; + + const Units::Milimetres Layer::water_equivalent() const + { + return Units::Milimetres{_height.value * _density.value / PhysConst::water_reference_density() * MM_PER_M}; + }; + + LayeredFirn::LayeredFirn(const Params* _p) : p(_p) {}; + + LayeredFirn::LayeredFirn(const Params* _p,const std::vector& L) : p(_p) { + layers.assign(L.begin(),L.end()); + }; + + LayeredFirn::LayeredFirn(const Params* _p, const std::deque& L) : p(_p), layers(L) {}; + + const Units::Milimetres LayeredFirn::water_equivalent() const + { + if (layers.empty()) + return Units::Milimetres{0.0}; + + auto sum = std::accumulate(layers.begin(),layers.end(),0.0, + [](double sum, const Layer& layer) { + return sum + layer.water_equivalent().value; + }); + + return Units::Milimetres{sum}; + }; + + void LayeredFirn::accumulate(const Units::Milimetres WE) + { + Layer new_layer(WE); + layers.push_back(new_layer); + }; + + const MeltInfo LayeredFirn::melt(Units::Milimetres input) + { + Units::Milimetres melt{0.0}; + Units::Milimetres init_layer_WE{0.0}; + while (input.value > 0.0 && !layers.empty()) + { + Layer& current = layers.back(); + + if ( input.value < current.water_equivalent().value ) + { + init_layer_WE = current.water_equivalent(); + current.remove(input); + melt += input; + input.value = 0.0; + if (current.water_equivalent().value < 1e-12) + layers.pop_back(); + } + else + { + const auto value = current.water_equivalent(); + melt += value; + input -= value; + layers.pop_back(); + }; + }; + + if (input.value > 0.0 && input.value < 1e-12) + { + input.value = 0.0; + layers.pop_back(); + }; + + return MeltInfo { + melt, + input + }; + + }; + + const LayeredFirn::LayerContainer& LayeredFirn::get_layers() const + { + return layers; + }; + + std::optional LayeredFirn::convert_to_ice() + { + std::optional result; + + if (layers.empty()) + return result; + + if (layers.front().density().value > p->firn_to_ice_density) + { + result.emplace(layers.front().water_equivalent()); + layers.pop_front(); + } + return result; + }; + + Ice::Ice(const Params* _p) : p(_p) {}; + + Ice::Ice(const Params* _p, const Units::Milimetres WE) : p(_p), _water_equivalent(WE) {}; + + void Ice::accumulate(const Units::Milimetres WE) + { + _water_equivalent += WE; + }; + + const MeltInfo Ice::melt(const Units::Milimetres input) + { + Units::Milimetres local_input = input; + auto melted = Units::Milimetres{0.0}; + + if (input > _water_equivalent) + { + melted += _water_equivalent; + _water_equivalent.value = 0.0; + } + else + { + melted += input; + _water_equivalent -= input; + } + + Units::Milimetres remaining_energy = local_input; + remaining_energy -= melted; + remaining_energy.value = std::max(remaining_energy.value,0.0); + + return MeltInfo{melted, + remaining_energy}; + }; + + const Units::Milimetres Ice::water_equivalent() const + { + return _water_equivalent; + } + + namespace Melt + { + Units::Milimetres convert_to_mass(const Units::Watts_per_m2 energy,const Params& p) + { + return Units::Milimetres{energy.value * p.seconds_per_step * MM_PER_M / (PhysConst::Lf() * PhysConst::water_reference_density() * p.thermal_factor)}; + }; + + MeltScenario get_scenario(const State& s, + Units::Milimetres melt_energy,const Units::Milimetres swe) + { + if (swe.value > 0.0 || melt_energy.value == 0.0) + return MeltScenario::NoMelt; + + if (s.firn.water_equivalent().value > 0.0 + && s.ice.water_equivalent().value == 0.0) + return MeltScenario::FirnMelt; + + if (s.firn.water_equivalent().value == 0.0 + && s.ice.water_equivalent().value > 0.0) + return MeltScenario::IceMelt; + + if (s.firn.water_equivalent().value > 0.0 + && s.ice.water_equivalent().value > 0.0) + { + if ( melt_energy.value > s.firn.water_equivalent().value) + return MeltScenario::FirnAndIceMelt; + else + return MeltScenario::FirnMelt; + } + + if (s.total_water_equiv().value == 0.0) + return MeltScenario::NoMelt; + + std::string err = + std::format("Unhandled melt scenario in glacier submodule for SWE of {} mm, firn of {} mm, and ice of {} mm",swe.value,s.firn.water_equivalent().value,s.ice.water_equivalent().value); + throw std::logic_error(err); + }; + + std::optional compute_melt(MeltScenario scenario,State& s,const Units::Milimetres melt_mass_max) + { + std::optional info; + switch(scenario) + { + case MeltScenario::NoMelt: + // No melt, so info remains uninitialized + break; + case MeltScenario::FirnMelt: + info.emplace(s.firn.melt(melt_mass_max)); + break; + case MeltScenario::FirnAndIceMelt: { + info.emplace(s.firn.melt(melt_mass_max)); + MeltInfo ice_info = s.ice.melt(info->remaining_energy); + info->melt += ice_info.melt; + info->remaining_energy += ice_info.remaining_energy; + } + break; + case MeltScenario::IceMelt: + info.emplace(s.ice.melt(melt_mass_max)); + + } + return info; + + }; + + + }; + + namespace Updates + { + void swe_to_firn(const Params& p, LayeredFirn& firn, Units::Milimetres& swe) + { + if (swe.value == 0.0) + return; + + firn.accumulate(swe); + + swe.value = 0.0; + + }; + + void firn_to_ice(const Params& p, LayeredFirn& firn, Ice& ice) + { + auto new_ice = firn.convert_to_ice(); + + if (new_ice) + ice.accumulate(*new_ice); + + }; + } + + State::State(const Params* p) : firn(p), ice(p) {}; + + State::State(const Params* p, const Ice& i) : firn(p), ice(i) {}; + + State::State(const LayeredFirn& f,const Ice& i) : firn(f), ice(i) {}; + + const Units::Milimetres State::total_water_equiv() const + { + auto result = ice.water_equivalent(); + result += firn.water_equivalent(); + + return result; + }; + + const Units::Metres State::total_depth() const + { + using namespace PhysConst; + Units::Metres sum{water_reference_density() / rho_ice() * + ice.water_equivalent().value / MM_PER_M}; + + auto& layers = firn.get_layers(); + + for (auto layer : layers) + { + sum += layer.height(); + } + + return sum; + }; +}; diff --git a/src/modules/submodules/Glacier.hpp b/src/modules/submodules/Glacier.hpp new file mode 100644 index 00000000..1625377a --- /dev/null +++ b/src/modules/submodules/Glacier.hpp @@ -0,0 +1,230 @@ +#pragma once + +#include "base_step.hpp" +#include +#include +#include +#include "PhysConst.h" + +namespace Units = PhysConst::units; +namespace Glacier +{ + enum class DensifyVersion { + HerronLangway, + Linear + }; + + struct Params + { + const double water_density = + PhysConst::water_reference_density(); + const double heat_capacity_air = + PhysConst::Cp(); + const double latent_heat_vapour = + PhysConst::Lv(); + const double gas_constant_dry = + PhysConst::RgasDry(); + const double gas_constant_vapour = + PhysConst::RgasVapour(); + DensifyVersion densify_version = DensifyVersion::HerronLangway; + double small_increment{25.0}; + double big_increment{50.0}; + double critical_density{550.0}; + double firn_to_ice_density{830.0}; + double thermal_factor{0.95}; + size_t seconds_per_step{3600}; + }; + + namespace Densification + { + // Low density + const double HerronLangway(const double depth, const double density, + const double glacier_temperature); + // Above density + const double HerronLangway(const double depth, const double density, + const double glacier_temperature, const double accumulation_rate, + const double critical_density = 550.0); + + const double Linear(const double density,const double density_increment); + }; + + + class Layer + { + public: + struct Height : public Units::Metres {}; + struct Density : public Units::DensitySI {}; + + Layer(const Units::Milimetres water_equivalent); + Layer(const Height,const Density); + + const Units::Milimetres water_equivalent() const; + void densifyLinear(const Density small, const Density big,const Density critical); + void densifyHerronLangway(const Units::Kelvin, double depth, + const double accumulation_rate, const Density critical_density); + void remove(const Units::Milimetres); + + const Height height() const; + const Density density() const; + + private: + Height _height; + Density _density; + + void update_height(const double old); + + }; + + struct MeltInfo + { + Units::Milimetres melt; + Units::Milimetres remaining_energy; + }; + + class LayeredFirn + { + public: + using LayerContainer = std::deque; + + private: + const Params* const p; + LayerContainer layers; + + public: + explicit LayeredFirn(const Params*); + LayeredFirn(const Params*,const std::vector&); + LayeredFirn(const Params*,const LayerContainer&); + std::optional convert_to_ice(); + + void accumulate(const Units::Milimetres); + + const MeltInfo melt(const Units::Milimetres); + const Units::Milimetres water_equivalent() const; + const LayerContainer& get_layers() const; + }; + + class Ice + { + const Params* const p; + Units::Milimetres _water_equivalent{0.0}; + public: + explicit Ice(const Params*); + Ice(const Params*, const Units::Milimetres); + + void accumulate(const Units::Milimetres); + const MeltInfo melt(const Units::Milimetres); + const Units::Milimetres water_equivalent() const; + }; + + struct State + { + LayeredFirn firn; + Ice ice; + + State() = delete; + State(const Params* p); + State(const Params* p, const Ice& i); + State(const Params* p, const LayeredFirn& f); + State(const LayeredFirn& f, const Ice& i); + const Units::Milimetres total_water_equiv() const; + const Units::Metres total_depth() const; + }; + + namespace Melt + { + enum class MeltScenario { + NoMelt, FirnMelt, FirnAndIceMelt, IceMelt}; + + Units::Milimetres convert_to_mass(const Units::Watts_per_m2,const Params&); + + MeltScenario get_scenario(const State&,const Units::Milimetres melt_energy,const Units::Milimetres swe); + + std::optional compute_melt(MeltScenario,State&,const Units::Milimetres); + }; + + namespace Updates + { + void swe_to_firn(const Params&, LayeredFirn&, Units::Milimetres& swe); + void firn_to_ice(const Params&, LayeredFirn&, Ice&); + }; + + template + concept GlacierData = requires(T& t) + { + { t.get_state() } -> std::same_as; + { t.swe() } -> std::same_as; + { t.melt_energy() } -> std::same_as; + { t.glacier_temperature() } -> std::same_as; + { t.update_now() } -> std::same_as; + + { t.glacier_water_equivalent(std::declval()) } -> std::same_as; + { t.total_depth(std::declval()) } -> std::same_as; + { t.firn_melt(std::declval()) } -> std::same_as; + { t.ice_melt(std::declval()) } -> std::same_as; + }; + + template + class Model : public base_step,data> + { + public: + void execute_impl(data& d) const + { + auto& s = d.get_state(); + auto SWE = d.swe(); + Units::Milimetres melt_energy_mm_per_dt = Melt::convert_to_mass(d.melt_energy(),p); + + if (melt_energy_mm_per_dt.value > 0.0) + { + using namespace Melt; + auto scenario = get_scenario(s, + melt_energy_mm_per_dt,SWE); + + std::optional state_copy; + + if (scenario == MeltScenario::FirnAndIceMelt) + state_copy.emplace(s); + + auto info = compute_melt(scenario,s,melt_energy_mm_per_dt); + + switch (scenario) { + case MeltScenario::NoMelt: + d.firn_melt(0.0); + d.ice_melt(0.0); + break; + case MeltScenario::FirnMelt: + d.firn_melt(info->melt.value); + d.ice_melt(0.0); + break; + case MeltScenario::IceMelt: + d.ice_melt(info->melt.value); + d.firn_melt(0.0); + break; + case MeltScenario::FirnAndIceMelt: + d.firn_melt(state_copy->firn.water_equivalent().value); + + d.ice_melt(state_copy->ice.water_equivalent().value + - s.ice.water_equivalent().value); + break; + } + } + + if (d.update_now()) + { + using namespace Updates; + + swe_to_firn(p,s.firn,SWE); + + firn_to_ice(p,s.firn,s.ice); + } + + d.glacier_water_equivalent(s.total_water_equiv().value); + d.total_depth(s.total_depth().value); + }; + Params& get_params() + { + return p; + }; + private: + Params p; + }; +}; diff --git a/src/modules/submodules/base_step.hpp b/src/modules/submodules/base_step.hpp new file mode 100644 index 00000000..bb632dd7 --- /dev/null +++ b/src/modules/submodules/base_step.hpp @@ -0,0 +1,68 @@ +// +// 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 + +/** + * @brief Base class for submodules using the Curiously Recurring Template Pattern (CRTP). + * + * This class enforces that derived classes implement an `execute_impl(Data&) const` method, + * and provides a public `execute(Data&) const` method that forwards to the derived implementation. + * + * The CRTP pattern allows compile-time polymorphism, enabling static dispatch and + * avoiding the overhead of virtual functions. This is useful for performance-critical + * code or when templates are preferred over inheritance. + * + * Usage example: + * @code + * struct MyStep : public base_step { + * void execute_impl(MyData& d) { + * // Step-specific logic here + * } + * }; + * MyStep step; + * MyData data; + * step.execute(data); // Calls MyStep::execute_impl(data) + * @endcode + */ + +template +concept GuaranteeImplementExecute = requires(const T& t, Data& d) { + { t.execute_impl(d) } -> std::same_as; +}; + +template +class base_step { +protected: + base_step() {}; +public: + void execute(Data& d) const { + + static_assert(GuaranteeImplementExecute, + "Derived class must implement: void execute_impl(Data&) const"); + + static_cast(this)->execute_impl(d); + } + ~base_step() {}; +}; diff --git a/src/modules/submodules/data_base.hpp b/src/modules/submodules/data_base.hpp new file mode 100644 index 00000000..eeae2568 --- /dev/null +++ b/src/modules/submodules/data_base.hpp @@ -0,0 +1,224 @@ +// +// 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 +// . +// + +// +// Created by Donovan Allum 2025 +// + +#pragma once + +#include "global.hpp" +#include "triangulation.hpp" +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/** + * @brief Base class for data management with lazy caching and output handling. + * + * This class provides a framework for managing data with optional caching, + * lazy evaluation, and output assignment. It uses CRTP-like patterns with + * concepts to enforce interface contracts at compile-time. + * + * Key features: + * - Automatic cache initialization and staleness checking + * - Lazy evaluation of values with cache-aware updates + * - Type-safe output assignment with runtime checks + * - Compile-time concept enforcement for derived cache types and value/output providers + * + * The class is designed for computational scenarios where: + * - Data may be expensive to compute and should be cached + * - Cache validity depends on external timestep counters + * - Values need to be fetched lazily when accessed + * - Outputs must be accumulated safely + * + * Usage example: + * @code + * struct MyCache : public cache_base { + * double temperature = default_value(); + * int iteration_count = default_value(); + * }; + * + * class MyData : public data_base { + * public: + * MyData(const mesh_elem& face_in, const std::shared_ptr param, + * const pt::ptree& cfg) : data_base(face_in, param, cfg) {} + * + * void compute_temperature() { + * update_value( + * [this]() -> auto& { return cache_->temperature; }, + * [this]() { return expensive_temperature_calculation(); } + * ); + * } + * + * void add_to_output(double contribution) { + * set_output( + * [this]() -> auto& { return output_variable; }, + * contribution + * ); + * } + * }; + * @endcode + * + * @tparam CacheType A type derived from cache_base that provides storage + * for cached values. Must satisfy the CacheRules concept. + */ + +struct cache_base +{ + cache_base() = default; + int64_t last_timestep = -1; + + template + static constexpr T default_value() { + if constexpr (std::is_floating_point_v) + return std::numeric_limits::quiet_NaN(); + else if constexpr (std::is_integral_v) + return std::numeric_limits::min(); + else + return T{}; + }; + +}; + +namespace data_base_concepts { + template + concept CacheRules = std::derived_from; + + template + concept ValueRules = + requires(V v) { + {v()} -> std::same_as>; + }; + + template + concept OutputRules = ValueRules && + requires(O o,const T t) + { + {o()} -> std::convertible_to; + {o() += t}; + }; +}; + +namespace pt = boost::property_tree; + +template +class data_base { + + template + bool constexpr is_unset(T t) + { + if constexpr (!std::is_floating_point_v) + return t == cache_base::default_value(); + else + return std::isnan(t); + }; + + bool is_stale(); + + void init_cache(); + +protected: + + data_base(const mesh_elem& face_in, const std::shared_ptr param, + const pt::ptree* cfg); + // Test constructor to skip checks of valid mesh_elem + data_base(const mesh_elem& face_in, const std::shared_ptr param, + const pt::ptree* cfg, bool istest); + ~data_base() {}; + + const mesh_elem face{nullptr}; + const std::shared_ptr global_param; + const pt::ptree* const cfg; + mutable std::optional cache_; + + template + void update_value(Value&& value, const Fetch& fetch); + + template Output> + void set_output(Output&& output,const T t); + +public: + void reset_cache() { cache_.reset(); }; + const std::optional& get_cache() const { return cache_; }; +}; + +template +void data_base::init_cache() { + if (!cache_ || is_stale()) { + cache_.emplace(); + cache_->last_timestep = global_param->timestep_counter; + } +} + +template +bool data_base::is_stale() +{ + return cache_->last_timestep != global_param->timestep_counter; +} + +template +data_base::data_base(const mesh_elem& face_in, const std::shared_ptr param, + const pt::ptree* cfg) : face(face_in), global_param(param), cfg(cfg) +{ + if (!face->is_valid()) + throw std::invalid_argument("Face handle points to an invalid face"); + + if (!global_param) + throw std::invalid_argument("global parameter holder is null"); +}; + +template +data_base::data_base(const mesh_elem& face_in, const std::shared_ptr param, + const pt::ptree* cfg, const bool istest) : face(face_in), global_param(param), cfg(cfg) +{ + if (!istest) + std::invalid_argument("data_base test constructor called with False istest flag. Should be true."); +}; +template +template +void data_base::update_value(Value&& value, const Fetch& fetch) { + init_cache(); + + auto& V = value(); + if ( is_unset(V) ) + { + V = fetch(); + } + +}; + +template +template Output> +void data_base::set_output(Output&& output,const T t) +{ + init_cache(); + + output() += t; +}; diff --git a/src/modules/submodules/katabatic_melt_energy.cpp b/src/modules/submodules/katabatic_melt_energy.cpp new file mode 100644 index 00000000..f587deb5 --- /dev/null +++ b/src/modules/submodules/katabatic_melt_energy.cpp @@ -0,0 +1,46 @@ +#include "katabatic_melt_energy.hpp" +#include "PhysConst.h" +#include +#include + +namespace katabatic_melt_energy +{ + const Units::m_per_s bulk_coefficient(const Params& p, const Units::TempDiff deficit, + const Units::LapseRateSI gamma, + const Units::Kelvin glacier_temperature) + { + if (gamma == Units::LapseRateSI{0.0} || + p.prandtl == 0.0 || + glacier_temperature == Units::Kelvin{0.0}) + { + std::string err = std::format("Divide by zero in katabatic_melt_energy::bulk_coefficient for gamma = {} K/m, Pr = {}, T = {} K",gamma.value,p.prandtl,glacier_temperature.value); + throw std::runtime_error(err); + } + auto result = Units::m_per_s{p.k * + std::pow(p.k2,2.0) * deficit.value * + std::sqrt(p.g / + (glacier_temperature.value * gamma.value * p.prandtl))}; + //result.value = (0.01 + result.value)/2; + return result; + }; + + const Units::Watts_per_m2 sensible_heat(const Params& p, const Units::m_per_s coeff,const Units::TempDiff deficit,const Units::DensitySI air_density) + { + auto m = Units::Watts_per_m2{air_density.value * p.heat_capacity_air + * coeff.value * deficit.value}; + + return m; + }; + + const Units::Watts_per_m2 latent_heat(const Params& p, const Units::m_per_s coeff,const Units::Pa vapour_pressure_deficit, + const Units::Kelvin glacier_temperature, const Units::DensitySI air_density,const Units::Pa air_pressure) + { + Units::Celsius Temp{glacier_temperature.value - 273.15}; + auto m = Units::Watts_per_m2{p.molecular_wt_ratio * air_density.value * + PhysConst::Lv(Temp) * coeff.value * vapour_pressure_deficit.value / air_pressure.value}; + + return m; + + }; + +}; diff --git a/src/modules/submodules/katabatic_melt_energy.hpp b/src/modules/submodules/katabatic_melt_energy.hpp new file mode 100644 index 00000000..4df33223 --- /dev/null +++ b/src/modules/submodules/katabatic_melt_energy.hpp @@ -0,0 +1,78 @@ +#pragma once + +#include "Atmosphere.h" +#include "base_step.hpp" +#include "PhysConst.h" +#include + +namespace Units = PhysConst::units; + +namespace katabatic_melt_energy +{ + struct Params + { + const double heat_capacity_air = PhysConst::Cp(); + const double molecular_wt_ratio = PhysConst::em(); + const double g = PhysConst::g(); + + double prandtl = 5.0; + double k = 4e-4; + double k2 = 1.0; + size_t seconds_per_step = 3600.0; + }; + + const Units::m_per_s bulk_coefficient( const Params&, const Units::TempDiff deficit, + const Units::LapseRateSI gamma, + const Units::Kelvin glacier_temperature); + + const Units::Watts_per_m2 sensible_heat(const Params&, const Units::m_per_s coeff,const Units::TempDiff deficit,const Units::DensitySI air_density); + + const Units::Watts_per_m2 latent_heat(const Params&, const Units::m_per_s coeff,const Units::Pa vapour_pressure_deficit, + const Units::Kelvin glacier_temperature, const Units::DensitySI air_density, const Units::Pa air_pressure); + + template + concept KatabaticData = requires(T& t) + { + { t.glacier_temperature() } -> std::same_as; + { t.air_temperature() } -> std::same_as; + { t.air_pressure() } -> std::same_as; + { t.vapour_pressure() } -> std::same_as; + { t.vapour_pressure_surface() } -> std::same_as; + { t.lapse_rate() } -> std::same_as; + + { t.latent_heat(std::declval()) } -> std::same_as; + { t.sensible_heat(std::declval()) } -> std::same_as; + }; + + template + class Model : public base_step,Data> + { + Params p; + public: + void execute_impl(Data& d) const + { + auto air_pressure = d.air_pressure(); + auto vapour_pressure = d.vapour_pressure(); + auto glacier_temperature = d.glacier_temperature(); + auto air_temperature = d.air_temperature(); + auto lapse_rate = d.lapse_rate(); + + auto rho_air = Atmosphere::air_density(air_pressure,glacier_temperature,vapour_pressure); + + Units::TempDiff T_deficit{air_temperature.value - glacier_temperature.value}; + + auto K = bulk_coefficient(p,T_deficit,lapse_rate,glacier_temperature); + + auto Q_sensible = sensible_heat(p,K,T_deficit,rho_air); + + Units::Pa vapour_pressure_deficit{vapour_pressure.value - d.vapour_pressure_surface().value}; + + auto Q_latent = latent_heat(p,K,vapour_pressure_deficit,glacier_temperature,rho_air,air_pressure); + + d.latent_heat(Q_latent.value); + d.sensible_heat(Q_sensible.value); + }; + + Params& get_params() { return p; }; + }; +}; diff --git a/src/modules/submodules/melt_routing_glacier.cpp b/src/modules/submodules/melt_routing_glacier.cpp new file mode 100644 index 00000000..18f07a8a --- /dev/null +++ b/src/modules/submodules/melt_routing_glacier.cpp @@ -0,0 +1,90 @@ +#include "melt_routing_glacier.hpp" +#include +#include + +namespace GlacierRouting +{ + LinearReservoir::LinearReservoir(const double k, const double dt) + : c0(get_c0(k,dt)), c1(get_c1(k,dt)) { + last_input = 0.0; + last_output = 0.0; + }; + static void verifyStability(const double k, const double dt) + { + if ( k > 0.0 && k <= dt / 2 ) + { + std::string err = std::format("In GlacierRouting, numerical stability conditions is not satisfied." + "k = {} s, and dt = {} s.",k,dt); + throw std::runtime_error(err); + } + }; + const double LinearReservoir::get_c0(const double k, const size_t dt) + { + verifyStability(k,dt); + return dt / (2 * k + dt); + }; + const double LinearReservoir::get_c1(const double k, const size_t dt) + { + return ( 2*k - dt ) / ( 2*k + dt ); + }; + + const double LinearReservoir::step(const double input) { + auto result = c0 * ( input + last_input ) + + c1 * last_output; + + last_input = input; + last_output = result; + + return result; + }; + + GlacierReservoir::GlacierReservoir(const Params* _p) : p(_p) + { + reservoir_chain = [&]{ + std::vector result; + result.reserve(p->chain_length); + for (size_t i = 0; i < p->chain_length-1; ++i) + { + // First chain_length - 1 elements are linear reservoirs but with k=0 + // This means that Q_J+1 = I_J+1 + result.emplace_back(0.0,p->seconds_per_step); + } + result.emplace_back(p->k,p->seconds_per_step); + return result; + }(); + + old_outputs.resize(p->chain_length - 1, 0.0); + }; + + const double GlacierReservoir::step(const double input) { + // Process reservoirs in sequence: each reservoir consumes the previous + // reservoir's output from the last timestep (old_outputs) and produces + // a new output that becomes the input for the next reservoir for the next timestep. + // The number of old_outputs is one less than reservoirs because the first + // reservoir takes the external input, not a previous output. + + auto old_output_upchain = old_outputs.begin(); + auto current_reservoir = reservoir_chain.begin(); + + double output = current_reservoir->step(input); + ++current_reservoir; + + // Note: Extra ++ on reservoir iterator moves it one ahead from the old_outputs iterator + // This is intentional and syncs up the right element of old_outputs with reservoir_chain + // as described in the first comment. + for (; current_reservoir != reservoir_chain.end() && old_output_upchain != old_outputs.end();) + { + auto output_temp = current_reservoir->step(*old_output_upchain); + *old_output_upchain = output; + output = output_temp; + ++old_output_upchain; + ++current_reservoir; + }; + + return output; + }; + + State::State(const Params* p) : snow_route(p), firn_route(p), ice_route(p) {}; + // Delete copy constructors + +}; diff --git a/src/modules/submodules/melt_routing_glacier.hpp b/src/modules/submodules/melt_routing_glacier.hpp new file mode 100644 index 00000000..c2024ebb --- /dev/null +++ b/src/modules/submodules/melt_routing_glacier.hpp @@ -0,0 +1,92 @@ +#pragma once + +#include "base_step.hpp" +#include "PhysConst.h" +#include + +namespace Units = PhysConst::units; + +namespace GlacierRouting +{ + struct Params { + size_t chain_length = 10u; + double k = 4e3; + size_t seconds_per_step = 3600u; + }; + + class LinearReservoir { + static const double get_c0(const double k, const size_t dt); + static const double get_c1(const double k, const size_t dt); + const double c0; + const double c1; + + double last_input; + double last_output; + public: + LinearReservoir(const double k, const double dt); + + const double step(const double input); + }; + + class GlacierReservoir + { + const Params* const p; + std::vector reservoir_chain; + std::vector old_outputs; + public: + GlacierReservoir(const Params* _p); + + const double step(const double input); + }; + + struct State + { + State(const Params* p); + // Delete copy constructors + State(const State&) = delete; + State& operator=(const State&) = delete; + GlacierReservoir snow_route; + GlacierReservoir firn_route; + GlacierReservoir ice_route; + }; + + template + concept GlacierRoutingData = requires(T& t) + { + { t.get_state() } -> std::same_as; + { t.snowmelt() } -> std::same_as; + { t.firnmelt() } -> std::same_as; + { t.icemelt() } -> std::same_as; + + { t.snowmelt_delayed(std::declval()) } -> std::same_as; + { t.firnmelt_delayed(std::declval()) } -> std::same_as; + { t.icemelt_delayed(std::declval()) } -> std::same_as; + { t.total_delayed(std::declval()) } -> std::same_as; + }; + + template + class Model : public base_step,Data> + { + Params p; + public: + void execute_impl(Data& d) const + { + auto snowmelt = d.snowmelt(); + auto firnmelt = d.firnmelt(); + auto icemelt = d.icemelt(); + auto& s = d.get_state(); + + auto snow = s.snow_route.step(snowmelt.value); + auto firn = s.firn_route.step(firnmelt.value); + auto ice = s.ice_route.step(icemelt.value); + + d.snowmelt_delayed(snow); + d.firnmelt_delayed(firn); + d.icemelt_delayed(ice); + + d.total_delayed(snow + firn + ice); + }; + + Params& get_params() { return p; } + }; +} diff --git a/src/modules/submodules/water_flux.cpp b/src/modules/submodules/water_flux.cpp new file mode 100644 index 00000000..591fc765 --- /dev/null +++ b/src/modules/submodules/water_flux.cpp @@ -0,0 +1,59 @@ +#include "water_flux.hpp" +#include "PhysConst.h" +namespace Units = PhysConst::units; +template +water_flux::water_flux(const double mm_per_s) : _value_mm_per_s(mm_per_s) {}; + +template water_flux water_flux::from_W_per_m_squared(const double Q,const Units::Celsius T) +{ + auto val_mm_per_s = Q / (PhysConst::water_reference_density() * PhysConst::Lv(T)) * mm_per_m; + + return water_flux(val_mm_per_s); +}; + +template water_flux water_flux::from_mm_per_dt(const double WE,const size_t dt) +{ + // mm / step (seconds/step) -> mm / si + auto val_mm_per_s = WE / dt; + return water_flux(val_mm_per_s); +}; + +template water_flux water_flux::from_mm_per_s(const double WE) +{ + return from_mm_per_dt(WE,1); +}; + +template water_flux water_flux::from_m_per_s(const double WE) +{ + constexpr auto MM_PER_M = 1000.0; + auto val_mm_per_s = WE * MM_PER_M; + return from_mm_per_s(val_mm_per_s); +}; + +template +const double water_flux::mm_per_dt(const size_t dt) const +{ + return _value_mm_per_s * dt; +}; + +template +const double water_flux::mm_per_s() const +{ + return mm_per_dt(1.0); +}; + +template +const double water_flux::m_per_s() const +{ + return _value_mm_per_s / 1000.0; +}; + +template +const double water_flux::W_per_m_squared(const Celsius T) const +{ + PhysConst::units::Celsius temp{T}; + return _value_mm_per_s * PhysConst::water_reference_density() * PhysConst::Lv(temp) / mm_per_m; +}; + +template +bool water_flux::empty() const { return _value_mm_per_s == 0.0; }; diff --git a/src/modules/submodules/water_flux.hpp b/src/modules/submodules/water_flux.hpp new file mode 100644 index 00000000..1e42c242 --- /dev/null +++ b/src/modules/submodules/water_flux.hpp @@ -0,0 +1,131 @@ +#pragma once +#include "PhysConst.h" +#include + +namespace Units = PhysConst::units; +enum class FluxType { +latent, +sensible, +Default +}; + +template +class water_flux +{ + static_assert(!std::constructible_from && + !std::constructible_from); + + double _value_mm_per_s; + explicit water_flux(const double mm_per_s); + + static constexpr auto mm_per_m = 1000.0; + + using Celsius = PhysConst::units::Celsius; +public: + + // Defaults to freezing temperature + static water_flux from_W_per_m_squared(const double Q, + const Celsius T = Celsius{0.0} ); + static water_flux from_mm_per_dt(const double WE,const size_t dt); + static water_flux from_mm_per_s(const double WE); + static water_flux from_m_per_s(const double WE); + + const double mm_per_dt(const size_t dt) const ; + const double W_per_m_squared(const Celsius T = Celsius{0.0}) const + requires (type == FluxType::latent);//std::is_same_v; + + const double mm_per_s() const; + const double m_per_s() const; + bool empty() const; + + water_flux operator+(const water_flux& other) const + { + return water_flux(_value_mm_per_s + other._value_mm_per_s); + }; + + water_flux operator-(const water_flux& other) const + { + return water_flux(_value_mm_per_s - other._value_mm_per_s); + }; + + void operator+=(const water_flux& other) + { + _value_mm_per_s += other._value_mm_per_s; + }; + + void operator-=(const water_flux& other) + { + _value_mm_per_s -= other._value_mm_per_s; + }; +}; + +template +water_flux::water_flux(const double mm_per_s) : _value_mm_per_s(mm_per_s) {}; + +template +water_flux water_flux::from_W_per_m_squared(const double Q) +{ + auto val_mm_per_s = 0.0; + if constexpr (type == FluxType::latent) + val_mm_per_s = Q / (PhysConst::water_reference_density() * PhysConst::Lf()) * mm_per_m; + else + { + static_assert(type != FluxType::latent && + "W/m^2 functions are not supported for sensible heat or other quantities" + "If a sensible heat conversion to mm/dt is required, it must be done manually" + "Carefully consider if this path is physical as sensible heat == temperature change, not a mass flux."); + } + + return water_flux(val_mm_per_s); +}; + +template +water_flux water_flux::from_mm_per_dt(const double WE,const size_t dt) +{ + // mm / step (seconds/step) -> mm / si + auto val_mm_per_s = WE / dt; + return water_flux(val_mm_per_s); +}; + +template +water_flux water_flux::from_mm_per_s(const double WE) +{ + return from_mm_per_dt(WE,1); +}; + +template +water_flux water_flux::from_m_per_s(const double WE) +{ + constexpr auto MM_PER_M = 1000.0; + auto val_mm_per_s = WE * MM_PER_M; + return from_mm_per_s(val_mm_per_s); +}; + +template +const double water_flux::mm_per_dt(const size_t dt) const +{ + return _value_mm_per_s * dt; +}; + +template +const double water_flux::mm_per_s() const +{ + return mm_per_dt(1.0); +}; + +template +const double water_flux::m_per_s() const +{ + return _value_mm_per_s / 1000.0; +}; + +template +const double water_flux::W_per_m_squared(const Celsius T) const +requires (type == FluxType::latent) +{ + PhysConst::units::Celsius temp{T}; + return _value_mm_per_s * PhysConst::water_reference_density() * PhysConst::Lv(temp) / mm_per_m; +}; + +template +bool water_flux::empty() const { return _value_mm_per_s == 0.0; }; diff --git a/src/physics/Atmosphere.cpp b/src/physics/Atmosphere.cpp index 380e49bd..64b4f581 100644 --- a/src/physics/Atmosphere.cpp +++ b/src/physics/Atmosphere.cpp @@ -22,6 +22,9 @@ // #include "physics/Atmosphere.h" +#include "PhysConst.h" +#include +#include namespace Atmosphere { @@ -79,4 +82,18 @@ namespace Atmosphere return Es; } + const Units::DensitySI air_density(const Units::Pa p, const Units::Kelvin T, + const Units::Pa e) + { + if (p < e) + { + std::string err = std::format("Air pressure, {} Pa, can never be less than the vapour pressure, {} Pa, which is a partial pressure",p.value,e.value); + throw std::runtime_error(err); + }; + + Units::DensitySI result{p.value / (T.value * PhysConst::RgasDry() ) * + (1 - (1 - PhysConst::em()) * e.value / p.value)}; + + return result; + }; } diff --git a/src/physics/Atmosphere.h b/src/physics/Atmosphere.h index c0035fcd..80287e8a 100644 --- a/src/physics/Atmosphere.h +++ b/src/physics/Atmosphere.h @@ -22,9 +22,11 @@ #include #include +#include "PhysConst.h" #pragma once +namespace Units = PhysConst::units; namespace Atmosphere { /********* Atmosphere ************/ @@ -43,6 +45,17 @@ namespace Atmosphere { double corr_precip_slope(double p, double slope); double saturatedVapourPressure(const double& T); + + /** + * @brief Calculation of air density + * @param p air pressure (Pa) + * @param T temperature (K) + * @param e vapour pressure (Pa) + * @return Air Density (kg/m^3) + */ + const Units::DensitySI air_density(const Units::Pa p, + const Units::Kelvin T, + const Units::Pa e); } diff --git a/src/physics/PhysConst.cpp b/src/physics/PhysConst.cpp new file mode 100644 index 00000000..8855bf1b --- /dev/null +++ b/src/physics/PhysConst.cpp @@ -0,0 +1,9 @@ +#include "PhysConst.h" + +namespace PhysConst +{ + const double Lv(const units::Celsius T) + { + return Lv() - 0.002361 * T.value * 1e6; + }; +}; diff --git a/src/physics/PhysConst.h b/src/physics/PhysConst.h index 52526bfd..2323205d 100644 --- a/src/physics/PhysConst.h +++ b/src/physics/PhysConst.h @@ -19,42 +19,150 @@ * along with Canadian Hydrological Model. If not, see * . */ +#pragma once #include +#include -#pragma once namespace PhysConst { - /********* Physical Constants ************/ - const double sbc = mio::Cst::stefan_boltzmann; // 5.670373e-8 (W m-2 K-4) - const double Rgas = mio::Cst::gaz_constant_dry_air; // (W m-2 K-4) - const double kappa = 0.4; // Von Karman constant - const double Ls = mio::Cst::l_water_sublimation; // 2.838e6; // Latent heat of sublimation // ((J kg-1) - const double Cp = mio::Cst::specific_heat_air; // 1004.67; // (J K-1), see Stull "Meteorology for scientists and engineers" p44 - const double Ci = mio::Cst::specific_heat_ice; // = 2100.0 (J K-1), at 0C + namespace units + { + template + struct base { + T value; + + static_assert(std::is_arithmetic_v,"PhysConst::units::base template must use a built-in integral or floating-point type"); - // From PBSM Pomeroy et al. 1993 - const double M = 18.01; // M is the molecular weight of water (18.01 kg kmol-l); - const double R = mio::Cst::gaz_constant; // R is the universal gas constant (8313 J mol-1 K-1); FYI: unit error in Pomeroy et al. 1993 - const double rho_ice = 917; //{density of ice (kg/m3)} + // Prevent derived classes from getting their own synthesized comparisons + void operator+=(const Derived& other) + { + static_cast(*this).value += other.value; + } + void operator-=(const Derived& other) + { + static_cast(*this).value -= other.value; + } + friend auto operator<=>(const Derived& lhs,const Derived& rhs) + { + return lhs.value <=> rhs.value; + } + friend bool operator==(const Derived& lhs,const Derived& rhs) + { + return lhs.value == rhs.value; + } + }; - // From CRHM_constants in common.h - const double Cs = 1.28E+06; // (J/m3/K) volumetric heat capacity of soil - //const double Cp = 1005; // (J/kg/K) volumetric heat capacity of dry air - //const double Rgas = 287.0; // Gas constant for dry air (J/kg/K) - //const double Tm = 273.15; // Melting point (K) + // Temperature + static inline constexpr double freeze_point = 273.15; + struct Celsius; + struct Kelvin + : public base { + Kelvin(const Celsius& c); + explicit Kelvin(const double in) { value = in; }; + }; + struct Celsius + : public base { + Celsius(const Kelvin& k) + { + value = k.value - freeze_point; + }; + explicit Celsius(const double in) {value = in;}; + }; + inline Kelvin::Kelvin(const Celsius& c) { value = c.value + freeze_point; }; + struct TempDiff + : public base { + TempDiff(const Kelvin &c) {value = c.value; }; + TempDiff(const Celsius &c) {value = c.value; }; + explicit TempDiff(const double in) { value = in; }; + }; - //const double Ls = 2.845e6; // Latent heat of sublimation (J/kg) - const double Lv = 2.50e6; // Latent heat of vaporization (J/kg) - const double Lf = 0.334e6; // Latent heat of fusion (J/kg) - //const double kappa = 0.4; - //const double sbc = 5.67E-8; // Stephan-Boltzmann constant W/m^2/k4 - const double SB = 4.899e-09; // Stephan-Boltzmann constant MJ/m^2-d + //depths + struct Milimetres + : public base {}; + struct Metres + : public base {}; - const double em = 0.622; // + // depth / time + struct m_per_s + : public base {}; + + // mass + struct Kg_per_m3 + : public base {}; + using DensitySI = Kg_per_m3; -} + // pressure + struct Pa + : public base {}; + + // Lapse rate + struct Degree_per_metre + : public base {}; + using LapseRateSI = Degree_per_metre; + // Energy rate + struct Watts_per_m2 + : public base {}; + }; + + + /********* Physical Constants ************/ + // Stefan-Boltzmann constant (W m-2 K-4) + inline const double sbc() { return mio::Cst::stefan_boltzmann; } + + // Gas constant for dry air (J kg-1 K-1) + inline const double RgasDry() { return mio::Cst::gaz_constant_dry_air; } + // Gas constant for water vapour (J kg-1 K-1) + inline const double RgasVapour() { return mio::Cst::gaz_constant_water_vapor; } + // Von Karman constant (dimensionless) + constexpr double kappa() { return 0.4; } + + // Latent heat of sublimation (J kg-1) + inline const double Ls() { return mio::Cst::l_water_sublimation; } + + // Specific heat of dry air at constant pressure (J kg-1 K-1) + inline const double Cp() { return mio::Cst::specific_heat_air; } + + // Specific heat of ice (J kg-1 K-1) + inline const double Ci() { return mio::Cst::specific_heat_ice; } + + inline const double Cw() { return mio::Cst::specific_heat_water; } + + // Molecular weight of water (kg kmol-1) + constexpr double M() { return 18.01; } + + // Universal gas constant (J mol-1 K-1) + inline const double R() { return mio::Cst::gaz_constant; } + + // Density of ice (kg m-3) + constexpr double rho_ice() { return 917.0; } + + // Volumetric heat capacity of soil (J m-3 K-1) + constexpr double Cs() { return 1.28E+06; } + + // Latent heat of vaporization, constant value (J kg-1) + constexpr double Lv() { return 2.501e6; } + + // Latent heat of vaporization with temperature dependence (J kg-1) + const double Lv(const units::Celsius T); + + // Latent heat of fusion (J kg-1) + constexpr double Lf() { return 0.334e6; } + + // Stefan-Boltzmann constant, day units (MJ m-2 d-1) + constexpr double SB() { return 4.899e-09; } + + // Ratio of molecular weights of water vapor to dry air (dimensionless) + constexpr double em() { return 0.622; } + + // Acceleration due to gravity (m s-2) + constexpr auto g() {return 9.81;} + + // Reference density of water (kg m-3) + constexpr double water_reference_density() { return 1000.0; } + +} diff --git a/src/tests/submoduletests/PhysConst.cpp b/src/tests/submoduletests/PhysConst.cpp new file mode 100644 index 00000000..e69de29b diff --git a/src/tests/submoduletests/test_bisection.cpp b/src/tests/submoduletests/test_bisection.cpp new file mode 100644 index 00000000..519c54f3 --- /dev/null +++ b/src/tests/submoduletests/test_bisection.cpp @@ -0,0 +1,51 @@ +#include "math/Bisection.hpp" +#include "gtest/gtest.h" +#include +#include + +using namespace math; + +class BisectionTest : public ::testing::Test +{ + +protected: + static constexpr auto foo = [](const double x) { + return x * x - 4.0; + }; +}; + +TEST_F(BisectionTest,NoRootGuaranteed) +{ + auto result = bisection(foo,-1.0,1.0); + + EXPECT_TRUE(std::isnan(result.value)); + EXPECT_EQ(result.root,Root::NoRootGuaranteed); +}; + +TEST_F(BisectionTest,DoesNotConverge) +{ + constexpr size_t matched_exponent = 3; + static constexpr auto slow = [](const double x) { + return std::log(x) - matched_exponent; + }; + + auto result = bisection(slow,1.0,200.0,1e-12,10); + + EXPECT_EQ(result.root,Root::DidNotConverge); + EXPECT_EQ(result.iterations,10); + + result = bisection(slow,std::exp(matched_exponent)*0.8,std::exp(matched_exponent)*1.2,1e-3); + + EXPECT_EQ(result.root,Root::Found); + + EXPECT_EQ(result.value,std::exp(matched_exponent)); +} + +TEST_F(BisectionTest,XSquaredRoot) +{ + + auto result = bisection(foo,-1.0,3.0); + + EXPECT_EQ(result.root,Root::Found); + EXPECT_NEAR(result.value,2.0,1e-10); +}; diff --git a/src/tests/submoduletests/test_data_base.cpp b/src/tests/submoduletests/test_data_base.cpp new file mode 100644 index 00000000..b27d2f90 --- /dev/null +++ b/src/tests/submoduletests/test_data_base.cpp @@ -0,0 +1,127 @@ +#include +#include "data_base.hpp" +#include "triangulation.hpp" +#include +// Mock Cache for Testing +struct MockCache : public cache_base { + double value = default_value(); +}; + +// Mock data class +class data : public data_base +{ +public: + data(mesh_elem face, std::shared_ptr param, pt::ptree& cfg) + : data_base(face,param,cfg,true) {}; + ~data() {}; + + // update_value is private to data and therefore must be access in a function + // rather than use directly in the test. + template + void update(V&& v, L&& l) + { + update_value(v, l); + }; + + // Same as update_value + template + void output(V&& output, const T& t) + { + set_output(output,t); + }; + + // Again, to access private member + double& get_value() + { + return cache_->value; + }; +}; + +// Test Fixture +class DataBaseTest : public ::testing::Test { +protected: + void SetUp() override { + mock_global = std::make_shared(); + } + + + // mock components for the data constructor + mesh_elem mock_face; + std::shared_ptr mock_global; + pt::ptree mock_cfg; + + // Lambdas for testing, standing in for face access + static inline auto A = []() -> auto& { static double a = 42.0; return a;}; + static inline auto B = []() -> auto& { static double b = 123.0; return b;}; + + +}; + + +TEST_F(DataBaseTest, CacheInitializesOnFirstUpdate) { + data db(mock_face, mock_global, mock_cfg); + + db.update([&]() -> auto& { return db.get_value();},A); + + EXPECT_FALSE(std::isnan(db.get_value())); + EXPECT_EQ(db.get_value(), A()); +} + +TEST_F(DataBaseTest, CacheRespectsTimestepChanges) { + data db(mock_face, mock_global, mock_cfg); + + // First call at timestep 0 + mock_global->timestep_counter = 0; + db.update([&]() -> auto& {return db.get_value();} ,B); + + // Second call at same timestep - should use cached value + db.update([&]() -> auto& {return db.get_value();} ,A); + EXPECT_EQ(db.get_value(), B()); // value remains B. Should not update from B to A + // at the same time step + + // Force cache reset by changing timestep + mock_global->timestep_counter = 1; + db.update([&]() -> auto& {return db.get_value();} ,A); + EXPECT_EQ(db.get_value(), A()); // Updates +} + +TEST_F(DataBaseTest, ManualCacheResetWorks) { + data db(mock_face, mock_global, mock_cfg); + + db.update([&]() -> auto& { return db.get_value(); },A); + db.reset_cache(); + + db.update([&]() -> auto& {return db.get_value(); },B); + EXPECT_EQ(db.get_value(), B()); +} + +TEST_F(DataBaseTest, SetOutputUpdatesValue) { + data db(mock_face, mock_global, mock_cfg); + double output = 0.0; + static constexpr double pi = 3.14; + EXPECT_FALSE(db.get_cache()); + db.output([&]() -> auto& {return output;}, pi); + EXPECT_EQ(output, pi); + EXPECT_TRUE(db.get_cache()); +} + +TEST_F(DataBaseTest, OnlyUpdatesNaNValues) { + data db(mock_face, mock_global, mock_cfg); + double test_value = 10.0; // Not NaN + auto L = [&]() -> auto& {return test_value;}; + db.update(L, A); + EXPECT_EQ(L(), 10.0); // Should remain unchanged +} + +TEST_F(DataBaseTest, ChecksDefaultValue) +{ + double double_var = MockCache::default_value(); + float float_var = MockCache::default_value(); + int int_var = MockCache::default_value(); + size_t size_t_var = MockCache::default_value(); + + EXPECT_TRUE(std::isnan(double_var)); + EXPECT_TRUE(std::isnan(float_var)); + EXPECT_TRUE(int_var == std::numeric_limits::min()); + EXPECT_TRUE(size_t_var == std::numeric_limits::min()); +}; diff --git a/src/tests/submoduletests/test_glacier.cpp b/src/tests/submoduletests/test_glacier.cpp new file mode 100644 index 00000000..eaf55fce --- /dev/null +++ b/src/tests/submoduletests/test_glacier.cpp @@ -0,0 +1,1668 @@ +#include +#include +#include +#include "Glacier.hpp" +#include "PhysConst.h" +#include +#include +#include + +using namespace Glacier; +inline constexpr auto MARGIN = 1e-12; +/* + * ParamsTest: Verifies the Params struct initialization. + */ + +TEST(ParamsTest,ParamsSetByPhysConst) +{ + Params p; + + EXPECT_EQ(p.water_density,PhysConst::water_reference_density()); + EXPECT_EQ(p.heat_capacity_air,PhysConst::Cp()); + EXPECT_EQ(p.latent_heat_vapour,PhysConst::Lv()); + EXPECT_EQ(p.gas_constant_dry,PhysConst::RgasDry()); + EXPECT_EQ(p.gas_constant_vapour,PhysConst::RgasVapour()); + EXPECT_EQ(p.molecular_wt_ratio,PhysConst::M()); + + auto _ = p.densify_version; +}; + +/* + * This function exists so that changes to the default values in the + * Glacier::Params struct can be changed without impacting the test success. + */ +[[maybe_unused]] static inline Params test_default_params() +{ + Params p; + p.densify_version = DensifyVersion::HerronLangway; + p.small_increment = 25.0; + p.big_increment = 50.0; + p.critical_density = 550.0; + p.firn_to_ice_density = 830.0; + p.thermal_factor = 0.95; + p.seconds_per_step = 3600.0; + return p; +}; + +/* + * DensificationTest: Verifies the calculations of densification of firn. + * + * Two methods: Linear method, and method based on Herron and Langway (1980) + */ + +class DensificationTest : public ::testing::Test +{ +protected: + Params p; + + const double critical_density = 550.0; + +}; + +/* + * DensificationTest - HerronLangway: Tests the Herron and Langway (1980) firn + * densification algorithm. + * + */ +class HerronLangwayTest : public DensificationTest +{ + static_assert(PhysConst::rho_ice() == 917.0, + "The original model assumes this specific value for ice density and cannot be modified.\n" + "Consider refactoring if PhysConst::rho_ice() must be modified.\n" + "I would suggest creating a ice_density constexpr function in Glacier.hpp in the Glacier namespace."); +}; + +TEST_F(HerronLangwayTest,ValidationAgainstPrecalculated) +{ + std::array heights{1.0,1.2,0.9,0.5}; + std::array densities{250.0,320.0,400.0,810.0}; + std::array glacier_temperature{273.15,250.7,250.8,200.1}; + std::array Accumulate_rate{0.11,0.265,0.4,0.022}; + const std::array new_density{ + 271.4597957366285, + 356.16647094729564, + 454.56304130714364, + 688.3589090630271 + }; + + double depth = 0.0; + for (size_t layer = 0; layer < densities.size(); ++layer) + { + depth += heights[layer]; + auto result = Densification::HerronLangway(depth,densities[layer],glacier_temperature[layer],Accumulate_rate[layer],critical_density); + + EXPECT_DOUBLE_EQ(result,new_density[layer]); + }; + +}; + +TEST_F(HerronLangwayTest,MaxDensityThrows) +{ + double height = 1.0; + double density = PhysConst::rho_ice() * 1.1; + + // density >= rho_ice = 917.0 kg/m^3 should not be permitted. + EXPECT_THROW( + Densification::HerronLangway(height,density,273.15,1.0,critical_density), + std::runtime_error); + EXPECT_THROW( + Densification::HerronLangway(height,PhysConst::rho_ice(),273.15,1.0,critical_density), + std::runtime_error); + +}; + +/* + * DensificationTest - Linear: Validates the linear densification algorithm used + * for testing and baseline comparisons. + * + */ + +TEST_F(DensificationTest,Linear) +{ + std::array densities{250.0,400.0,725.0}; + const double small_increment = 25.0; + const double big_increment = 55.0; + + auto result = Densification::Linear(densities[0],small_increment); + EXPECT_EQ(result,densities[0] + small_increment); + + result = Densification::Linear(densities[1],small_increment); + EXPECT_EQ(result,densities[1] + small_increment); + + result = Densification::Linear(densities[2],big_increment); + EXPECT_EQ(result,densities[2] + big_increment); +}; + +/*/ + * FirnLayerTest: Class containing the data nad behaviour of a single layer of firn + */ + +class FirnLayerTest : public ::testing::Test +{ +protected: + // Testing on a single layer, constructed in each test + std::unique_ptr layer; + Params p = test_default_params(); + + // placeholder to be allocated later and passed to the layer object + Units::Milimetres water_equivalent; + + // Height is metres, density is kg/m^3 + static_assert(std::derived_from); + static_assert(std::derived_from); +}; + +// Layer constructs from Units::Milimetres struct and returns the same value with the same type +TEST_F(FirnLayerTest,ConstructsFromMilimetreStruct) +{ + water_equivalent.value = 30.0; + layer = std::make_unique(water_equivalent); + + auto WE = layer->water_equivalent(); + + EXPECT_EQ(WE.value,water_equivalent.value); + +}; + +// Construct from Layer::Height and Layer::Density objects. +// Reproduce WE using WE formula: Dingman (2008), equation (5-12) +TEST_F(FirnLayerTest,ConstructsFromHeightandDensity) +{ + Layer::Height height{1.0}; + Layer::Density density{300.0}; + layer = std::make_unique(height,density); + + auto WE = layer->water_equivalent(); + + // returns water equivalent from original height/density + auto constexpr MM_PER_M = 1000.0; + EXPECT_EQ(WE.value,height.value * density.value / p.water_density * MM_PER_M); +}; + +TEST_F(FirnLayerTest,PartialRemove) +{ + water_equivalent.value = 350.0; + layer = std::make_unique(water_equivalent); + + Units::Milimetres to_remove{150.0}; + layer->remove(to_remove); + + // Layer uses a bisection method upon construction + // to convert the WE to height and density. There will + // definitely be a loss of precision. + // + // adding a layer is typically only once a year, so its acceptable + // uncertainty. + EXPECT_NEAR(layer->water_equivalent().value,water_equivalent.value - to_remove.value,MARGIN); +}; + +// Remove function is just for reducing a layer by an amount and should not be used for +// removal, thats what pop_back() is for! +TEST_F(FirnLayerTest,ExcessiveRemoveThrowsException) +{ + water_equivalent.value = 250; + Units::Milimetres to_remove{water_equivalent.value * 2.0}; + layer = std::make_unique(water_equivalent); + + EXPECT_THROW(layer->remove(to_remove),std::invalid_argument); +}; + +// Same as above, shouldn't be equal to total either. +// A layer with no WE should be deleted as an object, not just set to zero +TEST_F(FirnLayerTest,RemoveEqualToTotalThrows) +{ + water_equivalent.value = 903.0; + Units::Milimetres to_remove{water_equivalent.value}; + to_remove += water_equivalent; + layer = std::make_unique(water_equivalent); + + EXPECT_THROW(layer->remove(to_remove),std::invalid_argument); +}; + +// +TEST_F(FirnLayerTest,LinearDensification) +{ + water_equivalent.value = 50.0; + layer = std::make_unique(water_equivalent); + + layer->densifyLinear(Layer::Density{p.small_increment},Layer::Density{p.big_increment},Layer::Density{p.critical_density}); + + EXPECT_EQ(layer->water_equivalent().value,water_equivalent.value); + +}; + +TEST_F(FirnLayerTest,HerronLangwayDensification) +{ + water_equivalent.value = 500.0; + layer = std::make_unique(water_equivalent); + + Units::Kelvin T{273.15}; + layer->densifyHerronLangway(T,0.0,1.2,Layer::Density{550.0}); + + EXPECT_DOUBLE_EQ(layer->water_equivalent().value,water_equivalent.value); +}; + +/*/ + * LayeredFirnTest: Test fixture for the LayeredFirn class, which manages a stratified column + * of firn layers. + * + * In the language of C++ containers: the array "front" is the bottom of the firn and the back is the top of the firn. + * + * Uses a double-ended queue (deque) so that removal and adding to front/back are equally efficient. + * deque is more efficient at adding new elements to the ends than std::vector. deque does not + * require that data is layed out sequentially, but vector does. Therefore, adding or removing a firn + * layer from the front/back doesn't require moving the vector. It simply frees the bytes. + * + * There are few loops over the layers, so it should be ok. + * + */ +class LayeredFirnTest : public ::testing::Test +{ +protected: + Params p = test_default_params(); + std::unique_ptr firn; + + + static constexpr auto MM_PER_M = 1000.0; + + std::vector> init_h_rho = { + {Layer::Height{0.6},Layer::Density{250.0}}, + {Layer::Height{0.3},Layer::Density{350.0}}, + {Layer::Height{0.2},Layer::Density{840.0}}, + {Layer::Height{0.1},Layer::Density{900.0}}, + }; + + + std::deque old_layers() + { + std::deque result; + + auto expected_firn = 0.0; + + for (auto init : init_h_rho) + { + result.push_front(Layer(init.first,init.second)); + }; + + return result; + }; +}; + +TEST_F(LayeredFirnTest,DefaultConstructNoLayeredFirn) +{ + firn = std::make_unique(&p); + + EXPECT_EQ(firn->water_equivalent().value,0.0); +}; + +TEST_F(LayeredFirnTest,ConstructFromLayersVector) +{ + std::vector init_layers = [this]() { + auto layers = old_layers(); + std::vector result(layers.begin(),layers.end()); + return result; + }(); + + auto expected_firn = 0.0; + + for (auto init : init_h_rho) + { + expected_firn += init.first.value * init.second.value / PhysConst::water_reference_density() * MM_PER_M; + } + LayeredFirn firn2(&p,init_layers); + firn = std::make_unique(&p,init_layers); + + EXPECT_EQ(firn->water_equivalent().value,expected_firn); + +}; + +TEST_F(LayeredFirnTest,Accumulate) +{ + constexpr Units::Milimetres WE{350.0}; + + firn = std::make_unique(&p); + + firn->accumulate(WE); + + auto result = firn->water_equivalent(); + + EXPECT_DOUBLE_EQ(result.value,WE.value); + +}; + +TEST_F(LayeredFirnTest,ConvertToIce) +{ + std::deque layers = old_layers(); + + firn = std::make_unique(&p,layers); + Units::Milimetres init_WE = firn->water_equivalent(); + + auto result = firn->convert_to_ice(); + + auto& new_layers = firn->get_layers(); + + EXPECT_EQ(firn->water_equivalent().value + result->value, + init_WE.value); + + EXPECT_EQ(layers.size() - 1,new_layers.size()); + +}; + +TEST_F(LayeredFirnTest,NoConvertToIce) +{ + std::deque layers = old_layers(); + std::reverse(layers.begin(),layers.end()); + + firn = std::make_unique(&p,layers); + Units::Milimetres init_WE = firn->water_equivalent(); + + auto result = firn->convert_to_ice(); + + auto& new_layers = firn->get_layers(); + + EXPECT_EQ(firn->water_equivalent().value,init_WE.value); + + EXPECT_EQ(layers.size(),new_layers.size()); + + EXPECT_EQ(layers.front().water_equivalent().value, + new_layers.front().water_equivalent().value); + +}; + + +/*/ + * LayeredFirnMeltTest: Specialized test fixture for testing melt calculations of firn for the + * variety of pathways + */ + +class LayeredFirnMeltTest : public LayeredFirnTest +{ +protected: + static constexpr std::array WE{200.0,150.0,303.0,525.0,201.0,25.0}; + std::deque layers; + + void SetUp() override { + constexpr size_t N = WE.size(); + + for (size_t n = 0; n < N; ++n) + { + Units::Milimetres water_eq{WE[n]}; + Layer layer(water_eq); + layers.push_back(layer); + }; + firn = std::make_unique(&p,layers); + }; +}; + +/* + * Unrelated to melt, but built deque here. + */ +TEST_F(LayeredFirnMeltTest,ConstructFromDeque) +{ + const auto sum = std::accumulate(WE.begin(),WE.end(),0.0); + + EXPECT_DOUBLE_EQ(sum,firn->water_equivalent().value); + + size_t ind = 0; + for (const auto layer : firn->get_layers()) + { + EXPECT_DOUBLE_EQ(layer.water_equivalent().value, WE[ind]); + ++ind; + } +}; + + +TEST_F(LayeredFirnMeltTest,PartialMeltSingleLayer) +{ + firn = std::make_unique(&p); + constexpr Units::Milimetres input{150.0}; + constexpr double melt_flux{75.0}; + firn->accumulate(input); + + // second argument is dt in units of seconds / step + // Must use p.seconds_per_step because internals do the same + // auto F = Units::Milimetres{melt_flux * p.seconds_per_step}; + auto F = Units::Milimetres{melt_flux}; + + auto melt_info = firn->melt(F); + + EXPECT_DOUBLE_EQ(melt_info.melt.value,melt_flux); + + EXPECT_EQ(melt_info.remaining_energy.value,0.0); + + auto& layers = firn->get_layers(); + + EXPECT_EQ(layers.size(),1); + EXPECT_EQ(layers.front().water_equivalent(),layers.back().water_equivalent()); + EXPECT_NEAR(layers.at(0).water_equivalent().value,input.value - melt_info.melt.value,MARGIN); +}; + +TEST_F(LayeredFirnMeltTest,PartialMeltManyLayer) +{ + const auto melt_energy_mm_per_dt = Units::Milimetres{705.0}; + const auto initial_WE = std::accumulate(WE.begin(),WE.end(),0.0); + const auto expected_final_WE = Units::Milimetres{initial_WE - melt_energy_mm_per_dt.value}; + + MeltInfo info = firn->melt(melt_energy_mm_per_dt); + + auto sum = 0.0; + auto excess = 0.0; + for (auto it = WE.rbegin(); it != WE.rend(); ++it) + { + sum += *it; + if (sum >= melt_energy_mm_per_dt.value) + { + excess = sum - melt_energy_mm_per_dt.value; + break; + } + } + + // amount melted is equal to input flux + EXPECT_DOUBLE_EQ(info.melt.value,melt_energy_mm_per_dt.value); + + // final WE is correct + EXPECT_DOUBLE_EQ(firn->water_equivalent().value, + expected_final_WE.value); + + // layers removed at END only. + // This isn't quite right. It will fail at the last iteration + // num_layers()? Might also be useful to "predict" number of layers. + const auto& layers = firn->get_layers(); + auto layers_it = layers.begin(); + auto WE_it = WE.begin(); + + while(layers_it != layers.end()-1 && WE_it != WE.end()) + { + EXPECT_DOUBLE_EQ((*layers_it).water_equivalent().value,*WE_it); + ++layers_it; + ++WE_it; + } + + EXPECT_NEAR((*layers_it).water_equivalent().value,excess,MARGIN); + ++layers_it; + EXPECT_TRUE(layers_it == layers.end()); + EXPECT_FALSE(WE_it == WE.end()); + +}; + +TEST_F(LayeredFirnMeltTest,MeltAllLayersExcess) +{ + constexpr double initial_WE = std::accumulate(WE.begin(),WE.end(),0.0); + constexpr auto melt_energy_mm_per_dt = initial_WE * 1.25; + double expected_final_WE = initial_WE - melt_energy_mm_per_dt; + + MeltInfo info = firn->melt(Units::Milimetres{melt_energy_mm_per_dt}); + + auto excess = melt_energy_mm_per_dt - info.remaining_energy.value; + + EXPECT_DOUBLE_EQ(info.melt.value,initial_WE); + + EXPECT_DOUBLE_EQ(firn->water_equivalent().value, + 0.0); + + EXPECT_NEAR(info.remaining_energy.value, + melt_energy_mm_per_dt - initial_WE,MARGIN); +}; + +TEST_F(LayeredFirnMeltTest,MeltAllLayersExact) +{ + constexpr double initial_WE = std::accumulate(WE.begin(),WE.end(),0.0); + constexpr auto melt_energy_mm_per_dt = initial_WE; + double expected_final_WE = initial_WE - melt_energy_mm_per_dt; + + MeltInfo info = firn->melt(Units::Milimetres{melt_energy_mm_per_dt}); + + auto excess = melt_energy_mm_per_dt - info.remaining_energy.value; + + EXPECT_DOUBLE_EQ(info.melt.value,initial_WE); + + EXPECT_DOUBLE_EQ(firn->water_equivalent().value, + 0.0); + + EXPECT_EQ(info.remaining_energy.value, + 0.0); + +}; + +/* + * IceTest: Test fixture for the Ice class, modelling ice as a bulk, single layer. + */ + +class IceTest : public ::testing::Test +{ +protected: + Params p = test_default_params(); + std::unique_ptr ice; +}; + +TEST_F(IceTest,ConstructEmpty) +{ + ice = std::make_unique(&p); + + EXPECT_EQ(ice->water_equivalent().value,0.0); +}; + +TEST_F(IceTest,ConstructFromEmptyMilimeterStruct) +{ + Units::Milimetres WE{0.0}; + ice = std::make_unique(&p,WE); + + EXPECT_EQ(ice->water_equivalent().value,0.0); +}; +TEST_F(IceTest,ConstructsFromMilimetreStruct) +{ + Units::Milimetres WE{500.34}; + ice = std::make_unique(&p,WE); + + EXPECT_EQ(ice->water_equivalent(),WE); +}; + +TEST_F(IceTest,Accumulate) +{ + Units::Milimetres WE{125.0}; + ice = std::make_unique(&p,WE); + + constexpr auto incoming = 340.5; + WE += Units::Milimetres{incoming}; + ice->accumulate(Units::Milimetres{incoming}); + + auto current_WE = ice->water_equivalent(); + + EXPECT_EQ(current_WE.value,WE.value); +}; + +TEST_F(IceTest,MeltSmall) +{ + Units::Milimetres WE{3030.0}; + ice = std::make_unique(&p,WE); + Units::Milimetres to_melt{125.0}; + + auto melt_info = ice->melt(to_melt); + + EXPECT_EQ(melt_info.melt,to_melt); + + EXPECT_EQ(ice->water_equivalent().value,WE.value - to_melt.value); + +}; + +TEST_F(IceTest,MeltBig) +{ + Units::Milimetres WE{1245.0}; + ice = std::make_unique(&p,WE); + Units::Milimetres to_melt{WE.value * 1.35}; + + auto melt_info = ice->melt(to_melt); + + EXPECT_EQ(melt_info.melt,WE); + + EXPECT_DOUBLE_EQ(melt_info.remaining_energy.value,WE.value*0.35); +}; + +class StateTest : public ::testing::Test +{ +protected: + Params p = test_default_params(); + + Ice ice() + { + Units::Milimetres ice_WE{100.0}; + return Ice{&p,ice_WE}; + }; + + LayeredFirn firn() + { + struct Firn_WE + { + Layer layer; + Units::Milimetres mm; + }; + + std::vector init {{ + {100.0},{200.0},{350.0},{75.0} + }}; + + std::vector layers; + + for (auto i : init) + { + layers.push_back(i); + } + + + return LayeredFirn{&p,layers}; + } + +}; + +TEST_F(StateTest,ConstructEmpty) +{ + State state(&p); + + auto WE = state.total_water_equiv(); + auto depth = state.total_depth(); + + Units::Milimetres expected{0.0}; + + EXPECT_DOUBLE_EQ(WE.value,expected.value); + constexpr auto MM_PER_M = 1000.0; + EXPECT_DOUBLE_EQ(depth.value * MM_PER_M,expected.value); +}; + +TEST_F(StateTest,ConstructWithIceOnly) +{ + Ice i = ice(); + State state(&p,i); + + auto IWE = i.water_equivalent(); + auto state_WE = state.total_water_equiv(); + + constexpr auto M_PER_MM = 1e-3; + const auto water_density = PhysConst::water_reference_density(); + auto expected_depth = IWE.value * water_density + / PhysConst::rho_ice() * M_PER_MM; + auto depth = state.total_depth(); + + EXPECT_DOUBLE_EQ(i.water_equivalent().value, + state.total_water_equiv().value); + EXPECT_DOUBLE_EQ(depth.value,expected_depth); + +}; + +TEST_F(StateTest,ConstructWithIceAndFirn) +{ + auto i = ice(); + auto f = firn(); + + State s(f,i); // No p, not needed by s + + Units::Milimetres init = i.water_equivalent(); + init += f.water_equivalent(); + + constexpr auto M_PER_MM = 1e-3; + const auto water_density = PhysConst::water_reference_density(); + auto expected_depth = i.water_equivalent().value * water_density + / PhysConst::rho_ice() * M_PER_MM; + auto& layers = f.get_layers(); + expected_depth += std::accumulate(layers.begin(),layers.end(), + 0.0, + [](double acc,Layer layer) { + acc += layer.height().value; + return acc; + }); + auto depth = s.total_depth(); + + + EXPECT_DOUBLE_EQ(init.value, + s.total_water_equiv().value); + EXPECT_DOUBLE_EQ(expected_depth,depth.value); +}; + +/* + * DailyMeltTest: Testing of the functions in the Glacier::Melt namespace which contain the + * melt logic. Used by the Model class directly. + */ + +class DailyMeltTest : public ::testing::Test +{ + std::array WE {{ + {200.0}, + {500.0}, + {250.0}, + {325.0} + }}; +protected: + bool system_set = false; + Params p = test_default_params(); + LayeredFirn firn{&p}; + void set_firn() + { + for (auto we : WE) + { + firn.accumulate(we); + } + + double sum = std::accumulate(WE.begin(),WE.end(),0.0, + [](double acc, const Units::Milimetres mm) + { + return acc + mm.value; + }); + EXPECT_NEAR(sum,firn.water_equivalent().value,MARGIN); + + auto& layers = firn.get_layers(); + + for (auto layer : layers) + { + EXPECT_LT(layer.density().value,p.firn_to_ice_density); + }; + + }; + + Ice ice{&p}; + void set_ice() + { + Units::Milimetres temp{400.0}; + ice.accumulate(temp); + }; + + Units::Milimetres swe{0.0}; + void set_swe() + { + swe.value = 450.0; + } + + Units::Milimetres melt_energy{0.0}; + void SetEnergySmall() + { + double firn_WE = std::accumulate(WE.begin(),WE.end(),0.0, + [](double acc,Units::Milimetres mm) { + acc += mm.value; + return acc;}); + melt_energy += Units::Milimetres{firn_WE/2.0}; + }; + + void SetEnergyBig() + { + double firn_WE = std::accumulate(WE.begin(),WE.end(),0.0, + [](double acc,Units::Milimetres mm) { + acc += mm.value; + return acc;}); + melt_energy = Units::Milimetres{firn_WE*2.0}; + //melt_energy += value; + }; + + void NoGlacier() { + }; + + void FullGlacier() { + set_swe(); + set_firn(); + set_ice(); + }; + + void no_firn() + { + set_swe(); + set_ice(); + }; + + void only_snow() + { + set_swe(); + }; + + void no_snow() + { + set_firn(); + set_ice(); + }; + + void only_firn() + { + set_firn(); + }; + + void no_ice() + { + set_swe(); + set_firn(); + }; + + void only_ice() + { + set_ice(); + }; + + template + void test_scenario(GlacierState f, + Melt::MeltScenario withoutEnergy, + Melt::MeltScenario withEnergy) + { + f(); + + State s(firn,ice); + + EXPECT_EQ(Melt::get_scenario(s,melt_energy,swe), + withoutEnergy); + + SetEnergySmall(); + + + melt_energy += Units::Milimetres{melt_energy}; + + EXPECT_EQ(Melt::get_scenario(s,melt_energy,swe), + withEnergy); + }; + + template + void test_scenario(GlacierState f, + Melt::MeltScenario withoutEnergy) + { + test_scenario(f,withoutEnergy,withoutEnergy); + }; +}; + +TEST_F(DailyMeltTest,ZeroWperM2_To_0mm) +{ + const Units::Watts_per_m2 input{0.0}; + using namespace Melt; + EXPECT_EQ(convert_to_mass(input,p).value,0.0); +}; + +TEST_F(DailyMeltTest,PositiveWperM2_To_mm) +{ + const Units::Watts_per_m2 input{1234.0}; + using namespace Melt; + using namespace PhysConst; + + auto output = Units::Watts_per_m2{input.value / + ( 0.95 * water_reference_density() * Lf() ) * p.seconds_per_step * 1000.0 }; + + EXPECT_EQ(convert_to_mass(input,p).value,output.value); +}; + + +TEST_F(DailyMeltTest,FullGlacierIsNoMelt) +{ + test_scenario([this]() { FullGlacier(); }, + Melt::MeltScenario::NoMelt); +}; + +TEST_F(DailyMeltTest,SWEandFirnIsNoMelt) +{ + test_scenario([this]() { no_ice(); }, + Melt::MeltScenario::NoMelt); +}; + +TEST_F(DailyMeltTest,OnlySnowIsNoMelt) +{ + test_scenario([this]() { only_snow(); }, + Melt::MeltScenario::NoMelt); +}; + +TEST_F(DailyMeltTest,NoGlacierIsNoMelt) +{ + test_scenario([this]() { NoGlacier(); }, + Melt::MeltScenario::NoMelt); +}; + +TEST_F(DailyMeltTest,NoSnowGlacier) +{ + test_scenario([this]() { no_snow(); }, + Melt::MeltScenario::NoMelt, + Melt::MeltScenario::FirnMelt); + + SetEnergyBig(); + + State s(firn,ice); + + EXPECT_EQ(Melt::get_scenario(s,melt_energy,swe), + Melt::MeltScenario::FirnAndIceMelt); +}; + +TEST_F(DailyMeltTest,OnlyIce) +{ + test_scenario([this]() { only_ice(); }, + Melt::MeltScenario::NoMelt, + Melt::MeltScenario::IceMelt); +}; + +TEST_F(DailyMeltTest,SnowAndIceNoMelt) +{ + test_scenario([this]() { no_firn(); }, + Melt::MeltScenario::NoMelt); +}; + +TEST_F(DailyMeltTest,OnlyFirn) +{ + test_scenario([this]() { only_firn(); }, + Melt::MeltScenario::NoMelt, + Melt::MeltScenario::FirnMelt); +}; + +TEST_F(DailyMeltTest,NoMelt) +{ + auto scenario = Melt::MeltScenario::NoMelt; + + FullGlacier(); + + State s(firn,ice); + + melt_energy = Units::Milimetres{1000.0}; + auto result = Melt::compute_melt(scenario,s,melt_energy); + + EXPECT_FALSE(result); +}; + +TEST_F(DailyMeltTest,FirnMelt) +{ + auto scenario = Melt::MeltScenario::FirnMelt; + + no_snow(); + + State s(firn,ice); + auto init_FWE = s.firn.water_equivalent(); + auto melt_energy = Units::Milimetres{firn.water_equivalent().value / 2.0}; + auto result = Melt::compute_melt(scenario,s,melt_energy); + auto final_FWE = s.firn.water_equivalent(); + + EXPECT_DOUBLE_EQ(result->melt.value,melt_energy.value); + EXPECT_DOUBLE_EQ(final_FWE.value,init_FWE.value - melt_energy.value); +}; + +TEST_F(DailyMeltTest,FirnAndIceMelt) +{ + auto scenario = Melt::MeltScenario::FirnAndIceMelt; + + no_snow(); + + State s(firn,ice); + auto init_IWE = s.ice.water_equivalent(); + constexpr auto ice_melt_fraction = 0.5; + + auto melt_energy = Units::Milimetres{firn.water_equivalent().value +init_IWE.value * ice_melt_fraction}; + auto result = Melt::compute_melt(scenario,s,melt_energy); + auto final_FWE = s.firn.water_equivalent(); + auto& layers = s.firn.get_layers(); + auto final_IWE = s.ice.water_equivalent(); + + EXPECT_DOUBLE_EQ(result->melt.value,melt_energy.value); + EXPECT_DOUBLE_EQ(final_FWE.value,0.0); + EXPECT_EQ(layers.size(),0); + EXPECT_DOUBLE_EQ(final_IWE.value,init_IWE.value * ice_melt_fraction); + +}; + +TEST_F(DailyMeltTest,FirnAndIceMeltAll) +{ + auto scenario = Melt::MeltScenario::FirnAndIceMelt; + + no_snow(); + + State s(firn,ice); + auto init_IWE = s.ice.water_equivalent(); + constexpr auto ice_melt_fraction = 1.5; + + auto melt_energy = Units::Milimetres{firn.water_equivalent().value +init_IWE.value * ice_melt_fraction}; + auto result = Melt::compute_melt(scenario,s,melt_energy); + auto final_FWE = s.firn.water_equivalent(); + auto& layers = s.firn.get_layers(); + auto final_IWE = s.ice.water_equivalent(); + + EXPECT_DOUBLE_EQ(result->melt.value,melt_energy.value - init_IWE.value * 0.5); + EXPECT_DOUBLE_EQ(final_FWE.value,0.0); + EXPECT_EQ(layers.size(),0); + EXPECT_DOUBLE_EQ(final_IWE.value,0.0); +}; + +TEST_F(DailyMeltTest,IceMelt) +{ + auto scenario = Melt::MeltScenario::IceMelt; + + only_ice(); + + State s(firn,ice); + + auto init_IWE = s.ice.water_equivalent(); + + auto melt_energy = Units::Milimetres{s.ice.water_equivalent().value * 0.75}; + auto result = Melt::compute_melt(scenario,s,melt_energy); + + auto final_IWE = s.ice.water_equivalent(); + + EXPECT_EQ(result->melt.value,melt_energy.value); + EXPECT_EQ(final_IWE.value,init_IWE.value - melt_energy.value); + +}; +/* + * The Updater is responsible for the infrequent updating part + * of the glacier season. The following processes occur: + * + * 1. Converts SWE to firn. + * 2. Converts firn to ice. + * 3. Increases the density of firn. + * + * The following is a test for the Updater class, constructed only + * as needed on the stack, obtaining references to firn and ice + * objects, as well as the swe as Units::Milimeter instance. + * + * The best approach is to have a member function for each: + * + * 1. SWE_to_firn() + * 2. firn_to_ice() + * 3. densify_firn() + * + */ +class YearlyUpdaterTest : public ::testing::Test +{ +protected: + Params p = test_default_params(); + Ice ice; + + YearlyUpdaterTest() : ice(&p) {}; + + // Construct fresh firn, no layers high density + LayeredFirn fresh_firn() { + std::vector WE {{ + {200.0}, + {500.0}, + {250.0}, + {325.0} + }}; + + LayeredFirn firn(&p); + + for (auto we : WE) + { + firn.accumulate(we); + } + + double sum = std::accumulate(WE.begin(),WE.end(),0.0, + [](double acc, const Units::Milimetres mm) + { + return acc + mm.value; + }); + EXPECT_NEAR(sum,firn.water_equivalent().value,MARGIN); + + auto& layers = firn.get_layers(); + + for (auto layer : layers) + { + EXPECT_LT(layer.density().value,p.firn_to_ice_density); + }; + + return firn; + }; + + //construct old firn, several (but not all) layers with + //high density + LayeredFirn old_firn() + { + struct Local + { + Layer::Height height; + Layer::Density density; + }; + std::array local{{ + {{4.0},{250.0}},{{3.2},{420.0}},{{2.5},{750.0}}, + {{1.2},{p.firn_to_ice_density + 20.0}}, + {{0.6},{p.firn_to_ice_density + 45.0}}, + {{0.3},{p.firn_to_ice_density + 100.0}} + }}; + + std::deque layers; + + for (auto l : local) + { + Layer layer(l.height,l.density); + layers.push_front(layer); + } + + LayeredFirn firn_local(&p,layers); + + auto& layers_ref = firn_local.get_layers(); + EXPECT_EQ(layers_ref.size(),6); + + size_t count_over_critical = 0; + for (auto layer : layers_ref) + { + if (layer.density().value > p.firn_to_ice_density) + ++count_over_critical; + } + EXPECT_EQ(count_over_critical,3); + + return firn_local; + }; + + void some_ice() + { + Units::Milimetres ice_init_WE{6e3}; + ice.accumulate(ice_init_WE); + }; + +}; + +TEST_F(YearlyUpdaterTest,NoSWEFreshFirnUnchanged) +{ + LayeredFirn firn = fresh_firn(); + Units::Milimetres swe{0.0}; + auto swe_copy = swe; + auto layers_init = firn.get_layers(); + + Updates::swe_to_firn(p,firn,swe); + + auto layers_after = firn.get_layers(); + + EXPECT_EQ(layers_init.size(),layers_after.size()); + EXPECT_EQ(layers_init.back().water_equivalent() + ,layers_after.back().water_equivalent()); + EXPECT_EQ(swe_copy.value,swe.value); +}; + +TEST_F(YearlyUpdaterTest,SWEupdatesFirn) +{ + LayeredFirn firn = fresh_firn(); + Units::Milimetres swe{150.0}; + auto swe_copy = swe; + const size_t init_firn_layers = + firn.get_layers().size(); + + Updates::swe_to_firn(p,firn,swe); + + auto& layers = firn.get_layers(); + + EXPECT_EQ(layers.size(),init_firn_layers+1); + + EXPECT_NEAR(layers.back().water_equivalent().value, + swe_copy.value,MARGIN); + + EXPECT_EQ(swe.value,0.0); +}; + +TEST_F(YearlyUpdaterTest,FreshFirnNotConvertedToIce) +{ + LayeredFirn firn = fresh_firn(); + Units::Milimetres swe{0.0}; + + //LayeredFirn firn2{Params()}; + const size_t init_firn_layers = + firn.get_layers().size(); + const Units::Milimetres init_ice = ice.water_equivalent(); + + Updates::firn_to_ice(p,firn,ice); + + EXPECT_EQ(firn.get_layers().size(),init_firn_layers); + + EXPECT_EQ(Units::Milimetres{0.0},ice.water_equivalent()); + +}; + +TEST_F(YearlyUpdaterTest,OldFirnConvertsFrontToIce) +{ + LayeredFirn firn = old_firn(); + some_ice(); + Units::Milimetres swe{0.0}; + + auto init_layers = firn.get_layers(); + const size_t init_firn_layers = + init_layers.size(); + const Units::Milimetres init_ice = ice.water_equivalent(); + + Updates::firn_to_ice(p,firn,ice); + + auto& new_layers = firn.get_layers(); + + auto it_new = new_layers.rbegin(); + auto it_init = init_layers.rbegin(); + + while (it_new != new_layers.rend()) + { + EXPECT_EQ(it_new->water_equivalent().value, + it_init->water_equivalent().value); + ++it_new; + ++it_init; + } + EXPECT_TRUE(it_init == init_layers.rend()-1); + + EXPECT_EQ(init_layers.back().water_equivalent(), + new_layers.back().water_equivalent()); + + EXPECT_EQ(init_firn_layers-1,new_layers.size()); + + EXPECT_NE(init_layers.front().water_equivalent(), + new_layers.front().water_equivalent()); + + EXPECT_EQ(new_layers.front().water_equivalent(), + init_layers[1].water_equivalent()); + + EXPECT_EQ(ice.water_equivalent().value, + init_ice.value + + init_layers.front().water_equivalent().value); +}; + +TEST_F(YearlyUpdaterTest,FreshFirnUnchangedIceAndFirn) +{ + LayeredFirn firn = fresh_firn(); + some_ice(); + + LayeredFirn init_firn = firn; + const Units::Milimetres init_ice = ice.water_equivalent(); + + Updates::firn_to_ice(p,firn,ice); + + EXPECT_EQ(init_firn.water_equivalent().value, + firn.water_equivalent().value); + + EXPECT_EQ(init_ice.value,ice.water_equivalent().value); +}; + +TEST_F(YearlyUpdaterTest,ManyIceUpdates) +{ + LayeredFirn firn = old_firn(); + const Units::Milimetres init_ice = ice.water_equivalent(); + + auto& layers = firn.get_layers(); + size_t init_layer_num = layers.size(); + auto init_layers = layers; + Units::Milimetres sum{init_ice}; + size_t iter_tracker = 0; + while (iter_tracker < init_layer_num) + { + if (layers.front().water_equivalent().value < p.firn_to_ice_density) + break; + ++iter_tracker; + sum += layers.front().water_equivalent(); + Updates::firn_to_ice(p,firn,ice); + EXPECT_EQ(sum.value,ice.water_equivalent().value); + + } + +}; + +class TestData +{ + State _state; + + // Inputs (set before calling execute_impl) + Units::Milimetres _swe{0.0}; + Units::Watts_per_m2 _melt_energy{0.0}; + Units::Celsius _glacier_temperature{273.15}; + bool _update_now{false}; + + // Outputs (recorded by execute_impl) + double _glacier_water_equivalent{0.0}; + double _total_depth{0.0}; + double _firn_melt{0.0}; + double _ice_melt{0.0}; + +public: + explicit TestData(const Params* p) : _state(p) { check_params(p);} + TestData(const Params* p, const LayeredFirn& f) + : _state(p, f) {check_params(p);} // State(Params*, LayeredFirn&) + TestData(const Params* p, const Ice& i) + : _state(p, i) {check_params(p);} + TestData(const LayeredFirn& f, const Ice& i) + : _state(f, i) {} + + void check_params(const Params* p) + { + EXPECT_EQ(p->seconds_per_step,3600) << "See hard-coded value in set_melt_energy"; + }; + + // ── Setters for test configuration ── + void set_swe(double v) { _swe.value = v; } + void set_melt_energy(double v) { + _melt_energy.value = v; } + void set_glacier_temperature(double v) { _glacier_temperature.value = v; } + void set_update_now(bool v) { _update_now = v; } + + // ── GlacierData concept interface ── + const Units::Milimetres swe() { return _swe; } + const Units::Watts_per_m2 melt_energy() { return _melt_energy; } + const Units::Celsius glacier_temperature() { return _glacier_temperature; } + bool update_now() { return _update_now; } + + void glacier_water_equivalent(const double v) { _glacier_water_equivalent = v; } + void total_depth(const double v) { _total_depth = v; } + void firn_melt(const double v) { _firn_melt = v; } + void ice_melt(const double v) { _ice_melt = v; } + + // ── Extra methods used by execute_impl ── + State& get_state() { return _state; } + const State& get_state() const { return _state; } + Units::Milimetres snowmelt() { return Units::Milimetres{0.0}; } // unused by logic + Units::Milimetres rainfall() { return Units::Milimetres{0.0}; } // unused by logic + + // ── Getters for test assertions ── + double recorded_glacier_water_equivalent() const { return _glacier_water_equivalent; } + double recorded_total_depth() const { return _total_depth; } + double recorded_firn_melt() const { return _firn_melt; } + double recorded_ice_melt() const { return _ice_melt; } +}; + +static_assert(GlacierData); + +class ModelTest : public ::testing::Test +{ +protected: + Params p; + Model model; + + // Helper: run model on data + void run(TestData& d) { model.execute_impl(d); } + + // Helper: create firn with known WE + LayeredFirn make_firn(double we_mm) + { + LayeredFirn f(&p); + if (we_mm > 0.0) + f.accumulate(Units::Milimetres{we_mm}); + return f; + } + + // Helper: create ice with known WE + Ice make_ice(double we_mm) + { + return Ice(&p, Units::Milimetres{we_mm}); + } + + double convert_mm_to_wperm2(const double v,const Params& p) + { + using namespace PhysConst; + // 0.95 -> themal quality factor + // mm/step * J/kg * kg/m^3 * step/s * m/mm = W/m^2 + // This is a very coupled test, with the 0.95 + // TODO improve this testing interface + return v * 0.95 * Lf() * water_reference_density() /( p.seconds_per_step * 1000.0); + }; +}; + +TEST_F(ModelTest, SwePositive_EmptyGlacier_NoMelt) +{ + TestData d(&p); + d.set_swe(100.0); + d.set_melt_energy(convert_mm_to_wperm2(500.0,p)); // would melt if SWE were 0 + + run(d); + + EXPECT_DOUBLE_EQ(d.recorded_firn_melt(), 0.0); + EXPECT_DOUBLE_EQ(d.recorded_ice_melt(), 0.0); +} + +TEST_F(ModelTest, SwePositive_WithFirn_NoMelt) +{ + auto firn = make_firn(200.0); + TestData d(firn, make_ice(0.0)); + d.set_swe(50.0); + d.set_melt_energy(convert_mm_to_wperm2(1000.0,p)); + + run(d); + + EXPECT_DOUBLE_EQ(d.recorded_firn_melt(), 0.0); + EXPECT_DOUBLE_EQ(d.recorded_ice_melt(), 0.0); +} + +TEST_F(ModelTest, SwePositive_WithIce_NoMelt) +{ + TestData d(&p, make_ice(500.0)); + d.set_swe(10.0); + d.set_melt_energy(convert_mm_to_wperm2(200.0,p)); + + run(d); + + EXPECT_DOUBLE_EQ(d.recorded_firn_melt(), 0.0); + EXPECT_DOUBLE_EQ(d.recorded_ice_melt(), 0.0); +} + +TEST_F(ModelTest, SwePositive_WithFirnAndIce_NoMelt) +{ + auto firn = make_firn(300.0); + TestData d(firn, make_ice(400.0)); + d.set_swe(1.0); + d.set_melt_energy(convert_mm_to_wperm2(9999.0,p)); + + run(d); + + EXPECT_DOUBLE_EQ(d.recorded_firn_melt(), 0.0); + EXPECT_DOUBLE_EQ(d.recorded_ice_melt(), 0.0); +} + +TEST_F(ModelTest, ZeroMeltEnergy_EmptyGlacier_NoMelt) +{ + TestData d(&p); + d.set_swe(0.0); + d.set_melt_energy(convert_mm_to_wperm2(0.0,p)); + + run(d); + + EXPECT_DOUBLE_EQ(d.recorded_firn_melt(), 0.0); + EXPECT_DOUBLE_EQ(d.recorded_ice_melt(), 0.0); +} + +TEST_F(ModelTest, ZeroMeltEnergy_WithFirn_NoMelt) +{ + auto firn = make_firn(200.0); + TestData d(firn, make_ice(0.0)); + d.set_swe(0.0); + d.set_melt_energy(convert_mm_to_wperm2(0.0,p)); + + run(d); + + EXPECT_DOUBLE_EQ(d.recorded_firn_melt(), 0.0); + EXPECT_DOUBLE_EQ(d.recorded_ice_melt(), 0.0); +} + +TEST_F(ModelTest, ZeroMeltEnergy_WithIce_NoMelt) +{ + TestData d(&p, make_ice(500.0)); + d.set_swe(0.0); + d.set_melt_energy(convert_mm_to_wperm2(0.0,p)); + + run(d); + + EXPECT_DOUBLE_EQ(d.recorded_firn_melt(), 0.0); + EXPECT_DOUBLE_EQ(d.recorded_ice_melt(), 0.0); +} + +TEST_F(ModelTest, FirnOnlyMelt_PartialMelt) +{ + // Firn only, no ice → FirnMelt scenario + auto firn = make_firn(500.0); + TestData d(firn, make_ice(0.0)); + d.set_swe(0.0); + d.set_melt_energy(convert_mm_to_wperm2(100.0,p)); // less than firn WE + d.set_glacier_temperature(0.0); + + run(d); + + EXPECT_GT(d.recorded_firn_melt(), 0.0); + EXPECT_DOUBLE_EQ(d.recorded_ice_melt(), 0.0); + // Firn WE should decrease + EXPECT_LT(d.get_state().firn.water_equivalent().value, 500.0); +} + +TEST_F(ModelTest, IceOnlyMelt) +{ + // Ice only, no firn → IceMelt scenario + TestData d(&p, make_ice(500.0)); + d.set_swe(0.0); + d.set_melt_energy(convert_mm_to_wperm2(100.0,p)); + d.set_glacier_temperature(0.0); + + run(d); + + EXPECT_DOUBLE_EQ(d.recorded_firn_melt(), 0.0); + EXPECT_GT(d.recorded_ice_melt(), 0.0); + EXPECT_LT(d.get_state().ice.water_equivalent().value, 500.0); +} + +TEST_F(ModelTest, FirnAndIceMelt_MeltExceedsFirn) +{ + // Both firn and ice present, melt energy exceeds firn WE + // → FirnAndIceMelt scenario + // Possibly failing this test because of thermal melt factor?! + auto firn = make_firn(50.0); // small firn + auto ice = make_ice(500.0); + TestData d(firn, ice); + d.set_swe(0.0); + d.set_melt_energy(convert_mm_to_wperm2(200.0,p)); // well exceeds firn + d.set_glacier_temperature(0.0); + + double initial_firn_we = d.get_state().firn.water_equivalent().value; + double initial_ice_we = d.get_state().ice.water_equivalent().value; + + run(d); + + // In FirnAndIceMelt, firn_melt should equal initial firn WE (all firn melted) + EXPECT_NEAR(d.recorded_firn_melt(), initial_firn_we, MARGIN); + // Ice melt should be positive + EXPECT_GT(d.recorded_ice_melt(), 0.0); + // Ice melt = initial_ice_we - remaining ice + double expected_ice_melt = initial_ice_we - d.get_state().ice.water_equivalent().value; + EXPECT_NEAR(d.recorded_ice_melt(), expected_ice_melt, MARGIN); +} + +TEST_F(ModelTest, FirnAndIce_MeltDoesNotExceedFirn) +{ + // Both present but melt energy < firn WE → FirnMelt only + auto firn = make_firn(500.0); + auto ice = make_ice(500.0); + TestData d(firn, ice); + d.set_swe(0.0); + d.set_melt_energy(convert_mm_to_wperm2(50.0,p)); // less than firn + d.set_glacier_temperature(0.0); + + run(d); + + EXPECT_GT(d.recorded_firn_melt(), 0.0); + EXPECT_DOUBLE_EQ(d.recorded_ice_melt(), 0.0); +} + +TEST_F(ModelTest, UpdateNowTrue_SweConvertedToFirn) +{ + TestData d(&p); + d.set_swe(300.0); + d.set_melt_energy(convert_mm_to_wperm2(0.0,p)); + d.set_update_now(true); + + EXPECT_EQ(d.get_state().firn.get_layers().size(), 0u); + + run(d); + + // SWE should have been accumulated as a new firn layer + EXPECT_GT(d.get_state().firn.get_layers().size(), 0u); + EXPECT_GT(d.get_state().firn.water_equivalent().value, 0.0); +} + +TEST_F(ModelTest, UpdateNowFalse_SweNotConvertedToFirn) +{ + TestData d(&p); + d.set_swe(300.0); + d.set_melt_energy(convert_mm_to_wperm2(0.0,p)); + d.set_update_now(false); + + run(d); + + // No firn conversion should happen + EXPECT_EQ(d.get_state().firn.get_layers().size(), 0u); +} + +TEST_F(ModelTest, UpdateNowTrue_DenseFirnConvertsToIce) +{ + // Manually create a firn layer with density above firn_to_ice threshold + Layer::Height h{1.0}; + Layer::Density rho{p.firn_to_ice_density + 10.0}; // above threshold + Layer dense_layer(h, rho); + + LayeredFirn firn(&p, std::vector{dense_layer}); + TestData d(firn, Ice(&p)); + d.set_swe(0.0); + d.set_melt_energy(convert_mm_to_wperm2(0.0,p)); + d.set_update_now(true); + + double initial_firn_we = d.get_state().firn.water_equivalent().value; + EXPECT_GT(initial_firn_we, 0.0); + EXPECT_DOUBLE_EQ(d.get_state().ice.water_equivalent().value, 0.0); + + run(d); + + // Dense layer should have been converted to ice + EXPECT_GT(d.get_state().ice.water_equivalent().value, 0.0); +} + +TEST_F(ModelTest, UpdateNowFalse_DenseFirnDoesNotConvertToIce) +{ + Layer::Height h{1.0}; + Layer::Density rho{p.firn_to_ice_density + 10.0}; + Layer dense_layer(h, rho); + + LayeredFirn firn(&p, std::vector{dense_layer}); + TestData d(firn, Ice(&p)); + d.set_swe(0.0); + d.set_melt_energy(convert_mm_to_wperm2(0.0,p)); + d.set_update_now(false); + + run(d); + + // No conversion + EXPECT_DOUBLE_EQ(d.get_state().ice.water_equivalent().value, 0.0); + EXPECT_GT(d.get_state().firn.water_equivalent().value, 0.0); +} + +TEST_F(ModelTest, OutputsMatchStateAfterExecution_EmptyGlacier) +{ + TestData d(&p); + d.set_swe(0.0); + d.set_melt_energy(convert_mm_to_wperm2(0.0,p)); + + run(d); + + EXPECT_DOUBLE_EQ(d.recorded_glacier_water_equivalent(), + d.get_state().total_water_equiv().value); + EXPECT_DOUBLE_EQ(d.recorded_total_depth(), + d.get_state().total_depth().value); +} + +TEST_F(ModelTest, OutputsMatchStateAfterExecution_WithMelt) +{ + auto firn = make_firn(300.0); + TestData d(firn, make_ice(200.0)); + d.set_swe(0.0); + d.set_melt_energy(convert_mm_to_wperm2(50.0,p)); + d.set_glacier_temperature(0.0); + + run(d); + + EXPECT_DOUBLE_EQ(d.recorded_glacier_water_equivalent(), + d.get_state().total_water_equiv().value); + EXPECT_DOUBLE_EQ(d.recorded_total_depth(), + d.get_state().total_depth().value); +} + +TEST_F(ModelTest, OutputsMatchStateAfterExecution_WithUpdate) +{ + TestData d(&p); + d.set_swe(500.0); + d.set_melt_energy(convert_mm_to_wperm2(0.0,p)); + d.set_update_now(true); + + run(d); + + EXPECT_DOUBLE_EQ(d.recorded_glacier_water_equivalent(), + d.get_state().total_water_equiv().value); + EXPECT_DOUBLE_EQ(d.recorded_total_depth(), + d.get_state().total_depth().value); +} + +TEST_F(ModelTest, MeltExceedsTotalGlacier_FullDepletion) +{ + auto firn = make_firn(50.0); + auto ice = make_ice(50.0); + TestData d(firn, ice); + d.set_swe(0.0); + d.set_melt_energy(convert_mm_to_wperm2(99999.0,p)); // massive melt + d.set_glacier_temperature(0.0); + + run(d); + + // Everything should be melted + EXPECT_NEAR(d.get_state().total_water_equiv().value, 0.0, 1e-6); + EXPECT_DOUBLE_EQ(d.recorded_glacier_water_equivalent(), + d.get_state().total_water_equiv().value); +} + +TEST_F(ModelTest, IceMeltThenUpdate_NoFirnConversionOnEmptyFirn) +{ + // Ice-only glacier, melt some ice, then update + TestData d(&p, make_ice(500.0)); + d.set_swe(0.0); + d.set_melt_energy(convert_mm_to_wperm2(50.0,p)); + d.set_glacier_temperature(0.0); + d.set_update_now(true); + + run(d); + + EXPECT_GT(d.recorded_ice_melt(), 0.0); + EXPECT_DOUBLE_EQ(d.recorded_firn_melt(), 0.0); + // update_now with no SWE and no dense firn → no state change from updates + EXPECT_EQ(d.get_state().firn.get_layers().size(), 0u); +} diff --git a/src/tests/submoduletests/test_katabatic_melt_energy.cpp b/src/tests/submoduletests/test_katabatic_melt_energy.cpp new file mode 100644 index 00000000..d1d2ef08 --- /dev/null +++ b/src/tests/submoduletests/test_katabatic_melt_energy.cpp @@ -0,0 +1,320 @@ +#include +#include "katabatic_melt_energy.hpp" +#include "PhysConst.h" +#include "Atmosphere.h" +#include +#include + +using namespace katabatic_melt_energy; +namespace Units = PhysConst::units; + +inline Params test_params() +{ + Params p; + + p.prandtl = 5.0; + p.k = 4e-4; + p.k2 = 1.0; + p.seconds_per_step = 3600; + + return p; +}; + +TEST(ParamsKatabaticTest,ParamsSetByPhysConst) +{ + Params p; + + + EXPECT_EQ(p.heat_capacity_air,PhysConst::Cp()); + EXPECT_EQ(p.molecular_wt_ratio,PhysConst::M()); + EXPECT_EQ(p.g,PhysConst::g()); +} + +class AirDensityTest : public ::testing::Test +{ +protected: + struct Data + { + Units::Pa pressure; + Units::Pa vapour_pressure; + Units::Kelvin air_temperature; + Units::DensitySI expected_air_density; + + Data(const double p, const double e, const double T, const double rho) : + pressure{p}, vapour_pressure{e}, air_temperature{T}, expected_air_density{rho} {}; + }; + + std::array data = {{ + {1e5,1e4,273.15,1.227141074763899}, + {1.1e5,1e3,300.2,1.2720886528749895} + }}; +}; +TEST_F(AirDensityTest,Analytical) +{ + for (auto d : data) + { + auto result = Atmosphere::air_density(d.pressure,d.air_temperature,d.vapour_pressure); + EXPECT_DOUBLE_EQ(result.value,d.expected_air_density.value); + }; + +}; + +TEST_F(AirDensityTest,ZeroPressureThrows) +{ + EXPECT_THROW(Atmosphere::air_density( Units::Pa{0.0}, + Units::Kelvin{273.15}, + Units::Pa{1e7}),std::runtime_error); +}; + +class BulkCoefficientTest : public ::testing::Test +{ +protected: + Params p = test_params(); + + struct Input + { + double C; + double gamma; + double prandtl; + double glacier_temperature; + double solution; + }; + std::array inputs{{ + {-10,0.005,5.0,273.15,0.004791841256996453}, + {10.0,0.001,1.0,250.3,-0.025028948481418418} + }}; +}; + +TEST_F(BulkCoefficientTest,AnalyticalBulkCoefficent) +{ + for (auto input : inputs) + { + input.solution = p.k * std::pow(p.k2,2.0) * input.C * std::sqrt(p.g / (input.glacier_temperature * input.gamma * input.prandtl)); + p.prandtl = input.prandtl; + auto result = bulk_coefficient( p, + Units::TempDiff{input.C}, + Units::LapseRateSI{input.gamma}, + Units::Kelvin{input.glacier_temperature}); + + EXPECT_EQ(result.value,input.solution); + } + +}; + +TEST_F(BulkCoefficientTest,CertainZerosThrow) +{ + /* + * inputs[0] -> gamma = 0, divide by zero + * inputs[1] -> Prandtl number = 0, divide by zero + * inputs[2] -> glacier temperature = 0, divide by zero + */ + std::array inputs {{ + {1.0,0.0,1.0,200.0,0.0}, + {1.0,1.0,0.0,200.0,0.0}, + {1.0,1.0,1.0,0.0,0.0} + }}; + + for (auto input : inputs) + { + p.prandtl = 0.0; + EXPECT_THROW(bulk_coefficient( p, + Units::Kelvin{input.C}, + Units::LapseRateSI{input.gamma}, + Units::Kelvin{input.glacier_temperature});, + std::runtime_error); + } + +}; + +class KatabaticHeatTest : public ::testing::Test +{ +protected: + Params p = test_params(); + +}; + +TEST_F(KatabaticHeatTest,AnalyticalSensibleHeat) +{ + const Units::m_per_s coefficient{1e3}; + const Units::Celsius deficit{10.0}; + const Units::DensitySI air_density{1.5}; + const auto expected = air_density.value * p.heat_capacity_air * coefficient.value * deficit.value; + + auto result = sensible_heat(p,coefficient,deficit,air_density); + + EXPECT_DOUBLE_EQ(result.value,expected); + +}; + +TEST_F(KatabaticHeatTest,AnalyticlLatentHeat) +{ + const Units::m_per_s coefficient{45.0}; + const Units::DensitySI air_density{1.102}; + const Units::Celsius glacier_temperature{-1.3}; + const Units::Pa vapour_pressure_deficit{1e-3}; + const auto expected = p.molecular_wt_ratio * air_density.value * PhysConst::Lv(glacier_temperature) * coefficient.value * vapour_pressure_deficit.value; + + auto result = latent_heat(p,coefficient,vapour_pressure_deficit,glacier_temperature,air_density); + + EXPECT_EQ(result.value,expected); +}; + + +// Simple mock data class that satisfies KatabaticData concept +class MockKatabaticData { +public: + // Store the values that will be passed to latent_heat and sensible_heat + double last_latent_heat_value{0}; + double last_sensible_heat_value{0}; + + // Required getters + Units::Kelvin glacier_temperature() const { return glacier_temp; } + Units::Kelvin air_temperature() const { return air_temp; } + Units::Pa air_pressure() const { return air_press; } + Units::Pa vapour_pressure() const { return vapour_press; } + Units::Pa vapour_pressure_surface() const { return vapour_press_surface; } + Units::LapseRateSI lapse_rate() const { return lapse_rate_val; } + + // Methods called by Model + void latent_heat(double value) { last_latent_heat_value = value; } + void sensible_heat(double value) { last_sensible_heat_value = value; } + + // Setters for test configuration + void set_glacier_temperature(Units::Kelvin v) { glacier_temp = v; } + void set_air_temperature(Units::Kelvin v) { air_temp = v; } + void set_air_pressure(Units::Pa v) { air_press = v; } + void set_vapour_pressure(Units::Pa v) { vapour_press = v; } + void set_vapour_pressure_surface(Units::Pa v) { vapour_press_surface = v; } + void set_lapse_rate(Units::LapseRateSI v) { lapse_rate_val = v; } + +private: + Units::Kelvin glacier_temp{273.15}; + Units::Kelvin air_temp{273.15}; + Units::Pa air_press{101325}; + Units::Pa vapour_press{1000}; + Units::Pa vapour_press_surface{500}; + Units::LapseRateSI lapse_rate_val{0.0065}; +}; + +class KatabaticModelTest : public ::testing::Test { +protected: + void SetUp() override { + // Set up default parameters + params.seconds_per_step = 3600; // 1 hour + } + + Params params; +}; + +TEST_F(KatabaticModelTest, ExecuteImplComputesCorrectSensibleHeat) { + MockKatabaticData data; + data.set_air_temperature(Units::Kelvin{275.15}); // 2°C + data.set_glacier_temperature(Units::Kelvin{273.15}); // 0°C + data.set_air_pressure(Units::Pa{101325}); + data.set_vapour_pressure(Units::Pa{1000}); + data.set_lapse_rate(Units::LapseRateSI{0.0065}); + + Model model; + + // Manually compute expected values using the static functions + auto rho_air = Atmosphere::air_density( + data.air_pressure(), + data.glacier_temperature(), + data.vapour_pressure() + ); + + Units::Kelvin T_deficit{ + data.air_temperature().value - data.glacier_temperature().value + }; + + auto K = bulk_coefficient( + params, + T_deficit, + data.lapse_rate(), + data.glacier_temperature() + ); + + auto expected_sensible = sensible_heat(params, K, T_deficit, rho_air); + + // Execute the model + model.execute_impl(data); + + // Verify + EXPECT_DOUBLE_EQ(data.last_sensible_heat_value, expected_sensible.value); +} + +TEST_F(KatabaticModelTest, ExecuteImplComputesCorrectLatentHeat) { + MockKatabaticData data; + data.set_air_temperature(Units::Kelvin{275.15}); + data.set_glacier_temperature(Units::Kelvin{273.15}); + data.set_air_pressure(Units::Pa{101325}); + data.set_vapour_pressure(Units::Pa{1000}); + data.set_vapour_pressure_surface(Units::Pa{500}); + data.set_lapse_rate(Units::LapseRateSI{0.0065}); + + Model model; + + // Manually compute expected values + auto rho_air = Atmosphere::air_density( + data.air_pressure(), + data.glacier_temperature(), + data.vapour_pressure() + ); + + Units::TempDiff T_deficit{ + data.air_temperature().value - data.glacier_temperature().value + }; + + auto K = bulk_coefficient( + params, + T_deficit, + data.lapse_rate(), + data.glacier_temperature() + ); + + Units::Pa vapour_pressure_deficit{ + data.vapour_pressure().value - data.vapour_pressure_surface().value + }; + + auto expected_latent = latent_heat( + params, + K, + vapour_pressure_deficit, + data.glacier_temperature(), + rho_air + ); + + // Execute + model.execute_impl(data); + + // Verify + EXPECT_DOUBLE_EQ(data.last_latent_heat_value, expected_latent.value); +} + +TEST_F(KatabaticModelTest, ExecuteImplWithZeroTemperatureGradient) { + MockKatabaticData data; + data.set_air_temperature(Units::Kelvin{273.15}); // Same as glacier + data.set_glacier_temperature(Units::Kelvin{273.15}); + data.set_air_pressure(Units::Pa{101325}); + data.set_vapour_pressure(Units::Pa{1000}); + + Model model; + model.execute_impl(data); + + // With zero temperature gradient, sensible heat should be zero + EXPECT_DOUBLE_EQ(data.last_sensible_heat_value, 0.0); +} + +TEST_F(KatabaticModelTest, ExecuteImplWithZeroVapourPressureDeficit) { + MockKatabaticData data; + data.set_vapour_pressure(Units::Pa{1000}); + data.set_vapour_pressure_surface(Units::Pa{1000}); // Equal pressures + data.set_air_temperature(Units::Kelvin{275.15}); + data.set_glacier_temperature(Units::Kelvin{273.15}); + + Model model; + model.execute_impl(data); + + // With zero vapour pressure deficit, latent heat should be zero + EXPECT_DOUBLE_EQ(data.last_latent_heat_value, 0.0); +} diff --git a/src/tests/submoduletests/test_melt_routing.cpp b/src/tests/submoduletests/test_melt_routing.cpp new file mode 100644 index 00000000..9b23c180 --- /dev/null +++ b/src/tests/submoduletests/test_melt_routing.cpp @@ -0,0 +1,309 @@ +#include "gtest/gtest.h" +#include +#include +#include "melt_routing_glacier.hpp" +#include "PhysConst.h" + +namespace Units = PhysConst::units; +using namespace GlacierRouting; + +inline Params test_params() +{ + Params p; + + p.chain_length = 10u; + p.k = 4e3; + p.seconds_per_step = 3600u; + + return p; +}; + +class LinearReservoirTest : public ::testing::Test +{ +protected: + void run(const double k, const size_t dt, std::vector inputs) { + LinearReservoir res(k,dt); + auto c0 = dt / (2*k + dt); + auto c1 = (2*k - dt) / (2*k + dt); + + auto last_input = 0.0; + auto last_output = 0.0; + + for (auto input : inputs) + { + auto result = res.step(input); + auto expected_result = c0 * (input + last_input) + + c1 * last_output; + EXPECT_DOUBLE_EQ(result,expected_result); + + last_input = input; + last_output = expected_result; + } + } +}; + +TEST_F(LinearReservoirTest,ZeroKInIsOut) +{ + LinearReservoir res(0.0,1.0); + + std::vector Vec{105.0,250.0,1234.111,24591.0}; + for (auto v : Vec) + { + auto output = res.step(v); + EXPECT_EQ(output,v); + } +}; + +TEST_F(LinearReservoirTest,NonZeroKLessThanCriticalThrows) +{ + const auto dt = 1000.0; + const auto k = dt / 2 - 0.01; + auto func = [=]() { + LinearReservoir res(k,dt); + }; + + EXPECT_THROW(func(),std::runtime_error); +}; + +TEST_F(LinearReservoirTest,StepTestSimple) +{ + const auto dt = 1500.0; + const auto k = 3050.0; + + std::vector inputs{1.0,2.5,1.3,8.2}; + + run(k,dt,inputs); +}; + +TEST_F(LinearReservoirTest,StepTestLong) +{ + size_t N = 1000; + std::vector inputs; + auto dt = 3600u; + auto k = 5000.0; + auto amplitude = 250.5; + for (size_t i = 0; i < N; ++i) + { + auto val = amplitude * + ( std::sin( amplitude* (i * dt)) + 1 ); + + inputs.emplace_back(val); + } + + run(k,dt,inputs); + +} + +class GlacierReservoirTest : public ::testing::Test +{ +protected: + Params p = test_params(); + +}; + +TEST_F(GlacierReservoirTest,ZeroUntilChainLengthSteps) +{ + GlacierReservoir res(&p); + + double initial_input = 150.0; + + res.step(initial_input); + + for (size_t i = 1; i < p.chain_length -1; ++i) + { + auto result = res.step(0.0); + EXPECT_EQ(result,0.0); + } + + auto result = res.step(0.0); + EXPECT_GT(result,0.0); +}; + +TEST_F(GlacierReservoirTest,ChainLengthOneSameAsLinearReservoir) +{ + p.chain_length = 1; + LinearReservoir lin_res(p.k,p.seconds_per_step); + GlacierReservoir gla_res(&p); + + double input = 100.0; + + auto out_lin = lin_res.step(input); + auto out_gla = gla_res.step(input); + + EXPECT_EQ(out_lin,out_gla); + + // Since they are equal, this just makes sure they aren't both + // zero + EXPECT_GT(out_lin + out_gla,0.0); +}; + +TEST_F(GlacierReservoirTest,LongChainSameAsLinearReservoirButDelayed) +{ + p.chain_length = 75u; + LinearReservoir lin_res(p.k,p.seconds_per_step); + GlacierReservoir gla_res(&p); + + size_t N = 500; + std::vector inputs; + auto dt = p.seconds_per_step; + auto amplitude = 250.5; + for (size_t i = 0; i < N; ++i) + { + auto val = amplitude * + ( std::sin( amplitude* (i * dt)) + 1 ); + + inputs.emplace_back(val); + } + + for (size_t i = 0; i < N; ++i) + { + auto out_lin = 0.0; + + if ( i + 1 >= p.chain_length) + out_lin = lin_res.step(inputs[i - p.chain_length + 1]); + + auto out_gla = gla_res.step(inputs[i]); + + if (i < p.chain_length - 1) + { + EXPECT_DOUBLE_EQ(out_gla,0.0); + } + else + ASSERT_DOUBLE_EQ(out_gla,out_lin); + }; +}; + +TEST(GlacierRoutingStateTest,Construct) +{ + Params p = test_params(); + + State s(&p); +}; + +class Data +{ + State s; + struct Inout + { + double melt = 0.0; + double delayed = 0.0; + }; + Inout snow; + Inout firn; + Inout ice; + double _total_delayed; +public: + State& get_state() { return s; } + Data(const Params* _p) : s(_p) {}; + const Units::Milimetres snowmelt() + { return Units::Milimetres{snow.melt}; } + const Units::Milimetres firnmelt() + { return Units::Milimetres{firn.melt}; } + const Units::Milimetres icemelt() + { return Units::Milimetres{ice.melt}; } + + void snowmelt(const double T) { snow.melt = T; }; + void firnmelt(const double T) { firn.melt = T; }; + void icemelt(const double T) { ice.melt = T; }; + + + void snowmelt_delayed(const double T) + { snow.delayed = T; }; + void firnmelt_delayed(const double T) + { firn.delayed = T; }; + void icemelt_delayed(const double T) + { ice.delayed = T; }; + void total_delayed(const double T) + { _total_delayed = T; }; + + const double snowmelt_delayed() + { return snow.delayed; } + const double firnmelt_delayed() + { return firn.delayed; } + const double icemelt_delayed() + { return ice.delayed; } +}; + +class GlacierRoutingModelTest : public ::testing::Test +{ +protected: + Model model; + Params& p = model.get_params(); + std::vector snowmelts; + std::vector firnmelts; + std::vector icemelts; + const size_t N = 500u; + void set_up_inputs(const double amplitude,std::vector& melts) + { + auto dt = p.seconds_per_step; + for (size_t i = 0; i < N; ++i) + { + auto val = amplitude * + ( std::sin( amplitude* (i * dt)) + 1 ); + + melts.emplace_back(val); + } + }; + + void SetUp() override { + p = test_params(); + p.chain_length = 10u; + }; +}; + +TEST_F(GlacierRoutingModelTest,ExecuteImpl) +{ + set_up_inputs(150.0,snowmelts); + set_up_inputs(45.3,firnmelts); + set_up_inputs(12.3,icemelts); + Data d(&p); + + struct Res + { + LinearReservoir snow; + LinearReservoir firn; + LinearReservoir ice; + Res(const double k, const size_t dt) + : snow(k,dt), firn(k,dt), ice(k,dt) {}; + } res(p.k,p.seconds_per_step); + + for (size_t i = 0; i < N; ++i) + { + d.snowmelt(snowmelts[i]); + d.firnmelt(firnmelts[i]); + d.icemelt(icemelts[i]); + + model.execute(d); + + auto out_snow = 0.0; + auto out_firn = 0.0; + auto out_ice = 0.0; + + if ( i + 1 >= p.chain_length) + { + out_snow = res.snow.step(snowmelts[i - p.chain_length + 1]); + out_firn = res.firn.step(firnmelts[i - p.chain_length + 1]); + out_ice = res.ice.step(icemelts[i - p.chain_length + 1]); + } + + + if (i < p.chain_length - 1) + { + EXPECT_DOUBLE_EQ(d.snowmelt_delayed(),0.0); + EXPECT_DOUBLE_EQ(d.firnmelt_delayed(),0.0); + EXPECT_DOUBLE_EQ(d.icemelt_delayed(),0.0); + } + else + { + EXPECT_DOUBLE_EQ(d.snowmelt_delayed(),out_snow); + EXPECT_DOUBLE_EQ(d.firnmelt_delayed(),out_firn); + EXPECT_NE(d.snowmelt_delayed(),d.firnmelt_delayed()); + ASSERT_DOUBLE_EQ(d.icemelt_delayed(),out_ice); + } + + } +}; + + + + diff --git a/src/tests/submoduletests/test_submodule.cpp b/src/tests/submoduletests/test_submodule.cpp new file mode 100644 index 00000000..a2af35fd --- /dev/null +++ b/src/tests/submoduletests/test_submodule.cpp @@ -0,0 +1,99 @@ +#include +#include +#include "data_base.hpp" +#include "base_step.hpp" +#include "triangulation.hpp" + +template + +concept SubmoduleData = requires(T& t) +{ + {t.input1()} -> std::floating_point; + {t.input2()} -> std::floating_point; + + {t.output(std::declval())} -> std::same_as; +}; + +template +class submodule : public base_step,data> +{ +public: + submodule() {}; + ~submodule() {}; + + void execute_impl(data& d) const + { + auto input1 = d.input1(); + auto input2 = d.input2(); + + auto output = input1 * input2; + d.output(output); + } +}; + +struct Cache : public cache_base +{ + double input1 = default_value(); + double input2 = default_value(); + + double output = 0.0; +}; + +class data2 : public data_base +{ +public: + data2(mesh_elem face_in,boost::shared_ptr param,pt::ptree& cfg) + : data_base(face_in,param,cfg,true) {}; + + double input1() + { + update_value( + [this]() -> auto& { return cache_->input1;}, + []() { return 4.0; } + ); + + return cache_->input1; + }; + + double input2() + { + update_value( + [this]() -> auto& { return cache_->input2;}, + []() { return 2.5; } + ); + + return cache_->input1; + }; + + void output(const double t) + { + set_output([this]() -> auto& { return cache_->output; },t); + }; +}; + +class TestModule : public ::testing::Test +{ +protected: + TestModule() : d(face,param,cfg) {}; + mesh_elem face{nullptr}; + boost::shared_ptr param; + pt::ptree cfg; + + data2 d; + submodule submodule_instance; + void run() + { + submodule_instance.execute(d); + }; +}; + +TEST_F(TestModule,SimpleSubmoduleTest) +{ + const auto& cache = d.get_cache(); + EXPECT_FALSE(cache); + run(); + EXPECT_TRUE(cache); + EXPECT_EQ(cache->input1,4.0); + EXPECT_EQ(cache->input2,2.5); + EXPECT_DOUBLE_EQ(cache->output,cache->input1 * cache->input2); +}; diff --git a/src/tests/submoduletests/test_water_flux.cpp b/src/tests/submoduletests/test_water_flux.cpp new file mode 100644 index 00000000..4624e5b2 --- /dev/null +++ b/src/tests/submoduletests/test_water_flux.cpp @@ -0,0 +1,146 @@ +#include +#include "PhysConst.h" +#include "water_flux.hpp" + +TEST(MassFluxTest,ConstructFromWperMSquared) +{ + double Q = 1400.0; + water_flux m = water_flux::from_W_per_m_squared(Q); +}; + +TEST(MassFluxTest,ReturnsMMperDT) +{ + double Q = 1400.0; + PhysConst::units::Celsius T{0.0}; + water_flux m = water_flux::from_W_per_m_squared(Q,T); + + auto result = m.W_per_m_squared(T); + + EXPECT_EQ(Q,result); +}; + +TEST(MassFluxTest,WperMSquaredtoMMperDT) +{ + auto Q = 250.5; + PhysConst::units::Celsius T{10.0}; + water_flux m = water_flux::from_W_per_m_squared(Q,T); + + int dt = 3600; + double Q_massflux = m.mm_per_dt(dt); + + constexpr auto mm_per_m = 1000.0; + + // [Q] = J/m^2/s + // [Q_massflux] = mm/s + // [water density] = kg/m^3 + // [latent heat of vapourization] = J/kg + // [milimetre per metre] = mm/m + // [dt] = s/step + // + // J/m^2/s * m^3/kg * kg/J * mm/m * s/step = + // J*m/s * 1/kg * kg/J * mm/m * s/step = + // m/s * mm/m * s/step = + // mm/step + auto Q_expect = Q / (PhysConst::water_reference_density() * PhysConst::Lv(T)) * mm_per_m; + + Q_expect *= dt; + EXPECT_DOUBLE_EQ(Q_expect,Q_massflux); +}; + +TEST(MassFluxTest,GetBackWperMSquared) +{ + double Q = 1400.0; + water_flux m = water_flux::from_W_per_m_squared(Q); + + auto Q_expect = m.W_per_m_squared(); + + EXPECT_DOUBLE_EQ(Q,Q_expect); +}; + +TEST(MassFluxTest,ConstructFromMMperS) +{ + double val = 200.0; + water_flux m = water_flux<>::from_mm_per_s(val); +}; + +TEST(MassFluxTest,GetBackMMperS) +{ + const double val = 200.0; + water_flux m = water_flux::from_mm_per_s(val); + + auto result = m.mm_per_dt(1); + + EXPECT_DOUBLE_EQ(result,val); + + auto result2 = m.mm_per_s(); // alias for mm_per_dt + + //result = result2 is implicit + EXPECT_DOUBLE_EQ(result2,val); +}; + +TEST(MassFluxTest,ConstructFromMMperDT) +{ + auto val = 34.0; + auto dt = 1000; + water_flux m = water_flux<>::from_mm_per_dt(34.0,dt); + + auto result = m.mm_per_dt(dt); + + EXPECT_DOUBLE_EQ(result,val); + + auto result_per_s = m.mm_per_s(); + + EXPECT_DOUBLE_EQ(val,result_per_s * dt); + +}; + +TEST(MassFluxTest,ConstructFromMperS) +{ + auto input = 550.0; + + water_flux m = water_flux<>::from_m_per_s(input); + + auto result = m.mm_per_s(); + + constexpr auto MM_PER_M = 1000.0; + EXPECT_DOUBLE_EQ(result,input*MM_PER_M); +}; + +TEST(MassFluxTest,GetBackMperS) +{ + auto input = 550.0; + + water_flux m = water_flux<>::from_mm_per_s(input); + + auto result = m.m_per_s(); + + constexpr auto M_PER_MM = 1 / 1000.0; + EXPECT_DOUBLE_EQ(result,input * M_PER_MM); +}; + +TEST(MassFluxTest,PlusOperatorOverload) +{ + auto value1 = 50.4; + auto m1 = water_flux<>::from_mm_per_s(value1); + auto value2 = 250.3; + auto m2 = water_flux<>::from_mm_per_s(value2); + + auto m3 = m1+m2; + auto m4 = m2+m1; + EXPECT_EQ(m3.mm_per_s(),value1+value2); + EXPECT_EQ(m4.mm_per_s(),value1+value2); +}; + +TEST(MassFluxTest,MinusOperatorOverload) +{ + auto value1 = 50.4; + auto m1 = water_flux<>::from_mm_per_s(value1); + auto value2 = 250.3; + auto m2 = water_flux<>::from_mm_per_s(value2); + + auto m3 = m1-m2; + auto m4 = m2-m1; + + EXPECT_EQ(m3.mm_per_s(),value1 - value2); + EXPECT_EQ(m4.mm_per_s(), value2 - value1); +}; diff --git a/src/tests/test_units.cpp b/src/tests/test_units.cpp new file mode 100644 index 00000000..aac208b1 --- /dev/null +++ b/src/tests/test_units.cpp @@ -0,0 +1,94 @@ +#include "gtest/gtest.h" +#include "PhysConst.h" + +using namespace PhysConst::units; + +/* + * Construct units of each type + * + * Add them together, substract them, no impl mulitply + * + * > < == <= >= += -= + */ +template +concept EqualityComparable = + requires (A a, B b) { + { a == b } -> std::convertible_to; + }; + +template +concept ThreeWayComparable = + requires (A a, B b) { + { a <=> b }; + }; + + +class UnitsTest : public ::testing::Test +{ +protected: + struct Example : public base {}; + class Other : public base {}; + + + static_assert(!EqualityComparable); + static_assert(!ThreeWayComparable); +}; + +TEST_F(UnitsTest,ConstructEmptyIsZero) +{ + Example example; + + EXPECT_EQ(example.value,0.0); +}; + +TEST_F(UnitsTest,ConstructInitializerBracesFromDouble) +{ + constexpr double init = 4.0; + Example exam{init}; + + EXPECT_EQ(exam.value,init); +}; + +TEST_F(UnitsTest,PlusEqual) +{ + Example e1{5.0}, e2{6.5}; + + e1 += e2; + + EXPECT_EQ(e1.value,11.5); +}; + +TEST_F(UnitsTest,MinusEqual) +{ + Example e1{9.0}, e2{2.5}; + + e1 -= e2; + + EXPECT_EQ(e1.value,6.5); +}; + +TEST_F(UnitsTest,GreaterThan) +{ + Example e1{10.0}, e2{5.5}; + + EXPECT_TRUE(e1 > e2); + EXPECT_TRUE(e1 >= e2); + EXPECT_TRUE(e1 >= e1); +}; + +TEST_F(UnitsTest,LessThan) +{ + Example e1{1000.0}, e2{450.0}; + + EXPECT_TRUE(e1 < e2); + EXPECT_TRUE(e1 <= e2); + EXPECT_TRUE(e1 <= e1); + +}; + +TEST_F(UnitsTest,Equal) +{ + Example e1{3486.35}; + + EXPECT_TRUE(e1 == e1); +}; diff --git a/test_data/unit-test-example-data/unit-test_example.json b/test_data/unit-test-example-data/unit-test_example.json new file mode 100644 index 00000000..cead13cc --- /dev/null +++ b/test_data/unit-test-example-data/unit-test_example.json @@ -0,0 +1,212 @@ +{ + "option": { + "per_triangle_timeseries": "true", + "ui": "false", + "debug_level": "debug", + "point_mode": { + }, + "prj_name": "UNITTEST" + }, + "modules": [ + "Infil_All", + "soil_module", + "Evapotranspiration_All", + "point_mode" + ], + "remove_depency": { + "Simple_Canopy": "snobal", + "Richard_albedo": "snobal", + "scale_wind_vert": "snobal", + + "Simple_Canopy": "Lehning_snowpack", + "Richard_albedo": "Lehning_snowpack", + "scale_wind_vert": "Lehning_snowpack", + + "scale_wind_vert": "FSM", + "snow_slide": "snobal", + "PBSM3D": "snobal", + "WindNinja": "snobal", + "WindNinja": "FSM", + "WindNinja": "Lehning_snowpack", + + "Simple_Canopy": "FSM", + "scale_wind_vert": "FSM", + "snow_slide": "FSM", + "PBSM3D": "FSM", + "WindNinja": "FSM", + + "Evapotranspiration_All" : "soil_module", + "Infil_All" : "soil_module" + }, + "config": { + "point_mode": { + "provide": { + "t": "false", + "U_2m_above_srf": "false", + "iswr": "true", + "unit_test_Donovan": "true" + } + }, + "Evapotranspiration_All": { + "alpha": "1.26", + "stomatal_resistance_min": "25", + "Frac_to_ground": 0.1, + "wind_measurement_height": "10" + }, + + "Infil_All": { + "infDays": "6", + "min_swe_to_freeze": "25.0", + "major": "5", + "AllowPrior_inf": "true", + "thaw_type": "0", + "lens_temp": "-10.0" + }, + + "soil_module": { + "Ksaturated_rechr": "0.00176", + "Ksaturated_lower": "6.9E-6", + "Ksaturated_groundwater": "6.9E-6", + "Ksaturated_organic": "6.9E-6", + "Trigthrhld": "100.0", + "depth_per_layer": "0.5", + "theta_default": "0.5", + "uniform_conductivities": "true", + "dry_soil_k": "2.5", + "sat_soil_frozen_k": "1.98", + "sat_soil_thaw_k": "1.67", + "number_XG_layers": "10", + "moisture_content_min_per_layer": "0.0001", + "update_k_behind_front_freeze": "true", + "update_k_behind_front_thaw": "true", + "k_update": "1", + "Johansen_conductivity": "false" + }, + + "Richard_albedo": { + "min_swe_refresh": "10", + "init_albedo_snow": "0.8", + "a2": "3.6e5" + }, + "snobal": { + "param_snow_compaction": "1" + }, + "Thornton_p": { + "apply_cosine_correction": "true" + }, + "p_no_lapse": { + "apply_cosine_correction": "true" + }, + "snow_slide": { + "use_vertical_snow": "true" + }, + "PBSM3D": { + "nLayer": "10", + "smooth_coeff": "6500", + "rouault_diffusion_coef": "false", + "min_sd_trans": "0.1", + "enable_veg": "true", + "snow_diffusion_const": "0.3", + "settling_velocity": "0.5", + "do_fixed_settling": "true", + "do_sublimation": "true", + "use_tanh_fetch": "true", + "use_R94_lambda": "false", + "use_PomLi_probability": "false", + "use_subgrid_topo": "false", + "use_subgrid_topo_V2": "false", + "z0_ustar_coupling": "false", + "debug_output": "false" + }, + "fetchr": { + "incl_veg": "false", + "I": "0.06" + }, + "WindNinja": { + "Max_spdup": "3.5", + "ninja_recirc": "true", + "compute_Sx": "false", + "Sx_crit": "20.0", + "L_avg": "1000" + }, + "iswr": { + "already_cosine_corrected": "true" + } + }, + "meshes": { + + "mesh": "unit-test_example.mesh", + "parameters": { + "file": "unit-test_example.param" + } + + + }, + "parameter_mapping": { + "landcover": { + "1": { + "is_water": "false" + }, + "0": { + "is_water": "true" + } + }, + "soils": { + "soil_type": { + "0": "sand", + "1": "loam", + "2": "clay", + "3": "loamsand", + "4": "sandloam", + "5": "siltloam", + "6": "saclloam", + "7": "clayloam", + "8": "siclloam", + "9": "sandclay", + "10": "siltclay", + "99": "pavement" + }, + "soil_texture": { + "0": "coarse_over_coarse", + "1": "medium_over_medium", + "2": "medium_over_fine", + "3": "fine_over_fine", + "4": "soil_over_bedrock" + }, + "soil_groundcover": { + "0" : "bare_soil", + "1" : "row_crop", + "2" : "poor_pasture", + "3" : "small_grains", + "4" : "good_pasture", + "5" : "forested" + } + + } + }, + "output": { + "UpperClearing" : { + "longitude": -115.175362, + "latitude": 50.956547, + "file": "output_test.txt", + "type": "timeseries" + } + + }, + "forcing": + { + "UTC_offset":11, + "num_stations_to_use": "1", + "interpolant": "nearest", + + "UpperClearing": + { + "file": "combined.csv", + "longitude": -115.175362, + "latitude": 50.95654, + "elevation": 1845 + } + + + } +} diff --git a/test_data/unit-test-example-data/unit-test_example.mesh b/test_data/unit-test-example-data/unit-test_example.mesh new file mode 100644 index 00000000..cb6013e3 --- /dev/null +++ b/test_data/unit-test-example-data/unit-test_example.mesh @@ -0,0 +1,40 @@ +{ + "mesh": { + "vertex": [ + [ + 488479.5, + 6713528.0, + 1000.0 + ], + [ + 490527.5, + 6713528.0, + 1000.0 + ], + [ + 490527.5, + 6711480.0, + 1000.0 + ] + ], + "nvertex": 3, + "neigh": [ + [ + -1, + -1, + -1 + ] + ], + "elem": [ + [ + 1, + 0, + 2 + ] + ], + "is_geographic": 0, + "proj4": "+proj=utm +zone=8 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ", + "UTM_zone": 8, + "nelem": 1 + } +} diff --git a/test_data/unit-test-example-data/unit-test_example.param b/test_data/unit-test-example-data/unit-test_example.param new file mode 100644 index 00000000..bac9c318 --- /dev/null +++ b/test_data/unit-test-example-data/unit-test_example.param @@ -0,0 +1,116 @@ +{ + "soil_storage": [ + 375.0 + ], + "soil_rechr_storage": [ + 125.0 + ], + "soil_storage_max": [ + 750 + ], + "soil_rechr_max": [ + 250 + ], + "detention_snow_init": [ + 0.0 + ], + "detention_organic_init": [ + 0.0 + ], + "ground_water_storage": [ + 187.0 + ], + "ground_water_max": [ + 375.0 + ], + "excess_to_ssr": [ + 1 + ], + "detention_snow_max": [ + 0.0 + ], + "detention_organic_max": [ + 0.0 + ], + "depression_max": [ + 1.0 + ], + "depression_storage": [ + 0.0 + ], + "soil_type": [ + 5 + ], + "LAI": [ + 0.3 + ], + "LAImax": [ + 0.5 + ], + "CanopyHeight": [ + 1.5 + ], + "soil_texture": [ + 0 + ], + "soil_groundcover": [ + 5 + ], + "soil_groundcover_ET": [ + 2 + ], + "elevation": [ + 1845 + ], + "landcover": [ + 1 + ], + "local_slope": [ + 5 + ], + "PSD_K_estimator": [ + 2.55 + ], + "PSD_K_organic": [ + 0.252 + ], + "Ksaturated_rechr": [ + 0.00176 + ], + "Ksaturated_lower": [ + 6.9E-6 + ], + "Ksaturated_ground_water": [ + 6.9E-6 + ], + "Ksaturated_organic": [ + 6.9E-6 + ], + "Ksaturated_snow": [ + 6.9E-6 + ], + "soil_index": [ + 3.3 + ], + "snow_grain_diameter": [ + 3.0 + ], + "soil_depth": [ + 1.5 + ], + "allow_runoff_from_infiltration": [ + 1 + ], + "perma_frost_depth": [ + 2 + ], + "init_freeze_front_depth": [ + 0 + ], + "init_thaw_front_depth": [ + 0 + ], + "landcover": [ + 1 + ] +}