Commit f723c09f authored by Lukas Riedel's avatar Lukas Riedel

Implement Miller scaling in the new parameter interface

* Add ScalingType and ScalingFactor to parameter storage
* Apply scaling on top-level function return call
* Copy Miller scaling factors from old parameterization
* Adapt verification test
parent 688f03e8
......@@ -81,13 +81,17 @@ namespace Dune {
H5File& h5file; //!< HDF5 interface
const int verbose; //!< verbosity level
const Domain pf_extensions; //!< Extensions of the parameter field
public:
const bool millerSimilarity; //!< Switches miller similarity on or off
const RF variance; //!< Random field variance, needed for Miller Similarity
protected:
const Domain grid_extensions; //!< Spatial extensions of the model grid
const Domain scale; //!< Scaling factor of the random field
const Domain offset; //!< Position of the origin of the random field
public:
Parameter<T> millerField; //!< Random field used for Miller similarity
protected:
Vector gVector; //!< Unit length vector pointing towards gravitational force
Domain fieldCells; //!< Random field extensions in cells
......
......@@ -22,18 +22,24 @@ private:
using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<
LevelGridView
>;
// Define the storage and cache for parameterizations and scaling
using Index = typename Mapper::Index;
using ParameterStorage = std::map<
Index, std::shared_ptr<Dorie::RichardsParameterization<Traits>>
using Parameterization = Dorie::RichardsParameterization<Traits>;
using ScalingType = typename Parameterization::Scaling::Type;
using ScalingFactor = typename Parameterization::ScalingFactor;
using ParameterStorage = std::map<Index,
std::tuple<std::shared_ptr<Parameterization>,
ScalingType,
ScalingFactor
>
>;
/// Need this gruesome typedef because 'map::value_type' has const key
using Cache = std::pair<typename ParameterStorage::key_type,
typename ParameterStorage::mapped_type
>;
/// Typedef for accessing parameters
using RP = typename Dune::Dorie::RichardsParameterization<Traits>;
/// Configuration file tree
const Dune::ParameterTree& _config;
/// Grid view of the coarsest grid configuration (level 0)
......@@ -50,7 +56,7 @@ private:
/// Check if an entity has been cached
void verify_cache () const
{
if (not _cache.second){
if (not std::get<std::shared_ptr<Parameterization>>(_cache.second)) {
DUNE_THROW(Dune::Exception, "No parameterization cached."
<< "Call 'bind' before accessing the parameterization.");
}
......@@ -85,10 +91,17 @@ public:
std::function<RF(RF)> conductivity_f () const
{
verify_cache();
auto& par = _cache.second;
using Saturation = typename RP::Saturation;
return [&par](const RF saturation){
return par->conductivity_f()(Saturation{saturation}).value;
auto& par = std::get<std::shared_ptr<Parameterization>>(_cache.second);
RF scale = 1.0;
if (std::get<ScalingType>(_cache.second) == ScalingType::Miller) {
scale = std::get<ScalingFactor>(_cache.second).value;
}
using Saturation = typename Parameterization::Saturation;
return [&par, scale](const RF saturation){
return par->conductivity_f()(Saturation{saturation}).value
* scale * scale;
};
}
......@@ -96,10 +109,16 @@ public:
std::function<RF(RF)> saturation_f () const
{
verify_cache();
auto& par = _cache.second;
using MatricHead = typename RP::MatricHead;
return [&par](const RF matric_head){
return par->saturation_f()(MatricHead{matric_head}).value;
auto& par = std::get<std::shared_ptr<Parameterization>>(_cache.second);
RF scale = 1.0;
if (std::get<ScalingType>(_cache.second) == ScalingType::Miller) {
scale = std::get<ScalingFactor>(_cache.second).value;
}
using MatricHead = typename Parameterization::MatricHead;
return [&par, scale](const RF matric_head){
return par->saturation_f()(MatricHead{matric_head * scale}).value;
};
}
......@@ -121,6 +140,7 @@ public:
);
// read parameter values
using RP = Parameterization;
using MvG = MualemVanGenuchtenParameterization<Traits>;
typename RP::ResidualWaterContent theta_r{p_values.at(0)};
typename RP::SaturatedWaterContent theta_s{p_values.at(1)};
......@@ -140,9 +160,20 @@ public:
n
);
ScalingType type = RP::Scaling::None;
ScalingFactor scale{1.0};
if (p_base.millerSimilarity) {
type = ScalingType::Miller;
const auto val = get_param_value(&p_base.millerField);
scale = ScalingFactor{ std::exp(val - p_base.variance) };
}
// insert parameters into map
const auto index = _mapper.index(entity);
auto res = _param.emplace(index, mvg);
const auto res = _param.emplace(index,
std::make_tuple(mvg, type, scale)
);
if (not res.second) {
DUNE_THROW(Dune::Exception, "Parameterization for entity"
......
......@@ -40,6 +40,23 @@ public:
RF value;
};
/// Generic scaling factor
struct ScalingFactor
{
RF value;
};
/// Define scaling types for parameterization
struct Scaling
{
enum Type
{
Miller, //!< Geometric similarity
Warrick,
None //!< No scaling applied
};
};
/// Parameter defining the saturated hydraulic conductivity
struct SaturatedConductivity
{
......
......@@ -62,18 +62,21 @@ public:
auto sat_f = _param_new.saturation_f();
auto cond_f = _param_new.conductivity_f();
// do this 10 times to be sure
for (int i = 0; i < 10; ++i) {
const auto head = dist(rng);
// evaluate
const auto sat_old = this->param->saturation(head, pos);
const auto head_ref = this->param->headMillerToRef(head, pos);
const auto sat_old = this->param->saturation(head_ref, pos);
const auto sat_new = sat_f(head);
if (std::abs(sat_old - sat_new) > eps) {
DUNE_THROW(Dune::Exception, "Saturation does not match!");
}
const auto cond_old = this->param->K(pos)
const auto cond_ref = this->param->K(pos)
* this->param->condFactor(sat_old, pos);
const auto cond_old = this->param->condRefToMiller(cond_ref, pos);
const auto cond_new = cond_f(sat_new);
if (std::abs(cond_old - cond_new) > eps) {
DUNE_THROW(Dune::Exception, "Conductivity does not match!");
......
......@@ -18,6 +18,4 @@ grid.FEorder = 1
grid.extensions = 1 1
boundary.file = "{_asset_path}/bcs/infiltration_2d.dat"
parameters.arrayFile = {_asset_path}/parfields/fft_2d.h5
parameters.arrayFile = {_asset_path}/parfields/fft_2d.h5 , {_asset_path}/parfields/fft_2d_miller.h5 | expand
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment