...
  View open merge request
Commits (99)
......@@ -95,6 +95,11 @@ function(dorie_add_unit_test)
muparser::muparser hdf5 yaml-cpp spdlog)
# add_coverage_links(${UNIT_TEST_TARGET})
target_compile_definitions(${UNIT_TEST_TARGET}
PUBLIC
GTEST
)
if (UNIT_TEST_CUSTOM_MAIN)
target_link_libraries(${UNIT_TEST_TARGET} gtest)
else ()
......@@ -232,6 +237,11 @@ function(dorie_add_metaini_test)
# add_coverage_links(${created_targets})
add_dependencies(build_unit_tests ${created_targets})
target_compile_definitions(${created_targets}
PUBLIC
GTEST
)
if (SYSTEM_TEST_CUSTOM_MAIN)
target_link_libraries(${created_targets} gtest)
else ()
......
spatial_resolution_north 0
spatial_resolution_south 0
spatial_resolution_west -1
spatial_resolution_east -1
number_BC_change_times 1
0 neumann -5.55e-6 dirichlet 0
\ No newline at end of file
spatial_resolution_north 2 0.25 0.75
spatial_resolution_south -1
spatial_resolution_west -1
spatial_resolution_east -1
number_BC_change_times 2
0 neumann 0 dirichlet 1 neumann 0
1E3 neumann 0 neumann 0 neumann 0
spatial_resolution_north_we 0
spatial_resolution_north_fb 0
spatial_resolution_south_we 0
spatial_resolution_south_fb 0
spatial_resolution_west_sn -1
spatial_resolution_west_fb -1
spatial_resolution_east_sn -1
spatial_resolution_east_fb -1
spatial_resolution_front_sn -1
spatial_resolution_front_we -1
spatial_resolution_back_sn -1
spatial_resolution_back_we -1
number_BC_change_times 1
0 neumann -5.55e-6 dirichlet 0
\ No newline at end of file
spatial_resolution_north_we 2 0.25 0.75
spatial_resolution_north_fb 2 0.25 0.75
spatial_resolution_south_we -1
spatial_resolution_south_fb -1
spatial_resolution_west_sn -1
spatial_resolution_west_fb -1
spatial_resolution_east_sn -1
spatial_resolution_east_fb -1
spatial_resolution_front_sn -1
spatial_resolution_front_we -1
spatial_resolution_back_sn -1
spatial_resolution_back_we -1
number_BC_change_times 2
0 neumann 0 neumann 0 neumann 0 neumann 0 dirichlet 1 neumann 0 neumann 0 neumann 0 neumann 0
1E3 neumann 0 neumann 0 neumann 0 neumann 0 neumann 0 neumann 0 neumann 0 neumann 0 neumann 0
\ No newline at end of file
......@@ -117,9 +117,10 @@ function(create_default_config INPUT OUTPUT SOURCE_DIR CSS)
endfunction()
# copy BC & parameter files
file(COPY 2d_infiltr.bcdat 3d_infiltr.bcdat
2d_solute.bcdat 3d_solute.bcdat
richards_param.yml transport_param.yml
file(COPY infiltration.yml
solute.yml
richards_param.yml
transport_param.yml
DESTINATION .)
# Random field generator
......@@ -167,4 +168,4 @@ create_default_config(
"config.ini;parameters.html;parameters.rst"
${PROJECT_SOURCE_DIR}/dune/dorie
${CMAKE_CURRENT_SOURCE_DIR}/parameters.css
)
\ No newline at end of file
)
upper:
index: 1
conditions:
heavy_rainfall:
type: Neumann
time: 0
value: -5.55e-6
horizontal_projection: true
lower:
index: 0
conditions:
water_table:
type: Dirichlet
time: 0
value: 0
default:
type: Neumann
value: 0
......@@ -90,24 +90,10 @@ adding an empty line, make text **bold** or ``monospaced``.
<category name="boundary">
<parameter name="file">
<definition> Path to the boundary condition file. </definition>
<definition> Path to the YAML boundary condition data file.
</definition>
<values> path </values>
</parameter>
<parameter name="fileType">
<definition> Type of spatial segmentation of the boundaries specified in the BC file </definition>
<values> rectangularGrid </values>
<suggestion> rectangularGrid </suggestion>
<comment> Choose type of boundary segmentation: rectangularGrid </comment>
</parameter>
<parameter name="interpolateBCvalues">
<definition> Whether to interpolate between the boundary conditions
at different times linearly (``true``) or not at all (``false``). May require
different boundary condition files. </definition>
<values> true, false </values>
<suggestion> false </suggestion>
</parameter>
</category>
<category name="initial">
......
upper:
index: 1
conditions:
tracer:
type: Dirichlet
time: 0
value: 1
concentration_type: water_phase
default:
type: Outflow
value: 0
......@@ -101,30 +101,9 @@ adding an empty line, make text **bold** or ``monospaced``.
<category name="boundary">
<parameter name="file">
<definition> Path to the boundary condition file. </definition>
<values> path </values>
</parameter>
<parameter name="fileType">
<definition> Type of spatial segmentation of the boundaries specified in the BC file </definition>
<values> rectangularGrid </values>
<suggestion> rectangularGrid </suggestion>
<comment> Choose type of boundary segmentation: rectangularGrid </comment>
</parameter>
<parameter name="interpolateBCvalues">
<definition> Whether to interpolate between the boundary conditions
at different times linearly (``true``) or not at all (``false``). May require
different boundary condition files. </definition>
<values> true, false </values>
<suggestion> false </suggestion>
</parameter>
<parameter name="dirichletMode">
<definition> Type of the input value for the dirichlet condition.
<definition> Path to the YAML boundary condition data file.
</definition>
<values> soluteConcentration, totalSolute </values>
<suggestion> soluteConcentration </suggestion>
<values> path </values>
</parameter>
</category>
......
spatial_resolution_north 2 0.4 0.8
spatial_resolution_south 0
spatial_resolution_west -1
spatial_resolution_east -1
number_BC_change_times 1
0 neumann 0.0 neumann -5.55e-6 neumann 0.0 dirichlet 0
\ No newline at end of file
spatial_resolution_north_we 2 0.4 0.6
spatial_resolution_north_fb 2 0.4 0.6
spatial_resolution_south_we 0
spatial_resolution_south_fb 0
spatial_resolution_west_sn -1
spatial_resolution_west_fb -1
spatial_resolution_east_sn -1
spatial_resolution_east_fb -1
spatial_resolution_front_sn -1
spatial_resolution_front_we -1
spatial_resolution_back_sn -1
spatial_resolution_back_we -1
number_BC_change_times 1
0 neumann 0.0 neumann 0.0 neumann 0.0 neumann 0.0 neumann -5.55e-6 neumann 0.0 neumann 0.0 neumann 0.0 neumann 0.0 dirichlet 0
\ No newline at end of file
......@@ -25,9 +25,9 @@ Templates for the required files can conveniently be generated with the
Specifies the soil parameterization type(s) and parameters.
Template: ``param.yml``
* :doc:`Boundary Condition Data File </manual/bcfile>` (``.bcdat``):
* :doc:`Boundary Condition Data File </manual/bcfile>` (``.yml``):
Specifies the boundary segmentation and the boundary conditions.
Templates: ``2d_infiltr.bcdat``, ``3d_infiltr.bcdat``
Templates: ``infiltration.yml``, ``solute.yml``.
* :ref:`GMSH Mesh File <man-grid_gmsh>` (optional): Mesh with possibly
complicated geometry to build an unstructured grid from.
......
This diff is collapsed.
......@@ -68,6 +68,7 @@ successful. You can deactivate the ``venv`` at any time by executing
You are now ready to use the DORiE Command Line Interface!
.. _man-cli-commands:
CLI Commands
============
......
#ifndef DUNE_DORIE_COMMON_BOUNDARY_CONDITION_BASE_HH
#define DUNE_DORIE_COMMON_BOUNDARY_CONDITION_BASE_HH
#include <string>
#include <string_view>
#include <optional>
#include <memory>
#include <dune/common/exceptions.hh>
#include <dune/dorie/common/logging.hh>
#include <dune/dorie/common/typedefs.hh>
#include <dune/common/float_cmp.hh>
namespace Dune {
namespace Dorie {
/// A generic boundary condition storage
/** Stores the boundary condition time interval and values,
* and interpolates them linearly
*
* \tparam RF Type of the boundary codition value
* \ingroup Boundary
*/
template<typename RF>
class BoundaryCondition
{
public:
/// The type used for identifying time
using Time = double;
protected:
/// The time interval in which the BC is applied
const TimeInterval<Time> _time_interval;
/// The starting (or overall) value
const RF _value_begin;
/// The optional final value
const std::optional<RF> _value_end{};
/// The name given by the user in the data file
const std::string _name;
/// Check if the given time is inside the time interval
/** \param time The point in time to check
* \return True if the time is inside the interval of this BC
*/
bool in_interval (const Time time) const
{
const auto& t_begin = _time_interval.begin;
const auto& t_end = _time_interval.end;
using namespace Dune::FloatCmp;
if (lt(time, t_begin) or gt(time, t_end))
{
auto log = Dorie::get_logger(Dorie::log_base);
log->error("Querying boundary condition '{}' outside its time "
"interval [{}, {}] with time: {}",
_name, t_begin, t_end, time);
return false;
}
return true;
}
public:
/// Construct a static boundary condition
BoundaryCondition (const std::string_view name,
const TimeInterval<Time> interval,
const RF value)
:
_time_interval(interval),
_value_begin(value),
_name(name)
{ }
/// Construct an interpolated boundary condition
BoundaryCondition (const std::string_view name,
const TimeInterval<Time> interval,
const RF value_begin,
const std::optional<RF> value_end)
:
_time_interval(interval),
_value_begin(value_begin),
_value_end(value_end),
_name(name)
{ }
/// Destroy this object
virtual ~BoundaryCondition () = default;
/// Return the name of this boundary condition
std::string_view name () const { return _name; }
/// Return the time interval of this boundary condition
const TimeInterval<Time>& time_interval () const
{
return _time_interval;
}
/// Return the type of this boundary condition
virtual Operator::BCType type () const
{
return Operator::BCType::none;
}
/// Return the value for a given point in time
/** \param time The time at which to evaluate the BC
* \return The (possibly interpolated) value
* \throw InvalidStateException if the time is not inside the interval
* \todo Maybe use `assert` on in_interval
*/
virtual RF evaluate (const Time time) const
{
if (not in_interval(time)) {
DUNE_THROW(InvalidStateException, "Evaluating boundary condition "
"with invalid time");
}
// without interpolation, return the first value
if (not _value_end)
return _value_begin;
// else: interpolate values
const auto& t_begin = _time_interval.begin;
const auto& t_end = _time_interval.end;
const auto& v_end = _value_end.value();
return _value_begin
+ (time - t_begin)
* (v_end - _value_begin)
/ (t_end - t_begin);
}
};
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_COMMON_BOUNDARY_CONDITION_BASE_HH
#ifndef DUNE_DORIE_COMMON_BOUNDARY_CONDITION_DIRICHLET_HH
#define DUNE_DORIE_COMMON_BOUNDARY_CONDITION_DIRICHLET_HH
#include <dune/dorie/common/boundary_condition/boundary_condition.hh>
namespace Dune {
namespace Dorie {
/// Toggle for concentration represented by Dirichlet BC
/** Used for \ref TransportModel model only.
*
* The relation of water_phase and total is given by
* \f[
* C_w = C_t \theta_w \, ,
* \f]
* where \f$ \theta_w \f$ is the volumetric water content.
*/
enum class SoluteConcentration {
total, //!< Total concentration in the volume \f$ C_t \f$.
water_phase, //!< Concentration in water phase \f$ C_w \f$.
none //!< No type, used for \ref RichardsModel model.
};
/// Class representing a Dirichlet boundary condition (fixed solution)
/** In the \ref TransportModel model, the #_solute_concentration switch is used
* to determine which type of solute concentration the boundary condition
* refers to.
*
* \ingroup Boundary
*/
template<typename RF>
class DirichletBoundaryCondition : public BoundaryCondition<RF>
{
protected:
using Base = BoundaryCondition<RF>;
using Time = typename Base::Time;
/// Setting for the type of concentration encoded in this BC value
const SoluteConcentration _solute_concentration;
public:
/// Construct a static boundary condition
DirichletBoundaryCondition (
const std::string_view name,
const TimeInterval<Time> interval,
const RF value,
const SoluteConcentration concentration=SoluteConcentration::none)
:
Base(name, interval, value),
_solute_concentration(concentration)
{ }
/// Construct an interpolated boundary condition
DirichletBoundaryCondition (
const std::string_view name,
const TimeInterval<Time> interval,
const RF value_begin,
const std::optional<RF> value_end,
const SoluteConcentration concentration=SoluteConcentration::none)
:
Base(name, interval, value_begin, value_end),
_solute_concentration(concentration)
{ }
/// Destroy this object
virtual ~DirichletBoundaryCondition () override = default;
/// Return the type of this boundary condition
virtual Operator::BCType type () const override
{
return Operator::BCType::Dirichlet;
}
/// The type string of this boundary condition class
static constexpr std::string_view type_str = "dirichlet";
/// Return the quantity encoded in this Dirichlet value
SoluteConcentration concentration_type () const
{
return _solute_concentration;
}
};
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_COMMON_BOUNDARY_CONDITION_DIRICHLET_HH
This diff is collapsed.
This diff is collapsed.
// These groups should gather information on strongly related portions of
// the code. In them we aim to explain to other developers:
// - The rationale of the design
// - How careful someone has to be when modifying it
// - Specific naming or other conventions upheld in the code
/**
@defgroup Boundary Boundary Conditions
@{
@ingroup Common
@brief Representing boundary conditions for all models
## Overview
There are several types of boundary conditions. They all store the time
interval during which they can be evaluated. Typically, one static value,
or two values for interpolation, are stored. Since boundary conditions are
implemented as derived classed of Dune::Dorie::BoundaryCondition, they
may also store additional members and define additional methods (see
Dune::Dorie::NeumannBoundaryCondition for an example).
## Implementation
Boundary conditions are handled by the
Dune::Dorie::BoundaryConditionManager. It loads the user input data and
assembles the necessary data structures. Its behavior is the same for all
models, with the difference that models may not support all types of
boundary conditions. Therefore, the manager has a template parameter for
setting its model mode, which is propagated to the
Dune::Dorie::BoundaryConditionFactory. New boundary conditions have to be
"registered" there, making sure they can be used for a certain (or all)
models.
@}
**/
#ifndef DUNE_DORIE_COMMON_BOUNDARY_CONDITION_NEUMANN_HH
#define DUNE_DORIE_COMMON_BOUNDARY_CONDITION_NEUMANN_HH
#include <dune/dorie/common/boundary_condition/boundary_condition.hh>
namespace Dune {
namespace Dorie {
/// Class representing a Neumann boundary condition (fixed flux)
/** In the \ref RichardsModel model, this class typically models irrigation and
* evaporation from soil surfaces and therefore stores an additional switch
* which can be used to consider the surface tilt when evaluating the boundary
* condition.
*
* If #_horizontal_projection is set to true, the local operator scales the
* boundary condition \f$ j \f$ with
* \f[
* j = j_N \left| \mathbf{\hat{n}}_F \cdot \mathbf{\hat{g}}_F \right|
* \f]
* where the vectors are the unit normal vector of the boundary face, and
* the unit vector along the gravitational force, respectively.
*
* #_horizontal_projection is ignored in the \ref TransportModel model.
*
* \ingroup Boundary
*/
template<typename RF>
class NeumannBoundaryCondition : public BoundaryCondition<RF>
{
protected:
using Base = BoundaryCondition<RF>;
using Time = typename Base::Time;
/// Value of the surface adjustment switch
const bool _horizontal_projection;
public:
/// Construct a static boundary condition
NeumannBoundaryCondition (const std::string_view name,
const TimeInterval<Time> interval,
const RF value,
const bool horizontal_projection)
:
Base(name, interval, value),
_horizontal_projection(horizontal_projection)
{ }
/// Construct an interpolated boundary condition
NeumannBoundaryCondition (const std::string_view name,
const TimeInterval<Time> interval,
const RF value_begin,
const std::optional<RF> value_end,
const bool horizontal_projection)
:
Base(name, interval, value_begin, value_end),
_horizontal_projection(horizontal_projection)
{ }
/// Destroy this object
virtual ~NeumannBoundaryCondition () override = default;
/// Return the type of this boundary condition
virtual Operator::BCType type () const override
{
return Operator::BCType::Neumann;
}
/// Return true the BC is to be projected horizontally
bool horizontal_projection () const { return _horizontal_projection; }
/// The type string of this boundary condition class
static constexpr std::string_view type_str = "neumann";
};
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_COMMON_BOUNDARY_CONDITION_NEUMANN_HH
#ifndef DUNE_DORIE_COMMON_BOUNDARY_CONDITION_OUTFLOW_HH
#define DUNE_DORIE_COMMON_BOUNDARY_CONDITION_OUTFLOW_HH
#include <dune/dorie/common/boundary_condition/boundary_condition.hh>
namespace Dune {
namespace Dorie {
/// Class representing an Outflow BC (matter may leave through the boundary)
/** The exact implementation depends on the model.
*
* \ingroup Boundary
*/
template<typename RF>
class OutflowBoundaryCondition : public BoundaryCondition<RF>
{
protected:
using Base = BoundaryCondition<RF>;
using Time = typename Base::Time;
public:
/// Re-use base constructor
using Base::Base;
/// Destroy this object
virtual ~OutflowBoundaryCondition () override = default;
/// Return the type of this boundary condition
virtual Operator::BCType type () const override
{
return Operator::BCType::Outflow;
}
/// The type string of this boundary condition class
static constexpr std::string_view type_str = "outflow";
};
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_COMMON_BOUNDARY_CONDITION_OUTFLOW_HH
#ifndef DUNE_DORIE_UTIL_HH
#define DUNE_DORIE_UTIL_HH
#include <iostream>
#include <memory>
#include <vector>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#ifndef DUNE_DORIE_COMMON_GFS_HELPER_HH
#define DUNE_DORIE_COMMON_GFS_HELPER_HH
#include <dune/geometry/type.hh>
......@@ -19,40 +12,6 @@
namespace Dune{
namespace Dorie{
/**
* @brief Enum class for output policy.
* @details It defines when a model class should produce an output.
*/
enum class OutputPolicy {None,EndOfRichardsStep,EndOfTransportStep};
/**
* @brief Enum class for output policy.
* @details It defines which variable is the target in order to mark the grid
*/
enum class AdaptivityPolicy {None,WaterFlux,SoluteFlux};
/// Return the estimation of entries per matrix row for the spatial GridOperator
/** This supposedly decreases matrix assembly time.
* The values specify the *blocks* per row. DG assembles one block for the
* actual element and one for each of its neighbors.
* \param dim Spatial dimension
* \param geo Geometry type of grid entities
* \return Estimated number of blocks per matrix row
*/
template<typename R = std::size_t>
constexpr R estimate_mbe_entries (const int dim, const Dune::GeometryType::BasicType geo)
{
if (geo==Dune::GeometryType::BasicType::cube){
return 2*dim + 1;
}
else if (geo==Dune::GeometryType::BasicType::simplex){
return dim + 2;
}
throw std::runtime_error("Cannot provide matrix backend entry estimate "
"for given geometry type!");
}
/// Provide types and construction of the GFS for the linear solver
template<typename GridView, typename RF, Dune::GeometryType::BasicType GeometryType>
struct LinearSolverGridFunctionSpaceHelper
......@@ -195,38 +154,7 @@ public:
}
};
/// Traits struct defining basic types based on grid an geometry
template<class GridType, Dune::GeometryType::BasicType GeometryType>
struct BaseTraits
{
static constexpr int dim = GridType::dimension;
static constexpr Dune::GeometryType::BasicType GridGeometryType = GeometryType;
using TF = double;
using TimeField = TF;
using RF = double;
using RangeField = RF;
using Array = std::vector<RangeField>;
using Scalar = Dune::FieldVector<RF,1>;
using Vector = Dune::FieldVector<RF,dim>;
using Tensor = Dune::FieldMatrix<RF,dim,dim>;
using Index = unsigned int;
using IndexArray = Dune::FieldVector<Index,dim>;
using Grid = GridType;
using DomainField = typename Grid::ctype;
using DF = DomainField;
using Domain = Dune::FieldVector<DF,dim>;
using IntersectionDomain = Dune::FieldVector<DF,dim-1>;
using GV = typename Grid::LeafGridView;
using GridView = GV;
using Element = typename GV::Traits::template Codim<0>::Entity;
using Intersection = typename GV::Traits::Intersection;
};
}
}
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_UTIL_HH
#endif // DUNE_DORIE_COMMON_GFS_HELPER_HH
#ifndef DUNE_DORIE_TYPEDEFS_HH
#define DUNE_DORIE_TYPEDEFS_HH
#include <dune/pdelab/common/function.hh>
namespace Dune{
namespace Dorie{
//! Boundary condition types
struct BoundaryCondition
{
//! Enum for the Boundary condition type
enum Type {
Neumann, //!< Fixed flux at boundary
Dirichlet, //!< Fixed matric head at boundary
Outflow, //!< Fixed outflow boundary condition
Other //!< UNUSED
};
//! Test for Dirichlet boundary condition
static bool isDirichlet (Type i)
{
return (i == Dirichlet);
}
//! Test for Neumann boundary condition
static bool isNeumann (Type i)
{
return (i == Neumann);
}
//! Test for Outflow boundary condition
static bool isOutflow (Type i)
{
return (i == Outflow);
}
//! UNUSED: Test for other(?) boundary condition
static bool isOther (Type i)
{
return (i == Other);
}
};
/**
* @brief Homogeneous initial condition: Matric head h=0
*/
template<typename Traits>
class ZeroInitial
: public Dune::PDELab::AnalyticGridFunctionBase<Dune::PDELab::AnalyticGridFunctionTraits<typename Traits::GridTraits::GridView, typename Traits::GridTraits::RangeField,1>,
ZeroInitial<Traits> >
{
public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<typename Traits::GridTraits::GridView, typename Traits::GridTraits::RangeField,1> AGFTraits;
typedef Dune::PDELab::AnalyticGridFunctionBase<AGFTraits,ZeroInitial<Traits> > BaseT;
enum {dim = Traits::GridTraits::dim};
ZeroInitial (const Dune::ParameterTree& config, const typename Traits::GridTraits::GridView gv) : BaseT(gv) {}
void evaluateGlobal (const typename Traits::GridTraits::Domain& x, typename Traits::GridTraits::Scalar& y) const
{
y = 0.;
}
};
}
}
#include <vector>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/geometry/type.hh>
namespace Dune {
namespace Dorie {
namespace Operator {
//! BC types implemented in the local operator(s)
enum class BCType
{
Neumann, //!< Fixed flux at boundary
Dirichlet, //!< Fixed matric head at boundary
Outflow, //!< Free flow of solute through boundary
none //!< No BC specified
};
/// Modes for the Dirichlet BC of the solute transport local operators
struct DirichletMode
{
enum Type { SoluteConcentration, //!< Concentration in water phase
TotalSolute, //!< Concentration in total volume
};
};
} // namespace Operator
/**
* @brief Enum class for output policy. It defines when a simulation class
* should produce an output
*/
enum class OutputPolicy {None, EndOfTransportStep, EndOfRichardsStep};
/**
* @brief Enum class for output policy. It defines which variable
* is the target in order to mark the grid
*/
enum class AdaptivityPolicy {None, WaterFlux, SoluteFlux};
/// A wrapper for begin and end time stamps of a time interval
/** \tparam TF Arithmetic data type for time stamps
*/
template<typename TF>
struct TimeInterval {
TF begin;
TF end;
};
/// Traits struct defining basic types based on grid an geometry
template<class GridType, Dune::GeometryType::BasicType GeometryType>
struct BaseTraits
{
static constexpr int dim = GridType::dimension;
static constexpr Dune::GeometryType::BasicType GridGeometryType = GeometryType;
using TF = double;
using TimeField = TF;
using RF = double;
using RangeField = RF;
using Array = std::vector<RangeField>;
using Scalar = Dune::FieldVector<RF,1>;
using Vector = Dune::FieldVector<RF,dim>;
using Tensor = Dune::FieldMatrix<RF,dim,dim>;
using Index = unsigned int;
using IndexArray = Dune::FieldVector<Index,dim>;
using Grid = GridType;
using DomainField = typename Grid::ctype;
using DF = DomainField;
using Domain = Dune::FieldVector<DF,dim>;
using IntersectionDomain = Dune::FieldVector<DF,dim-1>;
using GV = typename Grid::LeafGridView;
using GridView = GV;
using Element = typename GV::Traits::template Codim<0>::Entity;
using Intersection = typename GV::Traits::Intersection;
//using GFS = GridFunctionSpaceHelper<GV,RF,order,GeometryType>;
//using AdaptivityGFS = GridFunctionSpaceHelper<GV,RF,0,GeometryType>;
//using LinearSolverGFS = GridFunctionSpaceHelper<GV,RF,0,GeometryType>;
};
} // namespace Dorie
} // namespace Dune
#endif
......@@ -6,6 +6,7 @@
#include <algorithm>
#include <cctype>
#include <memory>
#include <exception>
#include <utility>
#include <sys/types.h>
......@@ -18,11 +19,36 @@
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/geometry/type.hh>
#include <dune/grid/common/gridview.hh>
namespace Dune {
namespace Dorie {
/// Return the estimation of entries per matrix row for the spatial GridOperator
/** This supposedly decreases matrix assembly time.
* The values specify the *blocks* per row. DG assembles one block for the
* actual element and one for each of its neighbors.
* \param dim Spatial dimension
* \param geo Geometry type of grid entities
* \return Estimated number of blocks per matrix row
*/
template<typename R = std::size_t>
constexpr R estimate_mbe_entries (const int dim,
const Dune::GeometryType::BasicType geo)
{
if (geo==Dune::GeometryType::BasicType::cube){
return 2*dim + 1;
}
else if (geo==Dune::GeometryType::BasicType::simplex){
return dim + 2;
}
throw std::runtime_error("Cannot provide matrix backend entry estimate "
"for given geometry type!");
}
/// Retrieve the grid extensions as a position vector
/// Retrieve a pair of vectors spanning the grid domain
/** Iterate over all vertices of the coarsest level grid view, determine the
* locally largest and smallest extensions per dimension, and compute the
......
......@@ -2,11 +2,12 @@
#define DUNE_DORIE_MODEL_BASE_HH
#include <dune/dorie/common/logging.hh>
#include <dune/dorie/common/typedefs.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/float_cmp.hh>
#include <dune/dorie/common/util.hh>
#include <dune/dorie/common/typedefs.hh>
namespace Dune{
namespace Dorie{
......@@ -247,27 +248,27 @@ protected:
### Adaptivity:
For the purpose of having a generic `run()` method in this class, the sketch
of the adaptivity must be provided here. First, adaptivity must be never
done inside the method `step()`. Additionally, if the model can perform
adaptivity, it must provide the method `mark_grid()`, otherwise, it must
throw a `NotImplemented` exception. For completion, there must be a method
called `adapt_grid()` that pack all the known degrees of freedom and, the
grid function spaces and perform the adaptivity.
Similar to the idea of `OutputPolicy` there must be an `AdaptivityPolicy`
set by the methods `set_policy(AdaptivityPolicy)` and `adaptivity_policy()`
respectively which will be set by default to `AdaptivityPolicy::None`.
Having this separation between marking the grid and adaptivity of solutions
leave coupled models to decide which model will mark the grid inside its
method of `mark_grid()` and still handle the adaptivity implementing
For the purpose of having a generic `run()` method in this class, the sketch
of the adaptivity must be provided here. First, adaptivity must be never
done inside the method `step()`. Additionally, if the model can perform
adaptivity, it must provide the method `mark_grid()`, otherwise, it must
throw a `NotImplemented` exception. For completion, there must be a method
called `adapt_grid()` that pack all the known degrees of freedom and, the
grid function spaces and perform the adaptivity.
Similar to the idea of `OutputPolicy` there must be an `AdaptivityPolicy`
set by the methods `set_policy(AdaptivityPolicy)` and `adaptivity_policy()`
respectively which will be set by default to `AdaptivityPolicy::None`.
Having this separation between marking the grid and adaptivity of solutions
leave coupled models to decide which model will mark the grid inside its
method of `mark_grid()` and still handle the adaptivity implementing
`adapt_grid()` separately.
### Time control:
A minimal interface is only providing `begin_time()`, `end_time()` and
`current_time()`.
A minimal interface is only providing `begin_time()`, `end_time()` and
`current_time()`.
In order to keep solution vectors minimal at the time of adaptivity, DORiE
expect the several models synchronized in time. Therefore, a good and still
......@@ -283,12 +284,12 @@ protected:
Since adaptivity can modify the complete setup of the model, and since
the method `adapt_grid()` is not necessarily called by the model.
There is a need for the method `post_adapt_grid()` which is in charge of
setup everything again to the new grid. Same argument holds for
setup everything again to the new grid. Same argument holds for
`pre_adapt_grid()`.
### Run:
With all the definitions given before, a general algorithm to run the
With all the definitions given before, a general algorithm to run the
program can be given by following code:
```
......@@ -300,7 +301,7 @@ protected:
while( current_time() < end_time() )
{
step();
if (adaptivity_policy() != AdaptivityPolicy::None)
if ( current_time() < end_time() )
{
......
......@@ -6,6 +6,10 @@
#include <dune/dorie/model/richards/richards.hh>
#include <dune/dorie/model/transport/transport.hh>
#ifdef GTEST
#include <gtest/gtest.h>
#endif
namespace Dune{
namespace Dorie{
......@@ -207,6 +211,12 @@ public:
if (_transport->output_policy() != OutputPolicy::None)
_transport->write_data();
}
private:
#ifdef GTEST
FRIEND_TEST(ModelCoupling, TimeSteps);
#endif
};
......
......@@ -24,7 +24,7 @@
#include <dune/dorie/common/logging.hh>
#include <dune/dorie/common/utility.hh>
#include <dune/dorie/common/grid_creator.hh>
#include <dune/dorie/common/util.hh>
#include <dune/dorie/common/utility.hh>
#include <dune/dorie/model/base.hh>
......
......@@ -20,7 +20,7 @@
namespace Dune {
namespace Dorie {
namespace Operator {
namespace Operator {
/// Local operator for residual-based error estimation.
/*
......@@ -208,6 +208,14 @@ namespace Dune {
const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
R& r_s) const
{
// query the boundary condition
const auto bc = boundary.bc(ig.intersection());
// check boundary condition type
if (bc->type() != BCType::Dirichlet){
return;
}
// get polynomial degree
const int degree = lfsu_s.finiteElement().localBasis().order();
const int intorder = intorderadd + quadrature_factor * degree;
......@@ -229,12 +237,6 @@ namespace Dune {
// loop over quadrature points and integrate normal flux
for (const auto& it : quadratureRule(gtface,intorder))
{
// check boundary condition type
const auto bc = boundary.bc(ig.intersection(), it.position(), time);
if (!BoundaryCondition::isDirichlet(bc)){
continue;
}
// position of quadrature point in local coordinates of elements
const Domain local_s = ig.geometryInInside().global(it.position());
......@@ -248,7 +250,7 @@ namespace Dune {
u_s += x_s(lfsu_s,i) * phi_s[i];
// boundary condition at position
const RF g = boundary.g(ig.intersection(), it.position(), time);
const RF g = bc->evaluate(time);
// compute jump in solution
const RF h_jump = u_s - g;
......
// -*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_DORIE_RICHARDS_BOUNDARY_HH
#define DUNE_DORIE_RICHARDS_BOUNDARY_HH
#include <memory>
#include <dune/common/exceptions.hh>
#include <dune/common/parametertree.hh>
#include <dune/dorie/common/boundary_condition.hh>
namespace Dune{
namespace Dorie{
/// Boundary type and condition value queries
/** This class containts functions that return the type of boundary conditions
* and the values of the boundary condition parameters. Both types of queries
* can be time dependent.
*/
template<typename Traits>
class FlowBoundary
{
private:
typedef typename Traits::RangeField RF;
typedef typename Traits::Domain Domain;
typedef typename Traits::IntersectionDomain ID;
typedef typename Traits::Intersection Intersection;
enum {dim = Traits::dim};
//! Object for handling the BC data file and function queries
std::unique_ptr<BCReadoutInterface<Traits>> bcDataHandler;
public:
/// Read BC file type and create data handling object
/** \see Dune::Dorie::BCReadoutInterface
* \throw Dune::IOError BC File Type not supported
*/
FlowBoundary(const Dune::ParameterTree& config)
{
const auto log = Dorie::get_logger(log_richards);
const auto filename = config.get<std::string>("boundary.file");
std::string bcFileType = config.get<std::string>("boundary.fileType");
log->info("Reading boundary condition data file: {}", filename);
if(bcFileType == "rectangularGrid") {
log->debug(" Using 'rectangularGrid' BC interface");
bcDataHandler = std::make_unique<RectangularGrid<Traits>>(config);
}
else {
log->error("Unknown boundary condition file type: {}", bcFileType);
DUNE_THROW(IOError,"unknown bcFileType specified!");
}
}
/// Return next time at which any boundary condition changes. Return -1.0 by default
/** \param time Current time
*/
RF getNextTimeStamp (const RF& time) const {
if(bcDataHandler)
return bcDataHandler->getNextTimeStamp(time);
return -1.0;
}
/// Return Boundary Condition Type at certain position and time
/** \param is Intersection enitity
* \param x Position in local entity coordinates
* \param time Time value
* \param xGlobal Global coordinate
* \return Boundary condition type enumerator
*/
BoundaryCondition::Type bc (const Intersection& is, const ID& x, const RF& time) const
{
const Domain xGlobal = is.geometry().global(x);
if(is.boundary()) {
return bcDataHandler->getBCtype(xGlobal, time);
}
return BoundaryCondition::Other; // unknown
}
/**
* @brief Dirichlet boundary condition at certain position and time
* \param is Intersection enitity
* \param x Position in local entity coordinates
* \param time Time value
* \param xglobal Global coordinate
* \return Value of matric head
*/
RF g (const Intersection& is, const ID& x, const RF& time) const
{
const Domain xGlobal = is.geometry().global(x);
if(is.boundary()) {
if(bcDataHandler)
return bcDataHandler->getDirichletValue(xGlobal, time);
}
return 0.0;
}
/**
* @brief Neumann boundary condition at certain position and time
* \param is Intersection enitity
* \param x Position in local entity coordinates
* \param time Time value
* \param xglobal Global coordinate
* \return Value of flux
*/
RF j (const Intersection& is, const ID& x, const RF& time) const
{
const Domain xGlobal = is.geometry().global(x);
if(is.boundary()) {
if(bcDataHandler)