scaling.hh 9.21 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11
#ifndef DUNE_DORIE_PARAM_SCALING_HH
#define DUNE_DORIE_PARAM_SCALING_HH

#include <string>
#include <cctype>
#include <vector>
#include <memory>
#include <algorithm>

#include <yaml-cpp/yaml.h>

12
#include <dune/dorie/common/logging.hh>
13 14
#include <dune/dorie/common/interpolator.hh>
#include <dune/dorie/common/h5file.hh>
15 16 17 18

namespace Dune {
namespace Dorie {

19 20 21 22
/// Define scaling factors used for any scaling type.
/** The variables are initialized to their default values
 *  (no effective scaling).
 *  \tparam RF Numeric type of scaling factors
23 24
 *  \see ScalingAdapter returns this object when evaluated at a certain position
 *  \ingroup RichardsParam
25 26 27 28 29
 */
template<typename RF>
struct ScalingFactors
{
    //! Scaling factor of matric head
30
    RF scale_head = 1.0;
31
    //! Scaling factor of conductivity
32
    RF scale_cond = 1.0;
33
    //! Scaling factor of porosity (a summand, actually!)
34
    RF scale_por = 0.0;
35 36 37
};

/// Abstract base class for scaling adapters. Provides the interface.
38 39 40 41 42 43 44 45
/** A ScalingAdapter implements the evaluate method which returns a set of
 *  ScalingFactors.
 *
 *  \tparam Traits Type traits of this adapter.
 *  \ingroup RichardsParam
 *  \author Lukas Riedel
 *  \date 2018
 *  \see ScalingAdapterFactory conveniently creates derived objects
46 47 48 49 50 51 52 53 54 55 56
 */
template<typename Traits>
class ScalingAdapter
{
protected:
    using RF = typename Traits::RF;
    using Domain = typename Traits::Domain;
    using return_t = ScalingFactors<RF>;

    //! Interpolator storage
    std::vector<std::shared_ptr<Interpolator<RF, Traits>>> _interpolators;
57 58
    /// Logger of this instance
    const std::shared_ptr<spdlog::logger> _log;
59 60 61 62 63

public:
    /// Create the adapter and store config data
    /** Take a YAML node as unnamed parameter to enforce derived classes to
     *  take one as argument, too.
64
     *  \log The logger for this instance (defaults to the Richards logger)
65
     */
66 67 68 69 70
    explicit ScalingAdapter (
        const YAML::Node&,
        const std::shared_ptr<spdlog::logger> log
    ):
        _log(log)
71 72 73 74 75 76 77 78 79 80 81 82
    { }

    /// Evaluate this adapter.
    virtual return_t evaluate (const Domain& pos) const = 0;

    /// Destroy this adapter.
    virtual ~ScalingAdapter () = default;

protected:
    /// Add an interpolator based on the data in the config node
    /** \param cfg "data" node to read the interpolator data from
     *      This is typically a sub-node of the node given in the constructor.
83
     *  \param node_name Name of the ``cfg`` node, for nice error messages.
84
     */
85
    void add_interpolator (const YAML::Node& cfg, const std::string node_name)
86
    {
87 88 89
        _log->trace("Adding interpolator from data node: {}", node_name);
        const auto dim = Traits::dim;

90 91 92
        // open the H5 file
        const auto filename = cfg["file"].as<std::string>();
        const auto dataset = cfg["dataset"].as<std::string>();
93
        H5File h5file(filename, _log);
94 95 96 97 98

        // read the dataset
        std::vector<RF> data;
        std::vector<size_t> shape;
        h5file.read_dataset(dataset, H5T_NATIVE_DOUBLE, data, shape);
99
        std::reverse(begin(shape), end(shape));
100 101 102 103

        // read extensions and offset
        const auto ext = cfg["extensions"].as<std::vector<RF>>();
        const auto off = cfg["offset"].as<std::vector<RF>>();
104 105 106 107 108 109 110 111 112 113 114
        if (shape.size() != dim) {
            _log->error("Expected {}-dimensional dataset. File: {}, "
                        "Dataset: {}, Dimensions: {}",
                        dim, filename, dataset, shape.size());
            DUNE_THROW(Dune::IOError, "Invalid dataset for interpolator");
        }
        if (ext.size() != dim) {
            _log->error("Expected {}-dimensional sequence as 'extensions' "
                        "at scaling section node: {}",
                        dim, node_name);
            DUNE_THROW(Dune::IOError, "Invalid interpolator input");
115
        }
116 117 118 119 120
        if (off.size() != dim) {
            _log->error("Expected {}-dimensional sequence as 'offset' "
                        "at scaling section node: {}",
                        dim, node_name);
            DUNE_THROW(Dune::IOError, "Invalid interpolator input");
121 122 123 124 125 126 127 128 129 130
        }

        // copy into correct data structure
        Domain extensions;
        std::copy(begin(ext), end(ext), extensions.begin());
        Domain offset;
        std::copy(begin(off), end(off), offset.begin());

        // build interpolator
        this->_interpolators.emplace_back(
131
            InterpolatorFactory<Traits>::create(
132 133 134 135
                cfg["interpolation"].template as<std::string>(),
                data,
                shape,
                extensions,
136 137
                offset,
                _log
138 139 140 141 142
            )
        );
    }
};

143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168
/// A ScalingAdapter implementing Miller scaling.
/** Miller scaling assumes geometric similarity and only varies the local
 *  length scale of the porous medium.
 *
 *  The scaled values of matric head \f$ h_m \f$ and the conducitivity
 *  \f$ K \f$ with respect to the reference values \f$ h_m^* \f$ and
 *  \f$ K^* \f$ are calculated by applying the Miller scaling factor
 *  \f$ \xi \f$, which indicates a local length scale with respect to the
 *  reference length scale of the medium.
 *
 *  \f{align*}{
 *      h_m &= h_m^* \cdot \xi \\
 *      K &= K^* \cdot \xi^2
 *  \f}
 *
 *  When entering data, the single dataset is used as Miller scaling factor
 *  \f$ \xi \f$. This scaling adapter effectively sets
 *
 *  \f{align*}{
 *      \xi_{h_m} = \xi_K = \xi
 *  \f}
 *
 *  \ingroup RichardsParam
 *  \author Lukas Riedel
 *  \date 2018
 */
169 170 171 172 173 174 175 176 177
template<typename Traits>
class MillerScalingAdapter : public ScalingAdapter<Traits>
{
private:
    using Base = ScalingAdapter<Traits>;
    using Domain = typename Base::Domain;
    using return_t = typename Base::return_t;

public:
178
    /// Create this adapter and the internal interpolator.
179 180 181 182 183
    explicit MillerScalingAdapter (
        const YAML::Node& scale_data_cfg,
        const std::shared_ptr<spdlog::logger> log
    ):
        Base(scale_data_cfg, log)
184
    {
185
        this->add_interpolator(scale_data_cfg["scale_miller"], "scale_miller");
186 187 188 189 190 191 192 193
    }

