Commit 99e71e5d authored by Santiago Ospina De Los Ríos's avatar Santiago Ospina De Los Ríos

Merge branch '121-add-linear-interpolator-and-use-it-for-scaling-fields' into 'master'

Resolve "Add linear interpolator and use it for scaling fields"

Closes #121

See merge request !145
parents 9e404b86 3b3f6baa
......@@ -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
......
......@@ -153,8 +153,8 @@ adding an empty line, make text **bold** or ``monospaced``.
<parameter name="interpolation">
<definition> Interpolation type used for the data (``data`` type only).
</definition>
<values> nearest </values>
<suggestion> nearest </suggestion>
<values> nearest, linear </values>
<suggestion> linear </suggestion>
</parameter>
</category>
......
......@@ -170,7 +170,7 @@ adding an empty line, make text **bold** or ``monospaced``.
<parameter name="interpolation">
<definition> Interpolation type used for the data (``data`` type only).
</definition>
<values> nearest </values>
<values> nearest, linear </values>
<suggestion> nearest </suggestion>
</parameter>
</category>
......
......@@ -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::
......
......@@ -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 <man-config-file>`.
:doc:`Configuration File <config-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
......
.. _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 <parameter-file>`, and when loading an initial condition
from data with the options given in the
:doc:`Configuration File <config-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 <man-parameter_scaling>` 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.
......@@ -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 <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/
......@@ -61,7 +61,7 @@ public:
}
private:
std::shared_ptr<Interpolator<RF, T>> _interpolator;
std::shared_ptr<Interpolator<RF, dim>> _interpolator;
};
} // namespace Dorie
......
This diff is collapsed.
......@@ -56,7 +56,7 @@ protected:
using return_t = ScalingFactors<RF>;
//! Interpolator storage
std::vector<std::shared_ptr<Interpolator<RF, Traits>>> _interpolators;
std::vector<std::shared_ptr<Interpolator<RF, Traits::dim>>> _interpolators;
/// Grid view optionally used for instatiating interpolators
const GridView& _grid_view;
/// Logger of this instance
......
......@@ -12,6 +12,7 @@
#include <dune/common/exceptions.hh>
#include <dune/common/float_cmp.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/fvector.hh>
#include <dune/dorie/common/interpolator.hh>
......@@ -19,7 +20,7 @@ template<int dimension>
struct InterpolatorTraits
{
static constexpr int dim = dimension;
using Domain = std::vector<double>;
using Domain = Dune::FieldVector<double, dim>;
using DF = double;
using RF = double;
};
......@@ -39,11 +40,11 @@ void test_nearest_neighbor ()
std::vector<size_t> shape(dim);
std::fill(begin(shape), end(shape), 3);
std::vector<double> extensions(dim);
std::fill(begin(extensions), end(extensions), 3.0);
Dune::FieldVector<double, dim> extensions(3.0);
// std::fill(begin(extensions), end(extensions), 3.0);
std::vector<double> offset(dim);
std::fill(begin(offset), end(offset), 0.0);
Dune::FieldVector<double, dim> offset(0.0);
// std::fill(begin(offset), end(offset), 0.0);
// build interpolator
auto interp = Dune::Dorie::InterpolatorFactory<InterpolatorTraits<dim>>
......@@ -55,7 +56,7 @@ void test_nearest_neighbor ()
// check without offset
using Dune::FloatCmp::eq; // floating-point comparison
std::vector<std::vector<double>> corners;
std::vector<Dune::FieldVector<double, dim>> 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<InterpolatorTraits<dim>>
::create("nearest",
data,
......@@ -128,16 +129,130 @@ void test_nearest_neighbor ()
}
}
/// Test the linear interpolator
template<int dim>
void test_linear ();
/// Test the linear interpolator in 2D
template<>
void test_linear<2> ()
{
constexpr int dim = 2;
std::vector<double> data({0.0, 1.0, 1.0, 2.0});
std::vector<size_t> shape(dim);
std::fill(begin(shape), end(shape), 2);
Dune::FieldVector<double, dim> extensions(1.0);
Dune::FieldVector<double, dim> offset(0.0);
// build interpolator
auto interp = Dune::Dorie::InterpolatorFactory<InterpolatorTraits<dim>>
::create("linear",
data,
shape,
extensions,
offset);
// check without offset
using Dune::FloatCmp::eq; // floating-point comparison
std::vector<Dune::FieldVector<double, dim>> 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<InterpolatorTraits<dim>>
::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<double> data({0.0, 0.0, 0.0, 0.0,
1.0, 1.0, 2.0, 2.0});
std::vector<size_t> shape(dim);
std::fill(begin(shape), end(shape), 2);
Dune::FieldVector<double, dim> extensions(1.0);
Dune::FieldVector<double, dim> offset(0.0);
// build interpolator
auto interp = Dune::Dorie::InterpolatorFactory<InterpolatorTraits<dim>>
::create("linear",
data,
shape,
extensions,
offset);
// check without offset
using Dune::FloatCmp::eq; // floating-point comparison
std::vector<Dune::FieldVector<double, dim>> 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<InterpolatorTraits<dim>>
::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) {
......
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