interface.hh 9.85 KB
Newer Older
1 2 3
#ifndef DUNE_DORIE_PARAM_INTERFACE_HH
#define DUNE_DORIE_PARAM_INTERFACE_HH

4 5 6 7 8
#include <cmath>
#include <vector>
#include <map>
#include <utility>
#include <functional>
9 10 11 12 13
#include <memory>

#include <dune/common/exceptions.hh>
#include <dune/grid/common/mcmgmapper.hh>

14 15
#include <dune/dorie/model/richards/parametrization/richards.hh>
#include <dune/dorie/model/richards/parametrization/mualem_van_genuchten.hh>
16 17 18 19

namespace Dune {
namespace Dorie {

Lukas Riedel's avatar
Lukas Riedel committed
20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
/// Top-level parameterization interface.
/** This class stores the parameters, maps them to grid entities, and provides 
 *  functions to access the parameterization. For doing so, the
 *  FlowParameters::bind function has to be called first to cache the
 *  parameters of the given grid entity. Parameters are only stored for entities
 *  of codim 0.
 * 
 *  \detail The parameterization consists of three structures:
 *      - This top-level class serves as interface and storage. It also casts to
 *          data types used by the DUNE operators
 *      - The Dorie::RichardsParameterization interface that provides strong 
 *          data types for the base parameters and purely virtual functions that
 *          need to be defined by derived parameterizations
 *      - The actual parameterization defining all parameters and functions
 */
35 36 37 38 39 40 41 42 43 44
template <typename Traits>
class FlowParameters
{
private:
    using RF = typename Traits::RF;
    using Grid = typename Traits::Grid;
    using LevelGridView = typename Grid::LevelGridView;
    using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<
        LevelGridView
    >;
45
    using Vector = typename Traits::Vector;
46 47

    // Define the storage and cache for parameterizations and scaling
Lukas Riedel's avatar
Lukas Riedel committed
48
    /// Index type used by the element mapper
49
    using Index = typename Mapper::Index;
Lukas Riedel's avatar
Lukas Riedel committed
50
    /// Base class of the parameterization
51
    using Parameterization = Dorie::RichardsParameterization<Traits>;
Lukas Riedel's avatar
Lukas Riedel committed
52

53 54
    using ScalingType = typename Parameterization::Scaling::Type;
    using ScalingFactor = typename Parameterization::ScalingFactor;
Lukas Riedel's avatar
Lukas Riedel committed
55 56

    /// Storage for parameters and skaling factors
57 58 59 60 61
    using ParameterStorage = std::map<Index,
        std::tuple<std::shared_ptr<Parameterization>,
                   ScalingType,
                   ScalingFactor
        >
62
    >;
63

64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
    /// Need this gruesome typedef because 'map::value_type' has const key
    using Cache = std::pair<typename ParameterStorage::key_type,
        typename ParameterStorage::mapped_type
    >;

    /// Configuration file tree
    const Dune::ParameterTree& _config;
    /// Grid view of the coarsest grid configuration (level 0)
    const LevelGridView _gv;
    /// Mapper for mapping grid elements to indices
    Mapper _mapper;
    /// Map for storing parameterization information for each grid entity
    ParameterStorage _param;
    /// Currently cached parameterization (bound to element)
    mutable Cache _cache;
79 80
    /// Gravity vector
    Vector _gravity;
81 82 83 84 85 86

private:

    /// Check if an entity has been cached
    void verify_cache () const
    {
87
        if (not std::get<std::shared_ptr<Parameterization>>(_cache.second)) {
88
            DUNE_THROW(Dune::Exception, "No parameterization cached. "
89 90 91 92 93
                << "Call 'bind' before accessing the parameterization.");
        }
    }

public:
Lukas Riedel's avatar
Lukas Riedel committed
94 95 96 97 98 99
    /// Constructor.
    /** Create a level grid view of level 0, create a mapper on it, and store
     *  the config tree.
     *  \param config Configuration file tree
     *  \param grid Shared pointer to the grid
     */
100 101 102 103 104 105
    FlowParameters (
        const Dune::ParameterTree& config,
        std::shared_ptr<Grid> grid
    ):
        _config(config),
        _gv(grid->levelGridView(0)),
106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
        _mapper(_gv, Dune::mcmgElementLayout()),
        _gravity(0.0)
    {
        _gravity[Traits::dim - 1] = -1.0;
    }

    /// Return normalized gravity vector
    const Vector& gravity() const { return _gravity; }

