Commit a2cc481b authored by Lukas Riedel's avatar Lukas Riedel Committed by Santiago Ospina De Los Ríos

Add option to compute IC from stationary problem

In the Richards model, add the IC type option 'stationary'. In this
mode, the stationary problem is solved by assembling a separate Newton
solver and applying it to the spatial grid operator only.
parent 67d678fd
......@@ -6,7 +6,6 @@ doc/manual/config-file.rst
doc/default_files/config.ini
doc/cookbook/1-infiltration-sand/config.ini
python/dorie/wrapper/pf_from_file.py
python/dorie/wrapper/test_dorie.py
python/dorie/dorie/cli/cmds.py
test/maps/cell_ids.h5
......
......@@ -32,6 +32,7 @@
* Parameter XML files support tags for version changes !171
* Unit test for Parameter XML file parser !171
* CMake option to enable code coverage flags on all targets !172
* Steady state initial condition in Richards model !176
* Changes to config file parameters listed per version in user docs !175
### Changed
......
......@@ -65,6 +65,15 @@ Boundary conditions of the same boundary segment must not overlap in time.
If you wish one condition to follow exactly after the other, make sure that
the end of the first interval matches the start of the second interval exactly.
Boundary conditions whose time intervals start and end with the same time, or
lie completely outside the simulation time interval, are ignored. The only
exception occurs when computing an :doc:`initial condition </manual/initial>`
from the stationary problem. In this case, the earliest boundary condition
which overlaps with the simulation start time is chosen. In particular, its
time interval may only consist of the simulation start time stamp. See the docs
of the :ref:`stationary initial condition type <initial-type-stationary>` for
an example.
.. note:: DORiE automatically adjusts its time steps to match a change in
boundary conditions.
......
Initial Conditions
==================
.. contents::
:depth: 3
:local:
DORiE computes transient solutions, and hence needs a solution from which it
can start. There are multiple ways of specifying an initial condition.
All have in common that the data provided must be interpolated on the
......@@ -27,65 +31,112 @@ Initial condition input is controlled entirely via the
Input Types
-----------
This is an overview of all input types for initial conditions.
They are controlled by the ``initial.type`` key and available for every model.
They are controlled by the ``initial.type`` key and shared between all models
unless otherwise specified.
Analytic
^^^^^^^^
* ``type = analytic``
An analytic function :math:`f(\vec{p})` which depends on the physical
position :math:`\vec{p}`. The function must be defined via the key
``initial.equation``. For parsing the input expression, we use muparser_
which supports a set of common mathematical functions. Additionally, the
following variables can be used:
Available variables:
* ``x``: X-coordinate :math:`p_1 \, [\mathrm{m}]`.
* ``y``: Y-coordinate :math:`p_2 \, [\mathrm{m}]`.
* ``z``: Z-coordinate :math:`p_3 \, [\mathrm{m}]` (only in 3D).
* ``h``: Height above origin. Synonymous to ``y`` in 2D and ``z`` in 3D.
* ``pi``: Mathematical constant :math:`\pi`.
* ``dim``: Number of physical dimensions.
.. tip::
Assuming the target quantity is the matric head (see
:ref:`initial-transformation`), typical initial conditions for a
Richards model are
* Hydrostatic equilibrium: A vertical gradient of :math:`-1` and a
fixed value ``<v>`` at height :math:`h = 0 \, \mathrm{m}`::
initial.equation = -h + <v>
* Gravity flow: Constant value.
.. tip::
The expression for a gaussian pulse of solute concentration centered at
:math:`\vec{p} = [0.5, 0.5]^T \, \mathrm{m}` and variance of
:math:`\sigma^2 = \left( 0.1 \, \mathrm{m} \right)^2` is::
initial.equation = <m> * exp(-((x-0.5)^2+(y-0.5)^2)/(2.*0.1^2)) / (2*pi*0.1^2)
where ``<m>`` is the total solute mass of the pulse
:math:`m_s \, [\text{kg}]`.
.. object:: Analytic
Dataset
^^^^^^^
* ``type = analytic``
* ``type = data``
An analytic function :math:`f(\vec{p})` which depends on the physical
position :math:`\vec{p}`. The function must be defined via the key
``initial.equation``. For parsing the input expression, we use muparser_
which supports a set of common mathematical functions. Additionally, the
following variables can be used:
Load the initial condition from a data file ``initial.file`` by opening the
dataset ``initial.dataset`` inside this file. The data is interpreted as
function :math:`f(\mathbf{p})` of the physical position :math:`\mathbf{p}`
using one of the available :ref:`interpolators <man-interpolators>`, which
can be chosen using the setting ``initial.interpolation``. The input data
is automatically streched to match the grid extensions.
Available variables:
* ``x``: X-coordinate :math:`p_1 \, [\mathrm{m}]`.
* ``y``: Y-coordinate :math:`p_2 \, [\mathrm{m}]`.
* ``z``: Z-coordinate :math:`p_3 \, [\mathrm{m}]` (only in 3D).
* ``h``: Height above origin. Synonymous to ``y`` in 2D and ``z`` in 3D.
* ``pi``: Mathematical constant :math:`\pi`.
* ``dim``: Number of physical dimensions.
.. note:: For ``FEorder > 0``, linear interpolation is recommended.
.. tip::
Assuming the target quantity is the matric head (see
:ref:`initial-transformation`), typical initial conditions for a
Richards model are
Supported file extensions:
* Hydrostatic equilibrium: A vertical gradient of :math:`-1` and a
fixed value ``<v>`` at height :math:`h = 0 \, \mathrm{m}`::
* ``.h5``: H5_ data file. ``initial.dataset`` may be a file-internal path
to the target dataset.
initial.equation = -h + <v>
.. _initial-type-stationary:
* Gravity flow: Constant value.
Stationary (Richards only)
^^^^^^^^^^^^^^^^^^^^^^^^^^
.. tip::
The expression for a gaussian pulse of solute concentration centered at
:math:`\vec{p} = [0.5, 0.5]^T \, \mathrm{m}` and variance of
:math:`\sigma^2 = \left( 0.1 \, \mathrm{m} \right)^2` is::
* ``type = stationary``
initial.equation = <m> * exp(-((x-0.5)^2+(y-0.5)^2)/(2.*0.1^2)) / (2*pi*0.1^2)
Compute the stationary part of the Richards equation,
where ``<m>`` is the total solute mass of the pulse
:math:`m_s \, [\text{kg}]`.
.. math::
.. object:: Dataset
- \nabla \cdot \left[ K \left[ \nabla h_m - \mathbf{\hat{g}} \right]
\right] = 0 \ ,
* ``type = data``
and use the solution as initial condition for the transient simulation. The
boundary conditions at the simulation start time are used to solve the
problem. No further input is needed, but the boundary condition file might
need to be adjusted.
Load the initial condition from a data file ``initial.file`` by opening the
dataset ``initial.dataset`` inside this file. The data is interpreted as
function :math:`f(\mathbf{p})` of the physical position :math:`\mathbf{p}`
using one of the available :ref:`interpolators <man-interpolators>`, which
can be chosen using the setting ``initial.interpolation``. The input data
is automatically streched to match the grid extensions.
.. note:: As this computes the actual solution, no transformation is applied.
The config file key ``initial.quantity`` is ignored.
.. note:: For ``FEorder > 0``, linear interpolation is recommended.
.. tip:: If you want the simulation to start with different boundary
conditions than those used for computing the stationary initial
condition, specify boundary conditions that *end* with the simulation
*start time* and start at the same time or earlier, like so (assuming
the simulation starts at :math:`t=0\,\text{s}`):
Supported file extensions:
.. code:: yaml
* ``.h5``: H5_ data file. ``initial.dataset`` may be a file-internal path
to the target dataset.
upper:
index: 1
conditions:
# Applied for computing the initial condition
for_ic:
type: Neumann
value: -5e-6
time: [0.0, 0.0]
# Applied during the entire transient simulation
for_sim:
type: Dirichlet
value: -6
time: 0.0
.. _initial-transformation:
......
......@@ -37,6 +37,17 @@ namespace Dorie {
* \f$ t_2 \f$, a new time step must start, during which only boundary
* condition \f$ 2 \f$ is applied.
*
* After initializing this class, the first boundary condition is selected
* whose begin time is the same or earlier than the simulation start time
* and whose end time is the same or later than the simulation start time.
* This is in contrast to the selection mechanism of set_time(), where the
* cached BC must have an end time which is later than the current simulation
* time. Calling set_time() with the simulation start time might lead to a
* different BC being cached. This allows for separate boundary conditions
* right at the start of the simulation, which is useful for computing an
* initial condition from the stationary problem before starting the actual
* simulation.
*
* This class is intended to be used as follows:
* 1. Build storage by constructing the class.
* 2. At the beginning of a simulation time step, call set_time() with the
......@@ -108,7 +119,7 @@ public:
_log(log),
_bcf(_log),
_bc_storage(build_bc_storage(config, boundary_index_map)),
_bc_cache(build_bc_cache(_bc_storage))
_bc_cache(build_bc_cache(config, _bc_storage))
{ }
/// Return the maximum time step allowed with the current BCs
......@@ -195,15 +206,18 @@ public:
if (time_matches_interval(time, _bc_cache[i]))
continue;
// find new BC
// Find new BC
// NOTE: Use reverse iterators because we want to find the last
// element in the list for which the time matches. The list
// is sorted by begin *and* end times.
const auto& bc_list = _bc_storage[i];
const auto it = std::find_if(begin(bc_list),
end(bc_list),
const auto it = std::find_if(rbegin(bc_list),
rend(bc_list),
std::bind(time_matches_interval,
time,
std::placeholders::_1));
if (it == end(bc_list)) {
if (it == rend(bc_list)) {
_log->error("Time '{}' exceeds boundary conditions of "
"boundary segment index '{}'",
time, i);
......@@ -287,8 +301,14 @@ private:
std::sort(begin(list), end(list),
[](const auto bc1, const auto bc2)
{
return bc1->time_interval().begin
< bc2->time_interval().begin;
const auto begin1 = bc1->time_interval().begin;
const auto end1 = bc1->time_interval().end;
const auto begin2 = bc2->time_interval().begin;
const auto end2 = bc2->time_interval().end;
return begin1 < begin2
or (FloatCmp::eq(begin1, begin2)
and end1 < end2);
});
// insert an initial boundary condition, if needed
......@@ -374,25 +394,62 @@ private:
}
/// Build an initialized cache from a boundary condition storage
/** This caches the primary BC of every entry in the storage
* and checks if all storage lists are non-empty
/** This caches the first BC (according to the simulation start time) of
* every entry in the storage and checks if all storage lists are
* non-empty.
*
* \param config The configuration file tree
* \param bc_storage The storage to initialize the cache from
* \return A list caching the first BC for every BC index
*/
BCList build_bc_cache (const BCStorage& bc_storage) const
BCList build_bc_cache (const ParameterTree& config,
const BCStorage& bc_storage) const
{
// Retrieve the simulation start time (this is where we should be)
const auto t_start = config.get<double>("time.start");
// Define lambda for matching time and time interval of a BC
// NOTE: This allows boundary conditions with the same time as
// begin and end times, useful for the stationary IC.
auto time_matches_interval = [](const RF time,
const std::shared_ptr<BC> bc)
{
const auto& time_int = bc->time_interval();
return FloatCmp::ge(time, time_int.begin)
&& FloatCmp::le(time, time_int.end);
};
// Iterate over all boundary segments
BCList bc_cache(bc_storage.size());
for (size_t i = 0; i < bc_storage.size(); ++i)
for (size_t i = 0; i < _bc_storage.size(); ++i)
{
if (bc_storage[i].empty()) {
// Retrieve the BC list for this boundary segment
const auto& bc_list = _bc_storage[i];
if (bc_list.empty()) {
_log->error("Boundary condition storage has no entry for "
"boundary segment index '{}'",
i);
DUNE_THROW(InvalidStateException, "Encountered empty BC"
"storage list");
}
bc_cache[i] = bc_storage[i].front();
// Find first BC matching the start time
const auto it = std::find_if(begin(bc_list),
end(bc_list),
std::bind(time_matches_interval,
t_start,
std::placeholders::_1));
if (it == end(bc_list)) {
_log->error("Start time '{}' exceeds boundary conditions of "
"boundary segment index '{}'",
t_start, i);
DUNE_THROW(InvalidStateException, "No BCs defined at "
"simulation start time");
}
// Update cache
bc_cache[i] = *it;
}
return bc_cache;
......
......@@ -26,12 +26,22 @@ private:
public:
/// Create an initial condition for the Richards solver
/**
* \note In case the user wants the initial condition to be computed from
* the stationary problem, an empty pointer is returned.
*/
static std::unique_ptr<InitialCondition<T>> create (
const Dune::ParameterTree& config,
const typename T::GridView& grid_view,
const std::shared_ptr<spdlog::logger> log=get_logger(log_richards)
)
{
if (config.get<std::string>("initial.type") == "stationary") {
// NOTE: Return empty pointer to communicate that the IC is
// computed internally.
return std::unique_ptr<InitialCondition<T>>();
}
const auto quantity = config.get<std::string>("initial.quantity");
// matric head: no transformation necessary
if (quantity == "matricHead") {
......
......@@ -84,9 +84,27 @@ ModelRichards<Traits>::ModelRichards (
tlop = std::make_unique<TLOP>(fparam);
}
// --- Solution Vectors and Initial Condition ---
// --- Solution Vector --- //
u = std::make_shared<U>(*gfs,.0);
Dune::PDELab::interpolate(*finitial,*gfs,*u);
// --- Operator Setup --- //
operator_setup();
// --- Initial Condition --- //
// Initial condition exists. Just interpolate onto DOFs
if (finitial) {
Dune::PDELab::interpolate(*finitial, *gfs, *u);
}
// Comput initial condition from stationary problem
else {
try {
compute_initial_condition();
}
catch (...) {
this->_log->error("Computing the initial condition failed");
throw;
}
}
// --- Writer Setup --- //
if (output_policy() != OutputPolicy::None)
......@@ -115,9 +133,6 @@ ModelRichards<Traits>::ModelRichards (
adaptivity = adaptivity_fac.create();
}
// --- Operator Setup ---
operator_setup();
this->_log->info("Setup complete");
}
......@@ -367,5 +382,68 @@ int ModelRichards<Traits>::newton_iterations_max (const RF dt) const
+ newton_it_max );
}
template<typename Traits>
void ModelRichards<Traits>::compute_initial_condition()
{
this->_log->info("Computing stationary problem for initial condition");
// Build new Newton solver
using Newton = Dune::PDELab::Newton<GO0, LS, U>;
Newton newton(*go0, *ls);
newton.setParameters(inifile.sub("NewtonParameters"));
newton.setVerbosityLevel(0);
newton.setMaxIterations(50); // Just a large number to ensure convergence
// Set time in local operator (for BCs)
slop->setTime(this->_time);
// Use hydrostatic equilibrium for an initial guess
using T = typename Traits::BaseTraits;
auto init = std::make_unique<InitialConditionAnalytic<T>>(
this->_log, gv, "-h"
);
Dune::PDELab::interpolate(*init, *gfs, *u);
// Lambda for applying the solver:
// Return true if the solver converged, let DUNE errors pass
auto solve_stationary = [&newton](U& u) -> bool
{
try {
newton.apply(u);
return true;
}
// Ignore Dune errors. The calling function will try again.
catch (Exception&) {
// pass
}
return false;
};
// Try solving
auto solved = solve_stationary(*u);
// Use static initial guesses if that did not work
std::vector<double> init_values({0.0, -1.0, -10.0});
for(size_t i = 0;
not solved and i < init_values.size();
++i)
{
this->_log->warn("Computing the initial condition failed. "
"Retrying with different initial guess...");
u = std::make_shared<U>(*gfs, init_values[i]);
solved = solve_stationary(*u);
}
if (not solved) {
DUNE_THROW(Exception, "No initial guess worked for computing the IC");
}
// Report times
const auto& result = newton.result();
this->_log->trace(" Matrix assembly: {:.2e}s", result.assembler_time);
this->_log->trace(" Linear solver: {:.2e}s", result.linear_solver_time);
}
} // namespace Dorie
} // namespace Dune
......@@ -637,6 +637,9 @@ protected:
*/
int newton_iterations_max (const RF dt) const;
/// Compute the initial condition based on the stationary problem
void compute_initial_condition ();
private:
#ifdef GTEST
......
......@@ -3,7 +3,6 @@ message(STATUS "Handling unit tests")
# copy ALL the config files and stuff
file(COPY scaling.h5
scaling.yml
test-boundary-condition.yml
test-param-richards-scaled.yml
test-param-richards-unscaled.yml
test-param-transport.yml
......
......@@ -208,8 +208,19 @@ TEST_F(BCManager, Neumann)
const auto coords = is.geometry().center();
if (Dune::FloatCmp::eq(coords[1], 1.0))
{
using namespace Dune::Dorie;
// Check initial cache (before setting time)
const auto bc = _manager->bc(is);
EXPECT_EQ(bc->name(), "neumann_initial");
EXPECT_EQ(bc->type(), Operator::BCType::Neumann);
EXPECT_DOUBLE_EQ(bc->evaluate(0.0), 0.0);
auto &time_int = bc->time_interval();
EXPECT_DOUBLE_EQ(time_int.begin, 0.0);
EXPECT_DOUBLE_EQ(time_int.end, 0.0);
// Check BCs after setting times
for (auto time: {0.0, 1.0, 1.5, 1.9}) {
using namespace Dune::Dorie;
_manager->set_time(time);
const auto bc = _manager->bc(is);
EXPECT_EQ(bc->name(), "neumann");
......@@ -243,6 +254,14 @@ TEST_F(BCManager, Dirichlet)
if (Dune::FloatCmp::eq(coords[1], 0.0))
{
using namespace Dune::Dorie;
// Check initial cache (before setting time)
const auto bc = _manager->bc(is);
EXPECT_EQ(bc->name(), "default");
EXPECT_EQ(bc->type(), Operator::BCType::Neumann);
EXPECT_DOUBLE_EQ(bc->evaluate(0.0), 0.0);
// Check BCs after setting times
for (auto time: {0.5, 1.4}) {
_manager->set_time(time);
const auto bc = _manager->bc(is);
......@@ -285,11 +304,23 @@ TEST_F(BCManager, Outflow)
const auto coords = is.geometry().center();
if (Dune::FloatCmp::eq(coords[0], 0.0))
{
using namespace Dune::Dorie;
// Check initial cache (before setting time)
const auto bc = _manager->bc(is);
EXPECT_EQ(bc->name(), "outflow_initial");
EXPECT_EQ(bc->type(), Operator::BCType::Outflow);
EXPECT_DOUBLE_EQ(bc->evaluate(0.0), 1.0);
auto &time_int = bc->time_interval();
EXPECT_DOUBLE_EQ(time_int.begin, -1.0);
EXPECT_DOUBLE_EQ(time_int.end, 0.0);
// Check BCs after setting times
for (auto time: {0.0, 1.0, 1.5, 1.9}) {
_manager->set_time(time);
const auto bc = _manager->bc(is);
EXPECT_EQ(bc->name(), "outflow");
EXPECT_EQ(bc->type(), Dune::Dorie::Operator::BCType::Outflow);
EXPECT_EQ(bc->type(), Operator::BCType::Outflow);
EXPECT_DOUBLE_EQ(bc->evaluate(time), 0.0);
auto &time_int = bc->time_interval();
EXPECT_DOUBLE_EQ(time_int.begin, 0.0);
......
include ${CMAKE_BINARY_DIR}/doc/default_files/config.ini
__name = bc
_asset_path = "${PROJECT_SOURCE_DIR}/test"
_asset_path = "${CMAKE_CURRENT_LIST_DIR}"
[grid]
dimensions = 2
......@@ -27,6 +27,6 @@ output.fileName = bc | unique name
output.outputPath = bc | unique name
output.logLevel = warning
boundary.file = test-boundary-condition.yml
boundary.file = {_asset_path}/test-boundary-condition.yml
time.end = 2.0
upper:
index: 1
conditions:
neumann_initial:
type: Neumann
value: 0.0
time: [0.0, 0.0]
neumann:
type: Neumann
value: 1.0
......@@ -23,6 +27,10 @@ left:
type: Outflow
value: 0.0
time: 0.0
outflow_initial:
type: Outflow
value: 1.0
time: [-1.0, 0.0]
default:
type: Neumann
......
......@@ -25,9 +25,9 @@ def evaluate(iniinfo,runtime):
fe_order = int(iniinfo["richards.numerics.FEorder"])
do_flux_rt = iniinfo["richards.fluxReconstruction.enable"]
do_flux_grad = fe_order > 0
is_stationary = (iniinfo["richards.initial.type"] == "stationary")
if do_flux_rt:
flux_rt_key = "flux_RT" + str(max(0,fe_order-1))
# get parameter data
param_data = read_yml(iniinfo["richards.parameters.file"])
......@@ -123,6 +123,8 @@ def evaluate(iniinfo,runtime):
# Check matric head
head_tol = 1E-5 if not "_ode.head_abstol" in iniinfo else float(iniinfo["_ode.head_abstol"])
if is_stationary:
head_tol = head_tol * float(iniinfo["_ode.head_eps_stat"])
passed = bool(l2_head < head_tol)
# Check flux from gradient
......@@ -138,6 +140,8 @@ def evaluate(iniinfo,runtime):
plt.close()
flux_tol = abs(influx) * 1e-5 if not "_ode.flux_abstol" in iniinfo else float(iniinfo["_ode.flux_abstol"])
if is_stationary:
flux_tol = flux_tol * float(iniinfo["_ode.flux_eps_stat"])
passed = passed and all([l < flux_tol for l in l2_flux])
# Check reconstructed flux
......@@ -153,6 +157,8 @@ def evaluate(iniinfo,runtime):
plt.close()
flux_rt_tol = 1e-13 if not "_ode.flux_rt_abstol" in iniinfo else float(iniinfo["_ode.flux_rt_abstol"])
if is_stationary:
flux_rt_tol = flux_rt_tol * float(iniinfo["_ode.flux_rt_eps_stat"])
passed = passed and all([l < flux_rt_tol for l in l2_flux_rt])
......
configure_file(pf_from_file.py.in ${CMAKE_CURRENT_SOURCE_DIR}/pf_from_file.py)
configure_file(test_dorie.py.in ${CMAKE_CURRENT_SOURCE_DIR}/test_dorie.py)
\ No newline at end of file
......@@ -4,15 +4,23 @@ from __future__ import absolute_import
import sys
import os
import argparse
from dune.testtools.wrapper.argumentparser import get_args
from dune.testtools.parser import parse_ini_file
from dorie.testtools.wrapper import test_dorie
DORIE_WRAPPER = "dorie"
def get_args():
parser = argparse.ArgumentParser()
parser.add_argument('-i', '--ini', help='The inifile', required=True)
args, _ = parser.parse_known_args()
return args