diff --git a/CHANGELOG.md b/CHANGELOG.md index 4dfbf78e3ff345779c5e3c1f3526272c2c41d208..9771c70fd62a994ccff4c6a1e8b68af8c1b5b52b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -40,6 +40,7 @@ * Coupling between transient models for water flow and solute transport !96 * Initial conditions generated from H5 input data !130 * Generic Python VTK file reader !143 +* Linear interpolator for initial conditions and scaling fields !145 * Parameterizations for hydrodynamic dispersion in solute transport !141 * Generic Python VTK file reader !143, !150 diff --git a/doc/default_files/richards-parameters.xml b/doc/default_files/richards-parameters.xml index 5edd54f1e28efe4edde9a02b5d9f182a1b055b39..86948249a4e561f1b969461546472fcc183af5b5 100644 --- a/doc/default_files/richards-parameters.xml +++ b/doc/default_files/richards-parameters.xml @@ -153,8 +153,8 @@ adding an empty line, make text **bold** or ``monospaced``. Interpolation type used for the data (``data`` type only). - nearest - nearest + nearest, linear + linear diff --git a/doc/default_files/transport-parameters.xml b/doc/default_files/transport-parameters.xml index 35d95f2fb775339c1a7f4b3e1d532c9c273972f7..352381ca549d10fc1c56455c954827ba369d13f2 100644 --- a/doc/default_files/transport-parameters.xml +++ b/doc/default_files/transport-parameters.xml @@ -170,7 +170,7 @@ adding an empty line, make text **bold** or ``monospaced``. Interpolation type used for the data (``data`` type only). - nearest + nearest, linear nearest diff --git a/doc/index.rst b/doc/index.rst index b82b96ff3ea36b2f311d03723a34b93d284181eb..524e0601fc981d70cc2b9ddfacbc474c583988e8 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -33,6 +33,7 @@ assignment and increment are based on the :doc:`public-api`. manual/bcfile manual/initial manual/fluxes + manual/interpolators public-api .. toctree:: diff --git a/doc/manual/initial.rst b/doc/manual/initial.rst index f7f68fe3eb9205c6c27448ba5c57e56882e4119f..f535f76aa9dd2ded375c866ecda0345dad51a0ec 100644 --- a/doc/manual/initial.rst +++ b/doc/manual/initial.rst @@ -14,7 +14,7 @@ long as the respective quantity has a unique transformation to the solver solution quantity. Initial condition input is controlled entirely via the -:doc:`Configuration File `. +:doc:`Configuration File `. .. note:: The initial condition is projected onto the actual solution function space. @@ -74,6 +74,8 @@ They are controlled by the ``initial.type`` key and available for every model. can be chosen using the setting ``initial.interpolation``. The input data is automatically streched to match the grid extensions. + .. note:: For ``FEorder > 0``, linear interpolation is recommended. + Supported file extensions: * ``.h5``: H5_ data file. ``initial.dataset`` may be a file-internal path diff --git a/doc/manual/interpolators.rst b/doc/manual/interpolators.rst new file mode 100644 index 0000000000000000000000000000000000000000..d49fa6f245439f12dd8981f5920a2ce107186ef0 --- /dev/null +++ b/doc/manual/interpolators.rst @@ -0,0 +1,59 @@ +.. _man-interpolators: + +Interpolators +============= + +Interpolators are used at several points in DORiE to evaluate discrete input +data at every physical location on the grid. + +Currently, interpolators are used in the +:doc:`Parameter File `, and when loading an initial condition +from data with the options given in the +:doc:`Configuration File `. + +.. note:: The type of interpolator may change how the input data is + interpreted. In particular, different interpolators require + different input dataset shapes for the same grid configuration. + +Interpolator Types +------------------ + +Every interpolator assumes that the input data spans a rectangle (2D) or cuboid +(3D). In case of initial conditions, the respective volume is streched to cover +the entire grid. For scaling fields, users may specify extensions and offset of +the volume, or opt for it to cover the entire grid by omitting the respective +keys in the parameter file. + +.. object:: Nearest neighbor interpolator + + * ``interpolation: nearest`` + + Interprets dataset values as *cell values*. No extrapolation is applied. + The lowest supported dataset dimensions are ``(1, 1)`` (single value + everywhere). + + Use this interpolator for cell-centered data like + :ref:`scaling field values ` and initial conditions + for a finite volume solver. + + .. tip:: To provide the exact scaling factors for each cell of a grid + of :math:`1 \times 10` cells, use the ``nearest`` interpolator + with a dataset of shape ``(10, 1)`` (dimensions are inverted). + The dataset value are the function values at the cell barycenters. + +.. object:: Linear interpolator + + * ``interpolation: linear`` + + Interprets dataset values as *vertex values* and linearaly interpolates + between them. No extrapolation is applied. The lowest supported dataset + dimensions are ``(2, 2)`` (one value in each corner). + + Use this interpolator if input data should be interpreted as + continuous functions like initial conditions for a DG solver. + + .. tip:: To provide the initial condition for a DG solver with finite + element polynomials of order :math:`k = 1` on a grid of + :math:`1 \times 10` cells, use the ``linear`` interpolator + with a dataset of shape ``(11, 2)`` (dimensions are inverted). + The dataset values are the function values at the grid vertices. diff --git a/doc/manual/parameter-file.rst b/doc/manual/parameter-file.rst index 7153dd7237f2d62da93fbcf2ab9128723e80b24a..0df001f24b1e07c724fe213d98ae4abaee56a135 100644 --- a/doc/manual/parameter-file.rst +++ b/doc/manual/parameter-file.rst @@ -54,12 +54,11 @@ of scaling. Every ``scaling_section`` has the same keys: The H5 filepath is given by ``file``, and the internal dataset path by ``dataset``. You can then choose an -``interpolation`` method for the dataset, which requires you to insert values +``interpolation`` method for the dataset. You may optionally insert values for the physical ``extensions`` of the dataset and the ``offset`` of its -(front) lower left corner from the origin. The latter two keys are optional. -If at least one of them is omitted, the inserted field is automatically -scaled to match the maximum grid extensions. This also works for irregular grid -shapes. +(front) lower left corner from the origin. If at least one of them is omitted, +the inserted field is automatically scaled to match the maximum grid +extensions. This also works for irregular grid shapes. .. code-block:: yaml @@ -108,9 +107,10 @@ One or more datasets in possibly multiple HDF5_ files are required for creating a medium with small-scale heterogeneities. The datasets must consist of floating-point values and must have the same dimensionality as the grid. Other than that, there is no restriction on their shape because they are -inserted into interpolators. The interpolators supply scaling factor values -based on positions on the grid and thus create a grid-independent -representation. Scaling field input works the same for any supported grid type. +inserted into :doc:`Interpolators `. The interpolators supply +scaling factor values ased on positions on the grid and thus create a +grid-independent representation. Scaling field input works the same for any +supported grid type. During setup, the program reads the interpolator values at the **barycenter** of every grid cell on the **coarsest** grid configuration. These values are @@ -120,8 +120,7 @@ scaling factors apply in all its child cells. Supported Parameter File Types ------------------------------ This section lists the available types for parameterizations, scalings, and -interpolators in a understandable format. You can also dig through the code -documentation below, or the code itself, to retrieve this information. +interpolators in a understandable format. Richards Parameterizations ~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -544,17 +543,6 @@ Miller Scaling scale_miller: # ... -.. _man-interpolators: - -Interpolators -------------- - -.. object:: Nearest neighbor interpolator - - Interprets dataset values as cell values. - - * ``interpolation: nearest`` - .. _YAML: https://gettaurus.org/docs/YAMLTutorial/ .. _HDF5: https://www.h5py.org/ diff --git a/dune/dorie/common/initial_condition/h5.hh b/dune/dorie/common/initial_condition/h5.hh index 10a15942e8370388ddeb75d63575e40ac26c3dbe..7d87eb51c0d7c4ac55ab1dc709be4498c1ba6431 100644 --- a/dune/dorie/common/initial_condition/h5.hh +++ b/dune/dorie/common/initial_condition/h5.hh @@ -61,7 +61,7 @@ public: } private: - std::shared_ptr> _interpolator; + std::shared_ptr> _interpolator; }; } // namespace Dorie diff --git a/dune/dorie/common/interpolator.hh b/dune/dorie/common/interpolator.hh index d1c07d31d7ba8b43ad06dd6c2aa71f07bf83f0e9..6ff3aa1037091374e4b1327bd0e0a328ebf6a4d7 100644 --- a/dune/dorie/common/interpolator.hh +++ b/dune/dorie/common/interpolator.hh @@ -7,6 +7,8 @@ #include #include #include +#include +#include #include #include @@ -28,14 +30,16 @@ namespace Dorie { * \author Lukas Riedel * \date 2018 */ -template +template class Interpolator { protected: - using Domain = typename Traits::Domain; - static constexpr int dim = Traits::dim; + /// Coordinate Type + using DF = double; + /// Physical vector type + using Domain = Dune::FieldVector; - const std::vector _data; //!< contiguously stored data of the array + const std::vector _data; //!< contiguously stored data of the array const std::vector _shape; //!< inverted shape of the original array //! physical extensions of the dataset const Domain _extensions; @@ -100,7 +104,7 @@ public: /** \param pos Position to evaluate * \return Value of the interpolated function */ - virtual T evaluate (const Domain& pos) const = 0; + virtual DataType evaluate (const Domain& pos) const = 0; }; /// A nearest-neighbor interpolator in 2D and 3D @@ -110,34 +114,25 @@ public: * * \ingroup Interpolators */ -template -class NearestNeighborInterpolator: public Interpolator +template +class NearestNeighborInterpolator: public Interpolator { private: - using Base = Interpolator; + using Base = Interpolator; using Domain = typename Base::Domain; - static constexpr int dim = Base::dim; - using DF = typename Traits::DF; + using DF = typename Base::DF; static constexpr DF eps = 1e-9; public: - template - NearestNeighborInterpolator (Data&& data, - Shape&& shape, - Domain1&& extensions, - Domain2&& offset) : - Base(std::forward(data), - std::forward(shape), - std::forward(extensions), - std::forward(offset)) - { } + /// Re-use the base class constructor + using Base::Base; /// Export the type of this intepolator static inline std::string type = "nearest"; /// Evaluate the interpolator - T evaluate (const Domain& x) const override + DataType evaluate (const Domain& x) const override { const auto pos = this->sub_offset(x); if (not this->inside_extensions(pos)) { @@ -168,6 +163,158 @@ public: ~NearestNeighborInterpolator() override = default; }; +/// Linear interpolator for 2D and 3D. +/** The implementation is based on articles + * https://en.wikipedia.org/wiki/Bilinear_interpolation + * and https://en.wikipedia.org/wiki/Trilinear_interpolation + * + * \ingroup Interpolators + * \author Lukas Riedel + * \date 2019 + */ +template +class LinearInterpolator : public Interpolator +{ +private: + using Base = Interpolator; + using Domain = typename Base::Domain; + using DF = typename Base::DF; + using MultiIdx = std::array; + + /// Get the position vector in mesh units + /** \param pos Global position vector + * \return Position vector in mesh units + */ + Domain get_pos_unit (const Domain& pos) const + { + Domain pos_unit; + for (size_t i = 0; i < pos.size(); ++i) { + pos_unit[i] = pos[i] * (this->_shape[i] - 1) + / this->_extensions[i]; + } + + return pos_unit; + } + + /// Return the position difference vector in mesh units + /** The difference is calculated by using the (front) lower left corner + * as origin. + * \param pos_unit Position vector in mesh units, see get_pos_unit() + * \pararm indices The multi index of (front) lower left corner + * \return Position difference vector in mesh units + */ + Domain get_pos_unit_diff (const Domain& pos_unit, const MultiIdx& indices) + const + { + Domain pos_diff; + for (size_t i = 0; i < pos_unit.size(); ++i) { + pos_diff[i] = pos_unit[i] - indices[i]; + } + return pos_diff; + } + + /// Get the indices for accessing the data knots + /** \param pos_unit The position in mesh units + * \return Indices in every dimension of the mesh + */ + std::pair get_knot_indices ( + const Domain& pos_unit + ) const + { + // round mesh unit positions up and down + MultiIdx idx_lower, idx_upper; + for (size_t i = 0; i < idx_lower.size(); ++i) { + idx_lower[i] = std::max(std::floor(pos_unit[i]), 0.0); + idx_upper[i] = std::min(std::ceil(pos_unit[i]), + this->_shape[i]-1); + } + + return std::make_pair(idx_lower, idx_upper); + } + + /// Transform a multi index into a index for accessing the data array + /** \param indices The multi index to transform + * \return The index to query the data set #_data + */ + size_t to_index (const MultiIdx& indices) const + { + size_t index = indices[0] + indices[1] * this->_shape[0]; + if constexpr (dim == 3) { + index += indices[2] * this->_shape[1] * this->_shape[0]; + } + + return index; + } + +public: + /// Export the type of this intepolator + static inline constexpr std::string_view type = "linear"; + + /// Re-use the base class constructor + using Base::Base; + + /// Delete this interpolator + ~LinearInterpolator () override = default; + + /// Evaluate the data at a global position + DataType evaluate (const Domain& x) const override + { + const auto pos = this->sub_offset(x); + if (not this->inside_extensions(pos)) { + DUNE_THROW(Exception, "Querying interpolator data outside its " + "extensions!"); + } + + // compute indices and position vectors in mesh units + const auto pos_unit = get_pos_unit(pos); + const auto [idx_lower, idx_upper] = get_knot_indices(pos); + const auto pos_unit_diff = get_pos_unit_diff(pos_unit, idx_lower); + const auto& data = this->_data; + + + DataType c_00, c_01, c_10, c_11; + DF y_d, z_d; + + // actual computation + if constexpr (dim == 2) { + c_00 = data[to_index(idx_lower)]; + c_01 = data[to_index({idx_lower[0], idx_upper[1]})]; + c_10 = data[to_index({idx_upper[0], idx_lower[1]})]; + c_11 = data[to_index(idx_upper)]; + + // use y, z here for common computation after if-clause + y_d = pos_unit_diff[0]; // actually x + z_d = pos_unit_diff[1]; // actually y + } + else { + DataType c_000, c_001, c_010, c_100, c_011, c_110, c_101, c_111; + c_000 = data[to_index(idx_lower)]; + c_100 = data[to_index({idx_upper[0], idx_lower[1], idx_lower[2]})]; + c_010 = data[to_index({idx_lower[0], idx_upper[1], idx_lower[2]})]; + c_001 = data[to_index({idx_lower[0], idx_lower[1], idx_upper[2]})]; + c_011 = data[to_index({idx_lower[0], idx_upper[1], idx_upper[2]})]; + c_110 = data[to_index({idx_upper[0], idx_upper[1], idx_lower[2]})]; + c_101 = data[to_index({idx_upper[0], idx_lower[1], idx_upper[2]})]; + c_111 = data[to_index(idx_upper)]; + + DF x_d = pos_unit_diff[0]; + y_d = pos_unit_diff[1]; + z_d = pos_unit_diff[2]; + + c_00 = c_000 * (1 - x_d) + c_100 * x_d; + c_01 = c_001 * (1 - x_d) + c_101 * x_d; + c_10 = c_010 * (1 - x_d) + c_110 * x_d; + c_11 = c_011 * (1 - x_d) + c_111 * x_d; + } + + const DataType c_0 = c_00 * (1 - y_d) + c_10 * y_d; + const DataType c_1 = c_01 * (1 - y_d) + c_11 * y_d; + + const DataType c = c_0 * (1 - z_d) + c_1 * z_d; + return c; + } +}; + /// Factory for creating interpolators /** * \see Interpolator for the base class of created objects @@ -193,29 +340,38 @@ using Domain = typename Traits::Domain; * Dune::Dorie::log_base * \return Shared pointer to the interpolator */ -template +template static auto create ( const std::string& type, Data&& data, Shape&& shape, - Domain&& extensions, - Domain&& offset, + SpaceVector&& extensions, + SpaceVector&& offset, const std::shared_ptr log=get_logger(log_base) ) -> std::shared_ptr::value_type, - Traits>> + typename std::decay_t::value_type, + Traits::dim>> { log->debug("Creating interpolator of type: {}", type); - using data_t = typename std::remove_reference_t::value_type; - using NNInterp = NearestNeighborInterpolator; + using data_t = typename std::decay_t::value_type; + using NNInterp = NearestNeighborInterpolator; + using LinearInterp = LinearInterpolator; if (type == NNInterp::type) { return std::make_shared( std::forward(data), std::forward(shape), - std::forward(extensions), - std::forward(offset) + std::forward(extensions), + std::forward(offset) + ); + } + else if (type == LinearInterp::type) { + return std::make_shared( + std::forward(data), + std::forward(shape), + std::forward(extensions), + std::forward(offset) ); } else { @@ -244,8 +400,8 @@ static auto create ( const std::shared_ptr log=get_logger(log_base) ) -> std::shared_ptr::value_type, - Traits>> + typename std::decay_t::value_type, + Traits::dim>> { /// retrieve extensions and offset from grid extensions const auto level_gv = grid_view.grid().levelGridView(0); @@ -272,7 +428,7 @@ static auto create ( const GridView& grid_view, const std::shared_ptr log=get_logger(log_base) ) - -> std::shared_ptr> + -> std::shared_ptr> { const auto filename = config.get("file"); const auto dataset = config.get("dataset"); @@ -300,7 +456,7 @@ static auto create ( const GridView& grid_view, const std::shared_ptr log=get_logger(log_base) ) - -> std::shared_ptr> + -> std::shared_ptr> { log->trace("Creating interpolator from YAML cfg node: {}", node_name); const auto dim = GridView::Grid::dimension; diff --git a/dune/dorie/model/richards/parameterization/scaling.hh b/dune/dorie/model/richards/parameterization/scaling.hh index d65563a5937eec863b24811fde78441b9a81357d..5f5617142f7013e4b4f4ac1d95a4f0dcd322b365 100644 --- a/dune/dorie/model/richards/parameterization/scaling.hh +++ b/dune/dorie/model/richards/parameterization/scaling.hh @@ -56,7 +56,7 @@ protected: using return_t = ScalingFactors; //! Interpolator storage - std::vector>> _interpolators; + std::vector>> _interpolators; /// Grid view optionally used for instatiating interpolators const GridView& _grid_view; /// Logger of this instance diff --git a/dune/dorie/test/test-interpolators.cc b/dune/dorie/test/test-interpolators.cc index ff2e89c696a9d38e69bde5724f41c26ccd81749c..c34e8c38e3f79905720aa4441142edcce95d5ae0 100644 --- a/dune/dorie/test/test-interpolators.cc +++ b/dune/dorie/test/test-interpolators.cc @@ -12,6 +12,7 @@ #include #include #include +#include #include @@ -19,7 +20,7 @@ template struct InterpolatorTraits { static constexpr int dim = dimension; - using Domain = std::vector; + using Domain = Dune::FieldVector; using DF = double; using RF = double; }; @@ -39,11 +40,11 @@ void test_nearest_neighbor () std::vector shape(dim); std::fill(begin(shape), end(shape), 3); - std::vector extensions(dim); - std::fill(begin(extensions), end(extensions), 3.0); + Dune::FieldVector extensions(3.0); + // std::fill(begin(extensions), end(extensions), 3.0); - std::vector offset(dim); - std::fill(begin(offset), end(offset), 0.0); + Dune::FieldVector offset(0.0); + // std::fill(begin(offset), end(offset), 0.0); // build interpolator auto interp = Dune::Dorie::InterpolatorFactory> @@ -55,7 +56,7 @@ void test_nearest_neighbor () // check without offset using Dune::FloatCmp::eq; // floating-point comparison - std::vector> corners; + std::vector> corners; if constexpr (dim == 2) { corners.resize(4); corners[0] = {0.0, 0.0}; @@ -88,7 +89,7 @@ void test_nearest_neighbor () } // check with offset - std::fill(begin(offset), end(offset), -1.0); + std::fill(offset.begin(), offset.end(), -1.0); interp = Dune::Dorie::InterpolatorFactory> ::create("nearest", data, @@ -128,16 +129,130 @@ void test_nearest_neighbor () } } +/// Test the linear interpolator +template +void test_linear (); + +/// Test the linear interpolator in 2D +template<> +void test_linear<2> () +{ + constexpr int dim = 2; + std::vector data({0.0, 1.0, 1.0, 2.0}); + + std::vector shape(dim); + std::fill(begin(shape), end(shape), 2); + + Dune::FieldVector extensions(1.0); + Dune::FieldVector offset(0.0); + + // build interpolator + auto interp = Dune::Dorie::InterpolatorFactory> + ::create("linear", + data, + shape, + extensions, + offset); + + // check without offset + using Dune::FloatCmp::eq; // floating-point comparison + std::vector> points(5); + + points[0] = {0.0, 0.0}; + points[1] = {1.0, 0.0}; + points[2] = {0.0, 1.0}; + points[3] = {1.0, 1.0}; + points[4] = {0.75, 0.75}; + assert(eq(interp->evaluate(points[0]), 0.0)); + assert(eq(interp->evaluate(points[1]), 1.0)); + assert(eq(interp->evaluate(points[2]), 1.0)); + assert(eq(interp->evaluate(points[3]), 2.0)); + assert(eq(interp->evaluate(points[4]), 1.5)); + + + // check with offset + std::fill(offset.begin(), offset.end(), -1.0); + interp = Dune::Dorie::InterpolatorFactory> + ::create("linear", + data, + shape, + extensions, + offset); + + points.resize(2); + points[0] = {0.0, 0.0}; + points[1] = {-0.25, -0.25}; + assert(eq(interp->evaluate(points[0]), 2.0)); + assert(eq(interp->evaluate(points[1]), 1.5)); +} + +/// Test the linear interpolator in 3D +template<> +void test_linear<3> () +{ + constexpr int dim = 3; + std::vector data({0.0, 0.0, 0.0, 0.0, + 1.0, 1.0, 2.0, 2.0}); + + std::vector shape(dim); + std::fill(begin(shape), end(shape), 2); + + Dune::FieldVector extensions(1.0); + Dune::FieldVector offset(0.0); + + // build interpolator + auto interp = Dune::Dorie::InterpolatorFactory> + ::create("linear", + data, + shape, + extensions, + offset); + + // check without offset + using Dune::FloatCmp::eq; // floating-point comparison + std::vector> points(5); + + points[0] = {0.0, 0.0, 0.0}; + points[1] = {1.0, 0.0, 0.0}; + points[2] = {0.0, 0.0, 1.0}; + points[3] = {1.0, 1.0, 1.0}; + points[4] = {0.5, 0.5, 1.0}; + assert(eq(interp->evaluate(points[0]), 0.0)); + assert(eq(interp->evaluate(points[1]), 0.0)); + assert(eq(interp->evaluate(points[2]), 1.0)); + assert(eq(interp->evaluate(points[3]), 2.0)); + assert(eq(interp->evaluate(points[4]), 1.5)); + + // check with offset + std::fill(offset.begin(), offset.end(), -1.0); + interp = Dune::Dorie::InterpolatorFactory> + ::create("linear", + data, + shape, + extensions, + offset); + + points.resize(2); + points[0] = {0.0, 0.0, 0.0}; + points[1] = {-0.5, -0.5, 0.0}; + assert(eq(interp->evaluate(points[0]), 2.0)); + assert(eq(interp->evaluate(points[1]), 1.5)); +} + int main (int argc, char** argv) { try{ // initialize MPI if needed auto& helper = Dune::MPIHelper::instance(argc, argv); - Dune::Dorie::create_base_logger(helper); + auto log = Dune::Dorie::create_base_logger(helper); + log->set_level(spdlog::level::trace); // test the NearestNeighbor interpolator test_nearest_neighbor<2>(); test_nearest_neighbor<3>(); + + test_linear<2>(); + test_linear<3>(); } catch (Dune::Exception &e) {