    /// Return the current cache
    const Cache& cache() const
    {
        verify_cache();
        return _cache;
    }
121

Lukas Riedel's avatar
Lukas Riedel committed
122 123
    /// Bind to a grid entity. Required to call parameterization functions!
    /** \param entity Grid entity (codim 0) to bind to
124 125
     */
    template <class Entity>
126
    void bind (Entity entity) const
127
    {
128 129 130 131
        // retrieve index of top father element of chosen entity
        while(entity.hasFather()) {
            entity = entity.father();
        }
132
        const auto index = _mapper.index(entity);
Lukas Riedel's avatar
Lukas Riedel committed
133 134 135 136 137 138 139 140

        // done if index is already cached
        if (index == _cache.first
            && std::get<std::shared_ptr<Parameterization>>(_cache.second))
        {
            return;
        }

141
        // otherwise find the index in the stored map and cache the values
142 143
        const auto it = _param.find(index);
        if (it == _param.end()) {
144
            DUNE_THROW(Dune::Exception, "Could not retrieve parameterization "
145 146 147 148 149
                << "for entity " << index);
        }
        _cache = *it;
    }

Lukas Riedel's avatar
Lukas Riedel committed
150 151 152 153 154
    /// Return a scaled version of the Conductivity function
    /** Uses the function of the underlying parameterization and applies
     *  scaling as specified by SkalingType. Cast to operator-internal RF.
     *  \return Function: Saturation -> Conductivity
     */
155 156 157
    std::function<RF(RF)> conductivity_f () const
    {
        verify_cache();
Lukas Riedel's avatar
Lukas Riedel committed
158 159 160 161
        const auto& par =
            std::get<std::shared_ptr<Parameterization>>(_cache.second);
        const auto& type = std::get<ScalingType>(_cache.second);
        const auto cond_f = par->conductivity_f();
162
        using Saturation = typename Parameterization::Saturation;
163

164 165 166 167 168
        if (type == ScalingType::Miller) {
            auto& scale = std::get<ScalingFactor>(_cache.second).value;
            return [cond_f, &scale](const RF saturation){
                return cond_f(Saturation{saturation}).value * scale * scale;
            };
169 170
        }

171 172
        return [cond_f](const RF saturation){
            return cond_f(Saturation{saturation}).value;
173 174 175
        };
    }

Lukas Riedel's avatar
Lukas Riedel committed
176 177 178 179 180
    /// Return a scaled version of the Saturation function
    /** Uses the function of the underlying parameterization and applies
     *  scaling as specified by SkalingType. Cast to operator-internal RF.
     *  \return Function: Matric Head -> Saturation
     */
181 182 183
    std::function<RF(RF)> saturation_f () const
    {
        verify_cache();
Lukas Riedel's avatar
Lukas Riedel committed
184 185 186 187
        const auto& par = 
            std::get<std::shared_ptr<Parameterization>>(_cache.second);
        const auto& type = std::get<ScalingType>(_cache.second);
        const auto sat_f = par->saturation_f();
188
        using MatricHead = typename Parameterization::MatricHead;
189

190
        if (type == ScalingType::Miller) {
Lukas Riedel's avatar
Lukas Riedel committed
191
            const auto& scale = std::get<ScalingFactor>(_cache.second).value;
192 193 194
            return [sat_f, &scale](const RF matric_head){
                return sat_f(MatricHead{matric_head * scale}).value;
            };
195 196
        }

197 198
        return [sat_f](const RF matric_head){
            return sat_f(MatricHead{matric_head}).value;
199 200 201
        };
    }

Lukas Riedel's avatar
Lukas Riedel committed
202 203 204 205 206
    /// Return a scaled version of the Water Content function
    /** Uses the function of the underlying parameterization and applies
     *  scaling as specified by SkalingType. Cast to operator-internal RF.
     *  \return Function: Saturation -> Water Content
     */
207 208 209
    std::function<RF(RF)> water_content_f () const
    {
        verify_cache();
Lukas Riedel's avatar
Lukas Riedel committed
210 211 212
        const auto& par =
            std::get<std::shared_ptr<Parameterization>>(_cache.second);
        const auto wc_f = par->water_content_f();
213 214 215 216 217 218 219
        using Saturation = typename Parameterization::Saturation;

        return [wc_f](const RF saturation) {
            return wc_f(Saturation{saturation}).value;
        };
    }

Lukas Riedel's avatar
Lukas Riedel committed
220 221 222 223
    /// Copy parameters from the 'old' parameterization implementation
    /** \todo Remove once new interface is established
     *  \param p_base Old parameterization base class
     */
224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241
    template <class Param>
    void copy_from_old_parameters (const Param& p_base)
    {
        // iterate over grid elements
        for (auto&& entity: elements(_gv))
        {
            const auto pos = entity.geometry().center();
            auto get_param_value = [&p_base, &pos](auto param) {
                return p_base.interpolate(*param, pos);
            };

            std::vector<RF> p_values;
            std::transform(p_base.parameters.begin(), p_base.parameters.end(),
                std::back_inserter(p_values),
                get_param_value
            );

            // read parameter values
242
            using RP = Parameterization;
243
            using MvG = MualemVanGenuchtenParameterization<Traits>;
Lukas Riedel's avatar
Lukas Riedel committed
244 245
            typename RP::ResidualWaterContent theta_r{p_values.at(0)};
            typename RP::SaturatedWaterContent theta_s{p_values.at(1)};
246 247 248
            typename MvG::Alpha alpha{p_values.at(2)};
            typename MvG::N n{p_values.at(3)};
            typename MvG::Tau tau{p_values.at(4)};
Lukas Riedel's avatar
Lukas Riedel committed
249
            typename RP::SaturatedConductivity k0{p_values.at(5)};
250

251
            // can create parameters out of order by using tuple
252 253 254
            auto mvg =
                std::make_shared<MualemVanGenuchtenParameterization<Traits>>
            (
255
                std::make_tuple(theta_s, theta_r, k0, alpha, tau, n)
256 257
            );

258 259 260 261 262 263 264 265 266
            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) };
            }

267 268
            // insert parameters into map
            const auto index = _mapper.index(entity);
269 270 271
            const auto res = _param.emplace(index,
                                            std::make_tuple(mvg, type, scale)
            );
272 273 274 275 276 277 278 279 280 281 282 283 284 285

            if (not res.second) {
                DUNE_THROW(Dune::Exception, "Parameterization for entity"
                    << index << " already exists."
                );
            }
        }
    }
};

} // namespace Dune
} // namespace Dorie

#endif // DUNE_DORIE_PARAM_INTERFACE_HH