Commit 6cbc598b authored by Santiago Ospina De Los Ríos's avatar Santiago Ospina De Los Ríos

Merge branch 'add-extendedMiller-ScalingAdapter' into 'master'

add extended miller scaling adapter

See merge request !122
parents 834e62b4 ff486b35
......@@ -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
......@@ -3,8 +3,6 @@ include ${CMAKE_BINARY_DIR}/doc/default_files/config.ini
__name = scaled_param
_asset_path = "${PROJECT_SOURCE_DIR}/test"
adaptivity.useAdaptivity = false
[grid]
initialLevel = 1
gridType = rectangular
......
......@@ -24,7 +24,7 @@ volumes:
tau: 0.0
scaling:
type: Miller
type: MilPor
data:
scale_miller:
file: scaling.h5
......@@ -32,3 +32,9 @@ scaling:
interpolation: nearest
extensions: [1, 1]
offset: [0, 0]
scale_porosity:
file: scaling.h5
dataset: porosity
interpolation: nearest
extensions: [1, 1]
offset: [0, 0]
#include <gtest/gtest.h>
#ifdef HAVE_CONFIG_H
#include "config.h"
#include "config.h"
#endif
#include <cassert>
#include <exception>
#include <vector>
#include <random>
#include <functional>
#include <random>
#include <string>
#include <vector>
#include <yaml-cpp/yaml.h>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/float_cmp.hh>
#include <dune/common/fvector.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/dorie/common/grid_creator.hh>
#include <dune/dorie/common/setup.hh>
#include <dune/dorie/model/richards/parameterization/scaling.hh>
/// Compare two scaling factor structs
template<typename T>
std::enable_if_t<std::is_arithmetic_v<T>, bool> operator==
(const Dune::Dorie::Parameterization::ScalingFactors<T>& a,
const Dune::Dorie::Parameterization::ScalingFactors<T>& b)
int argc = 0;
char** argv;
struct Traits
{
if (a.scale_cond == b.scale_cond
and a.scale_head == b. scale_head
and a.scale_por == b.scale_por)
static constexpr int dim = 2;
using RF = double;
using DF = double;
using Grid = Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<DF, dim>>;
using Domain = Dune::FieldVector<DF, dim>;
using GridView = typename Grid::LeafGridView;
};
/// Fixture storing scaling settings, grid, and sampling positions
class ScalingFixture : public ::testing::Test
{
protected:
using Grid = typename Traits::Grid;
using Domain = typename Traits::Domain;
using SAF = Dune::Dorie::Parameterization::ScalingAdapterFactory<Traits>;
using Scaling =
Dune::Dorie::Parameterization::ScalingFactors<typename Traits::RF>;
std::shared_ptr<spdlog::logger> log;
std::shared_ptr<typename Traits::Grid> grid;
YAML::Node scale_file;
/// Positions where the scaling adapter will be sampled
std::vector<Domain> pos_sample;
void SetUp() override
{
auto [cfg, lg, helper] = Dune::Dorie::Setup::init(argc, argv);
log = lg;
scale_file = YAML::LoadFile(cfg["_scaling_file_name"]);
// Create the grid
Dune::Dorie::GridCreator<Grid> gc(cfg, helper);
grid = gc.grid();
// Refine the grid so we sample more than the coarsest grid vertices
grid->globalRefine(1);
// Store the sampling positions
for (auto&& vertex : vertices(grid->leafGridView()))
{
return true;
pos_sample.emplace_back(vertex.geometry().center());
}
}
return false;
void TearDown() override { spdlog::drop_all(); }
};
// Inject into correct namespace for Google Test
namespace Dune::Dorie::Parameterization
{
/// Compare two scaling factor structs
template<typename T>
std::enable_if_t<std::is_arithmetic_v<T>, bool>
operator==(const ScalingFactors<T>& a, const ScalingFactors<T>& b)
{
using namespace Dune::FloatCmp;
if (eq(a.scale_cond, b.scale_cond) and eq(a.scale_head, b.scale_head)
and eq(a.scale_por, b.scale_por))
{
return true;
}
return false;
}
/// Traits struct to be plugged in as traits mock-up
template<int dimension>
struct Traits
} // namespace Dune::Dorie::Parameterization
/// Create a scaling adapter
template<class GV>
decltype(auto)
create_scaling(const YAML::Node& settings,
const GV& grid_view,
const std::shared_ptr<spdlog::logger> log)
{
static constexpr int dim = dimension;
using RF = double;
using DF = double;
using Grid = Dune::YaspGrid<dim,
Dune::EquidistantOffsetCoordinates<DF, dim>>;
using Domain = Dune::FieldVector<DF, dim>;
using GridView = typename Grid::LeafGridView;
};
using SAF = Dune::Dorie::Parameterization::ScalingAdapterFactory<Traits>;
const auto type = settings["type"].as<std::string>();
return SAF::create(type, settings["data"], grid_view, log);
}
/// Test a single scaling setting
template<int dim>
void test_scaling_setting (const YAML::Node& cfg,
const typename Traits<dim>::GridView& grid_view)
// --- Test cases --- //
TEST_F(ScalingFixture, DummyScaling)
{
// create the scaling
const auto type = cfg["type"].as<std::string>();
using SAF = Dune::Dorie::Parameterization::ScalingAdapterFactory<Traits<dim>>;
const auto log = Dune::Dorie::get_logger(Dune::Dorie::log_base);
auto scaling = SAF::create(type, cfg["data"], grid_view, log);
// initialize random positions
using Domain = typename Traits<dim>::Domain;
std::vector<Domain> positions;
for(auto&& vertex : vertices(grid_view)) {
positions.push_back(vertex.geometry().center());
}
auto scaling = create_scaling(scale_file["dummy"], grid->leafGridView(), log);
const Scaling reference{ 1.0, 1.0, 0.0 };
for (auto&& pos : pos_sample)
{
const auto value = scaling->evaluate(pos);
EXPECT_EQ(value, reference) << "Position : " << pos[0] << ", " << pos[1];
}
}
evaluate_scaling(type, scaling, positions);
TEST_F(ScalingFixture, MillerAutomaticExtension)
{
auto scaling = create_scaling(
scale_file["miller_automatic_extension"], grid->leafGridView(), log);
const Scaling reference{ 2.0, 2.0, 0.0 };
for (auto&& pos : pos_sample)
{
const auto value = scaling->evaluate(pos);
EXPECT_EQ(value, reference) << "Position : " << pos[0] << ", " << pos[1];
}
}
/// Evaluate the given scaling adapter for a set of positions
template<typename ScalingAdapter, typename Domain>
void evaluate_scaling (const std::string& type,
const std::shared_ptr<ScalingAdapter> scaling,
const std::vector<Domain>& positions)
TEST_F(ScalingFixture, MillerCustomExtension)
{
// create reference scaling
using Scaling = Dune::Dorie::Parameterization::ScalingFactors<double>;
Scaling reference;
if (type == "none") {
reference = Scaling{1.0, 1.0, 0.0};
}
else if (type == "Miller") {
reference = Scaling{2.0, 2.0, 0.0};
}
else {
throw std::runtime_error("Scaling type not supported "
"by this test");
}
auto scaling = create_scaling(
scale_file["miller_custom_extension"], grid->leafGridView(), log);
const Scaling reference{ 2.0, 2.0, 0.0 };
for (auto&& pos : pos_sample)
{
const auto value = scaling->evaluate(pos);
EXPECT_EQ(value, reference) << "Position : " << pos[0] << ", " << pos[1];
}
}
// evaluate the scaling
for (auto&& pos : positions) {
const auto value = scaling->evaluate(pos);
assert(value == reference);
}
TEST_F(ScalingFixture, MillerShift)
{
auto scaling =
create_scaling(scale_file["miller_shift"], grid->leafGridView(), log);
const Scaling ref_low{ 0.0, 0.0, 0.0 }, ref_high{ 1.0, 1.0, 0.0 };
for (auto&& pos : pos_sample)
{
const auto value = scaling->evaluate(pos);
auto reference = ref_low;
if (int(pos[1]) > 0)
reference = ref_high;
EXPECT_EQ(value, reference) << "Position : " << pos[0] << ", " << pos[1];
}
}
int main (int argc, char** argv)
TEST_F(ScalingFixture, MilPor)
{
try{
// initialize MPI if needed
auto& helper = Dune::MPIHelper::instance(argc, argv);
// build the logger
const auto log = Dune::Dorie::create_base_logger(helper,
spdlog::level::info);
constexpr int dim = 2;
// Read YAML file
std::string filename;
if (argc == 1) {
filename = "scaling.yml";
}
else if (argc == 2) {
filename = argv[1];
}
else {
throw std::runtime_error("Call this program with single argument <config>");
}
const YAML::Node param_file = YAML::LoadFile(filename);
log->info("Creating grid");
// get extensions from cfg
using Domain = typename Traits<dim>::Domain;
using retrieve_t = std::vector<double>;
const auto off = param_file["grid"]["offset"].template as<retrieve_t>();
const auto ext = param_file["grid"]["extensions"].template as<retrieve_t>();
Domain extensions;
std::copy(begin(ext), end(ext), extensions.begin());
Domain lower_left;
std::copy(begin(off), end(off), lower_left.begin());
const auto upper_right = lower_left + extensions;
// build scaling from grid extensions
using Grid = typename Traits<dim>::Grid;
std::array<int, dim> cells;
std::fill(begin(cells), end(cells), 10);
Grid grid(lower_left, upper_right, cells);
const auto grid_view = grid.leafGridView();
// test the scaling settings
for (auto&& scaling_cfg : param_file["scalings"]) {
log->info("Creating scaling: {}", scaling_cfg.first.as<std::string>());
test_scaling_setting<dim>(scaling_cfg.second, grid_view);
}