    /// Export type of this scaling adapter
    static inline std::string type = "Miller";

    /// Evaluate this interpolator
    return_t evaluate (const Domain& pos) const override
    {
194
        const auto scaling_value = this->_interpolators[0]->evaluate(pos);
195
        return return_t{scaling_value, scaling_value};
196 197 198 199 200 201
    }

    /// Delete this adapter
    ~MillerScalingAdapter () override = default;
};

202 203 204 205 206 207 208 209 210
/// A ScalingAdapter that does not implement any scaling.
/** This adapter returns a default-initialized instance of ScalingFactors
 *  when evaluated.
 *
 *  \see ScalingFactors
 *  \ingroup RichardsParam
 *  \author Lukas Riedel
 *  \date 2018
 */
211 212 213 214 215 216 217 218 219 220
template<typename Traits>
class DummyScalingAdapter : public ScalingAdapter<Traits>
{
private:
    using Base = ScalingAdapter<Traits>;
    using Domain = typename Base::Domain;
    using return_t = typename Base::return_t;

public:
    /// Create this adapter. Don't initialize members.
221 222 223 224 225
    explicit DummyScalingAdapter (
        const YAML::Node& scale_data_cfg,
        const std::shared_ptr<spdlog::logger> log
    ):
        Base(scale_data_cfg, log)
226 227 228 229 230 231 232 233
    { }

    /// Export type of this scaling adapter
    static inline std::string type = "none";

    /// Evaluate this scaling (always returns no scaling)
    return_t evaluate (const Domain& pos) const override
    {
234
        return return_t{};
235 236 237 238 239 240 241 242
    }

    /// Delete this adapter
    ~DummyScalingAdapter () override = default;
};


/// Factory class for creating ScalingAdapters.
243 244 245 246 247 248
/**
 *  \ingroup RichardsParam
 *  \author Lukas Riedel
 *  \date 2018
 *  \see ScalingAdapter is the base class of instances created by this factory
 */
249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272
template<typename Traits>
struct ScalingAdapterFactory
{
private:
    using RF = typename Traits::RF;

    /// Turn a string into lowercase characters
    /** \param in Input string
     *  \return Lowercase output string
     */
    static std::string to_lower (const std::string& in)
    {
        std::string out = in;
        std::transform(in.begin(), in.end(), out.begin(),
            [](unsigned char c){ return std::tolower(c); }
        );
        return out;
    }

public:

    /// Create a scaling adapter using a YAML config node
    /** \param type Type of the scaling adapter to be created
     *  \param config YAML config node specifying the scale adapter settings
273 274
     *  \log The logger to use in the created object. Defaults to
     *      Dune::Dorie::log_richards.
275 276 277 278
     *  \return Shared pointer to the scaling adapter
     */
    static std::shared_ptr<ScalingAdapter<Traits>> create (
        const std::string& type,
279 280
        const YAML::Node& config,
        const std::shared_ptr<spdlog::logger> log=get_logger(log_richards)
281 282
    )
    {
283 284
        log->debug("Creating scaling adapter of type: {}", type);

285
        if (type == MillerScalingAdapter<Traits>::type) {
286
            return std::make_shared<MillerScalingAdapter<Traits>>(config, log);
287 288
        }
        else if (to_lower(type) == DummyScalingAdapter<Traits>::type) {
289
            return std::make_shared<DummyScalingAdapter<Traits>>(config, log);
290 291
        }
        else {
292 293
            log->error("Unknown scaling type: {}", type);
            DUNE_THROW(Dune::NotImplemented, "Unknown scaling type");
294 295 296 297 298 299 300 301 302
        }
    }
};


} // namespace Dorie
} // namespace Dune

#endif // DUNE_DORIE_PARAM_SCALING_HH