Commit ff486b35 authored by Hannes Bauser's avatar Hannes Bauser Committed by Santiago Ospina De Los Ríos
Browse files

Update scaling adapter unit test

* Port test to Google Test.
* Add meta-ini config file to build grid in test fixture.
* Add Python script for creating scaling datasets before testing.
* Add test case for non-cubic dataset.
parent 834e62b4
......@@ -2,6 +2,10 @@
## Unreleased
### Added
* Extended Miller scaling adapter with additional porosity scaling !122
### Changed
* Linerar solver for finite volumes changed from `AMG_4_DG` to
......
......@@ -20,6 +20,7 @@ add_custom_target(example_tests
# create the mapping datafile as preparation for all tests
add_custom_target(prepare_testing
COMMAND ${DUNE_PYTHON_VIRTUALENV_EXECUTABLE} ${PROJECT_SOURCE_DIR}/test/maps/create_param_maps.py ${PROJECT_SOURCE_DIR}/test/maps/cell_ids.h5
COMMAND ${DUNE_PYTHON_VIRTUALENV_EXECUTABLE} ${PROJECT_SOURCE_DIR}/dune/dorie/test/create_scaling.py ${PROJECT_BINARY_DIR}/dune/dorie/test/scaling.h5
)
add_dependencies(system_tests prepare_testing)
add_dependencies(unit_tests prepare_testing)
......
......@@ -126,18 +126,19 @@ Richards Parameterizations
~~~~~~~~~~~~~~~~~~~~~~~~~~
As the *soil porosity* :math:`\phi \, [-]` and the *residual water content*
:math:`\theta_s \, [-]` vary between soil types, it is convenient to define the
:math:`\theta_r \, [-]` vary between soil types, it is convenient to define the
*soil water saturation* :math:`\Theta \, [-]` by
.. math::
\Theta = \frac{\theta_w - \theta_r}{\phi - \theta_r},
\Theta (\theta_w) = \frac{\theta_w - \theta_r}{\phi - \theta_r},
\quad \Theta \in \left[ 0, 1 \right],
where :math:`\theta_w \, [-]` is the volumetric soil water content.
One typically assumes that the entire pore space can be filled up with water,
hence we set the *saturated water content* :math:`\theta_s \, [-]` equal to the
porosity, :math:`\theta_s = \phi`.
This relation can be manipulated in certain :ref:`local scalings <man-scaling>`.
Mualem–van Genuchten
""""""""""""""""""""
......@@ -492,6 +493,7 @@ Effective Hydromechanic Dispersion Tensor for Isotropic Media
lambda_t:
# <Deff_type> parameters ...
.. _man-scaling:
Scalings
~~~~~~~~
......@@ -507,6 +509,7 @@ Every ``scaling_section`` has the following layout:
extensions: <sequence>
offset: <sequence>
The setting ``interpolation`` accepts any of the :doc:`implemented Interpolators <interpolators>`.
The values of ``extensions`` and ``offset`` are sequences containing the
coordinates in the respective spatial dimensions.
......@@ -545,6 +548,34 @@ Miller Scaling
scale_miller:
# ...
Miller Porosity Scaling
"""""""""""""""""""""""
Applies porosity scaling in addition to regular :ref:`Miller scaling <man-miller_scaling>`, violating its geometric similarity assumption.
.. math::
\Theta &= \Theta (h_m \cdot \xi_M) \\
K &= K (\Theta) \cdot \xi_M^2 \\
\phi &= \theta_s - \xi_\theta
* ``type: MilPor``
Scaling sections:
* ``scale_miller``: Scaling factor :math:`\xi_M` applied onto matric head and conductivity simultaneously.
* ``scale_porosity``: Scaling summand :math:`\xi_\theta` subtracted from the saturated conductivity.
YAML template:
.. code-block:: yaml
type: MilPor
data:
scale_miller:
# ...
scale_porosity:
# ...
.. _YAML: https://gettaurus.org/docs/YAMLTutorial/
.. _HDF5: https://www.h5py.org/
......@@ -72,7 +72,7 @@ which are applied to the parameterization functions as follows:
@f{align*}{
\Theta &= \Theta(h_m \cdot \xi_{h_m}) \\
\phi &= \theta_s + \xi_{\phi} \\
\phi &= \theta_s - \xi_{\phi} \\
K &= K(\Theta) \cdot \xi_K^2
@f}
......
......@@ -5,6 +5,8 @@
#include <string>
#include <map>
#include <dune/common/float_cmp.hh>
namespace Dune {
namespace Dorie {
namespace Parameterization {
......@@ -123,9 +125,11 @@ public:
const RF& scale_por
) const
{
verify_porosity_scaling(scale_por);
return [this, &scale_por](const Saturation sat) {
return WaterContent{
sat.value * (_theta_s.value - _theta_r.value + scale_por)
sat.value * (_theta_s.value - _theta_r.value - scale_por)
+ _theta_r.value
};
};
......@@ -140,13 +144,15 @@ public:
const RF& scale_por
) const
{
verify_porosity_scaling(scale_por);
return [this, &scale_por](const WaterContent theta)
{
// Small value to avoid numerical issues
constexpr double eps_min = 1E-6;
const double sat = (theta.value - _theta_r.value)
/ (_theta_s.value + scale_por - _theta_r.value);
/ (_theta_s.value - scale_por - _theta_r.value);
return Saturation{std::clamp(sat, eps_min, 1.0)};
};
}
......@@ -189,6 +195,24 @@ public:
*/
virtual std::unique_ptr<Richards<Traits>> clone () const = 0;
private:
/// Check if the porosity scaling factor conforms to the parameterization
/** \param scaling The porosity scaling factor to check
*/
void verify_porosity_scaling (const RF scaling) const
{
if (FloatCmp::eq(scaling, 0.0))
return;
if (scaling >= _theta_s.value - _theta_r.value) {
DUNE_THROW(Dune::IOError, "Invalid porosity scaling "
"(negative water content).");
}
else if(_theta_s.value - scaling > 1.0){
DUNE_THROW(Dune::IOError, "Invalid porosity scaling "
"(water content exceeds 1).");
}
}
};
......
......@@ -160,6 +160,62 @@ public:
~MillerScalingAdapter () override = default;
};
/// ScalingAdapter implementing Miller scaling with additional porosity scaling.
/**
* For the Description of Miller scaling, see MillerScalingAdapter.
* In addition, the porosity is scaled through an additive scaling
* \f$ \xi_\theta \f$:
* \f[
* \Theta (\theta) =
* \frac{\theta - \theta_r}{\theta_s - \xi_\theta - \theta_r}
* \f]
*
* The MilPor scaling violates the geometric similarity assumption of Miller
* scaling. The choice of \f$ -\xi_\theta \f$ (instead of \f$ +\xi_\theta \f$)
* is due to a material with smaller length scale (Miller scaling factor
* \f$ \xi < 1 \f$) having a larger porosity.
*
* \ingroup RichardsParam
* \author Hannes Bauser
* \date 2018
*/
template<typename Traits>
class MilPorScalingAdapter : public ScalingAdapter<Traits>
{
private:
using Base = ScalingAdapter<Traits>;
using Domain = typename Base::Domain;
using GridView = typename Base::GridView;
using return_t = typename Base::return_t;
public:
/// Create this adapter and the internal interpolator.
explicit MilPorScalingAdapter (
const YAML::Node& scale_data_cfg,
const GridView& grid_view,
const std::shared_ptr<spdlog::logger> log
):
Base(scale_data_cfg, grid_view, log)
{
this->add_interpolator(scale_data_cfg["scale_miller"], "scale_miller");
this->add_interpolator(scale_data_cfg["scale_porosity"], "scale_porosity");
}
/// Export type of this scaling adapter
static inline std::string type = "MilPor";
/// Evaluate this interpolator
return_t evaluate (const Domain& pos) const override
{
const auto miller_scaling_value = this->_interpolators[0]->evaluate(pos);
const auto porosity_scaling_value = this->_interpolators[1]->evaluate(pos);
return return_t{miller_scaling_value, miller_scaling_value, porosity_scaling_value};
}
/// Delete this adapter
~MilPorScalingAdapter () override = default;
};
/// A ScalingAdapter that does not implement any scaling.
/** This adapter returns a default-initialized instance of ScalingFactors
* when evaluated.
......@@ -241,6 +297,11 @@ public:
grid_view,
log);
}
else if (type == MilPorScalingAdapter<Traits>::type) {
return std::make_shared<MilPorScalingAdapter<Traits>>(config,
grid_view,
log);
}
else if (to_lower(type) == DummyScalingAdapter<Traits>::type) {
return std::make_shared<DummyScalingAdapter<Traits>>(config,
grid_view,
......
message(STATUS "Handling unit tests")
# copy ALL the config files and stuff
file(COPY scaling.h5
scaling.yml
test-param-richards-scaled.yml
file(COPY test-param-richards-scaled.yml
test-param-richards-unscaled.yml
test-param-transport.yml
test-param-factory.yml
DESTINATION .)
# scaling adapter test
dorie_add_unit_test(SOURCES test-model-base.cc NAME test-model-base)
dorie_add_unit_test(SOURCES test-grid-function-container.cc NAME test-grid-function-container)
dorie_add_unit_test(SOURCES test-interpolators.cc NAME test-interpolators)
dorie_add_unit_test(SOURCES test-scaling-adapters.cc NAME test-scaling-adapters)
dorie_add_unit_test(SOURCES test-param-factory.cc NAME test-param-factory)
dorie_add_unit_test(SOURCES test-boundary-condition.cc
NAME test-boundary-condition)
# scaling adapter test
dorie_add_metaini_test(
UNIT_TEST CUSTOM_MAIN
SOURCE test-scaling-adapters.cc
BASENAME test-scaling-adapters
METAINI test-scaling-adapters.mini.in
CREATED_TARGETS sa_target
)
# boundary condition test
dorie_add_metaini_test(
UNIT_TEST CUSTOM_MAIN
......
import h5py
import numpy as np
import os
import argparse
def create_and_write_datasets(args):
zeros = np.zeros((2, 2))
ones = np.ones((2, 2))
twos = np.ones((2, 2)) * 2
shift = np.zeros((2, 1))
shift[1, ...] = 1
porosity = np.ones((2, 2)) * 0.1
f = h5py.File(args.path, "w")
f.create_dataset("zeros", data=zeros)
f.create_dataset("ones", data=ones)
f.create_dataset("twos", data=twos)
f.create_dataset("shift", data=shift)
f.create_dataset("porosity", data=porosity)
f.close()
def parse_cli():
parser = argparse.ArgumentParser(
description="Write parameter mapping datasets for all tests"
)
parser.add_argument(
"path", type=os.path.realpath, help="File path to write the data"
)
return parser.parse_args()
if __name__ == "__main__":
args = parse_cli()
create_and_write_datasets(args)
# macro grid on which the scaling adapters are evaluated
grid:
extensions: [1, 1]
offset: [0, 0]
scalings:
# Loads zeros as Miller scaling factors
scale1:
type: Miller
data:
scale_miller:
file: scaling.h5
dataset: twos
interpolation: nearest
extensions: [1, 1]
offset: [0, 0]
# Loads zeros as Miller scaling factors (automatic extension deduction)
scale2:
type: Miller
data:
scale_miller:
file: scaling.h5
dataset: twos
interpolation: nearest
# Loads Dummy scale adapter returning 1.0 as factors
scale3:
type: none
......@@ -19,7 +19,7 @@ file = none
type = {_initial_type}
quantity = {_quantity}
file = ${CMAKE_CURRENT_LIST_DIR}/scaling.h5
file = scaling.h5
dataset = twos
interpolation = nearest
......@@ -34,6 +34,6 @@ quantity = soluteConcentration
equation = -1+2*y
file = ${CMAKE_CURRENT_LIST_DIR}/scaling.h5
file = scaling.h5
dataset = twos
interpolation = nearest
#include <gtest/gtest.h>
#ifdef HAVE_CONFIG_H
#include "config.h"
#include "config.h"
#endif
#include <string>
#include <random>
#include <exception>
#include <random>
#include <string>
#include <vector>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/fvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/float_cmp.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/dorie/common/grid_creator.hh>
#include <dune/dorie/common/setup.hh>
......@@ -21,98 +18,203 @@
#include "test-parameterization.hh"
/// Compare a scaled and an unscaled parameterization.
/** Retrieve the actual scaling factors of the scaled parameterization
* and apply them to the functions of the unscaled one. Then compare
* the results.
*/
template<class FlowParameters, class Grid>
void test_scaled_parameterization (FlowParameters& param_scaled,
FlowParameters& param_unscaled,
std::shared_ptr<Grid> grid)
int argc = 0;
char** argv;
class ScaledParamTest : public ::testing::Test
{
// initialize RNG
std::mt19937 rng(1);
std::uniform_real_distribution<double> dist_head(-50.0, 0.0);
std::uniform_real_distribution<double> dist_sat(0.0, 1.0);
auto gv = grid->leafGridView();
for (auto&& cell : elements(gv))
{
// attach to element
param_scaled.bind(cell);
param_unscaled.bind(cell);
// retrieve actual scaling
using Scale = Dune::Dorie::Parameterization::ScalingFactors<Traits<2>::RF>;
auto scale = std::get<Scale>(param_scaled.cache());
// check results of scaled functions
namespace Cmp = Dune::FloatCmp;
// matric head -> saturation
auto sat_f_s = param_scaled.saturation_f();
auto sat_f_u = param_unscaled.saturation_f();
const auto scale_head = scale.scale_head;
const auto head = dist_head(rng);
assert(Cmp::eq(sat_f_s(head), sat_f_u(head * scale_head)));
// saturation -> conductivity
auto cond_f_s = param_scaled.conductivity_f();
auto cond_f_u = param_unscaled.conductivity_f();
const auto scale_cond = scale.scale_cond;
const auto sat = dist_sat(rng);
assert(Cmp::eq(cond_f_s(sat), cond_f_u(sat) * scale_cond * scale_cond));
}
protected:
using T = Traits<2>;
using Grid = typename T::Grid;
using Domain = typename T::Domain;
using Scaling = Dune::Dorie::Parameterization::ScalingFactors<typename T::RF>;
std::shared_ptr<spdlog::logger> log;
Dune::ParameterTree config;
std::shared_ptr<typename T::Grid> grid;
std::vector<int> index_map;
std::mt19937 rng;
std::uniform_real_distribution<double> dist_head;
std::uniform_real_distribution<double> dist_sat;
std::uniform_real_distribution<double> dist_wc;
void SetUp() override
{
auto [cfg, lg, helper] = Dune::Dorie::Setup::init(argc, argv);
log = lg;
config = cfg;
const auto log_level = config.get<std::string>("richards.output.logLevel");
log = Dune::Dorie::create_logger(
Dune::Dorie::log_richards, helper, spdlog::level::from_str(log_level));
// Create the grid
Dune::Dorie::GridCreator<Grid> gc(cfg, helper);
grid = gc.grid();
index_map = gc.element_index_map();
// Set up RNG
rng = std::mt19937(1);
dist_head = std::uniform_real_distribution<double>(-50.0, 0.0);
dist_sat = std::uniform_real_distribution<double>(0.0, 1.0);
// NOTE: Depends on parameter set and porosity scaling in test!
dist_wc = std::uniform_real_distribution<double>(0.03, 0.2);
}
void TearDown() override { spdlog::drop_all(); }
};
using Parameters = Dune::Dorie::FlowParameters<Traits<2>>;
template<class Grid>
Parameters
create_flow_parameters(const std::string& parameters_file,
Dune::ParameterTree config,
const std::shared_ptr<Grid> grid,
const std::vector<int>& index_map)
{
config["parameters.file"] = parameters_file;
return Parameters(config, grid, index_map);
}
int main (int argc, char** argv)
// --- Test cases --- //
TEST_F(ScaledParamTest, HeadSaturation)
{
try{
// Initialize ALL the things!
auto [inifile, log, helper] = Dune::Dorie::Setup::init(argc, argv);
// create the grid
Dune::Dorie::GridCreator<Traits<2>::Grid> gc(inifile, helper);
auto grid = gc.grid();
const auto index_map = gc.element_index_map();
inifile = Dune::Dorie::Setup::prep_ini_for_richards(inifile);
// create the Richards logger because FlowParameters expects it
const auto log_level = inifile.get<std::string>("output.logLevel");
log = Dune::Dorie::create_logger(Dune::Dorie::log_richards,
helper,
spdlog::level::from_str(log_level));
inifile["parameters.file"] = "test-param-richards-scaled.yml";
Dune::Dorie::FlowParameters<Traits<2>> param_scaled(inifile,
grid,
index_map);
// use unscaled params now
inifile["parameters.file"] = "test-param-richards-unscaled.yml";
Dune::Dorie::FlowParameters<Traits<2>> param_unscaled(inifile,
grid,
index_map);
test_scaled_parameterization(param_scaled, param_unscaled, grid);
return 0;
}
catch (Dune::Exception &e) {
std::cerr << "Dune reported error: " << e << std::endl;
return 1;
}
catch (std::exception &e) {
std::cerr << "Exception thrown: " << e.what() << std::endl;
return 1;
}
catch (...) {
std::cerr << "Unknown exception!" << std::endl;
return 1;
}
auto param_unscaled = create_flow_parameters(
"test-param-richards-unscaled.yml", config, grid, index_map);
auto param_scaled = create_flow_parameters(
"test-param-richards-scaled.yml", config, grid, index_map);
auto gv = grid->leafGridView();
for (auto&& cell : elements(gv))
{
// Attach to element
param_scaled.bind(cell);
param_unscaled.bind(cell);
// Retrieve stored scaling information for element
auto scale = std::get<Scaling>(param_scaled.cache());
// Retrieve parameterization functions
auto sat_f_s = param_scaled.saturation_f();
auto sat_f_u = param_unscaled.saturation_f();
const auto scale_head = scale.scale_head;
// Perform test
const auto head = dist_head(rng);
EXPECT_DOUBLE_EQ(sat_f_s(head), sat_f_u(head * scale_head));
// Retrieve inverse functions
auto head_f_s = param_scaled.matric_head_f();
auto head_f_u = param_unscaled.matric_head_f();
// Perform test
const auto sat = dist_sat(rng);
EXPECT_DOUBLE_EQ(head_f_s(sat), head_f_u(sat) / scale_head);
}
}
TEST_F(ScaledParamTest, SaturationToConductivity)
{
auto param_unscaled = create_flow_parameters(
"test-param-richards-unscaled.yml", config, grid, index_map);
auto param_scaled = create_flow_parameters(
"test-param-richards-scaled.yml", config, grid, index_map);
auto gv = grid->leafGridView();
for (auto&& cell : elements(gv))
{
// Attach to element
param_scaled.bind(cell);
param_unscaled.bind(cell);
// Retrieve stored scaling information for element
auto scale = std::get<Scaling>(param_scaled.cache());
// Retrieve parameterization functions
auto cond_f_s = param_scaled.conductivity_f();
auto cond_f_u = param_unscaled.conductivity_f();
const auto scale_cond = scale.scale_cond;
const auto sat = dist_sat(rng);
EXPECT_DOUBLE_EQ(cond_f_s(sat), cond_f_u(sat) * scale_cond * scale_cond);