...
 
Commits (25)
......@@ -269,7 +269,7 @@ adding an empty line, make text **bold** or ``monospaced``.
</parameter>
<parameter name="checkTol">
<definition> Whenever ``checkJumps`` is activated, it check that flux
<definition> Whenever ``checkJumps`` is activated, it checks that flux
reconstruction engine is creating conforming normal fluxes up to ``jumpTol``.
</definition>
<values> float &gt; 0 </values>
......@@ -277,6 +277,24 @@ adding an empty line, make text **bold** or ``monospaced``.
</parameter>
</category>
<category name="conservationOfMass">
<parameter name="check">
<definition> Check that mass is conserved in the domain with respect to
boundary fluxes.
</definition>
<suggestion> none </suggestion>
<values> none, warn, error </values>
</parameter>
<parameter name="checkTol">
<definition> Whenever ``checkJumps`` is activated, it check that mass is
conserved in the domain with respect to boundary fluxes. ``checkTol``
is the persentage of mass that can be lost without a warning/error.
</definition>
<values> float &gt; 0 </values>
<suggestion> 1E-10 </suggestion>
</parameter>
</category>
<category name="NewtonParameters">
<parameter name="ForceIteration" hidden="true">
......
......@@ -249,6 +249,13 @@ adding an empty line, make text **bold** or ``monospaced``.
</category>
<category name="numerics">
<parameter name="penaltyFactor">
<definition> Penalty factor to be used in the Discontinuous Galerkin scheme
</definition>
<values> float </values>
<suggestion> 10 </suggestion>
</parameter>
<parameter name="timestepMethod">
<definition> Numerical scheme to perform time steps in the simulation.
``alex2`` and ``implicit_euler`` are implicit methods.
......@@ -297,7 +304,7 @@ adding an empty line, make text **bold** or ``monospaced``.
</parameter>
<parameter name="checkTol">
<definition> Whenever ``checkJumps`` is activated, it check that flux
<definition> Whenever ``checkJumps`` is activated, it checks that flux
reconstruction engine is creating conforming normal fluxes up to ``jumpTol``.
</definition>
<values> float &gt; 0 </values>
......@@ -305,6 +312,25 @@ adding an empty line, make text **bold** or ``monospaced``.
</parameter>
</category>
<category name="conservationOfMass">
<parameter name="check">
<definition> Check that mass is conserved in the domain with respect to
boundary fluxes.
</definition>
<suggestion> none </suggestion>
<values> none, warn, error </values>
</parameter>
<parameter name="checkTol">
<definition> Whenever ``checkJumps`` is activated, it check that mass is
conserved in the domain with respect to boundary fluxes. ``checkTol``
is the persentage of mass that can be lost without a warning/error.
</definition>
<values> float &gt; 0 </values>
<suggestion> 1E-10 </suggestion>
</parameter>
</category>
<category name="solverParameters">
<parameter name="reduction">
<definition> Required relative defect reduction.
......@@ -322,4 +348,41 @@ adding an empty line, make text **bold** or ``monospaced``.
<values> float </values>
</parameter>
</category>
<category name="numerics.experimental">
<parameter name="method">
<definition> DG discretization method for skeleton terms.
**SIPG:** Symmetric Interior Penalty
**NIPG:** Non-Symmetric Interior Penalty
**OOB:** Oden, Babuska, Baumann: no penalty term
**IIP:** Incomplete Interior Penalty: no symmetry term
</definition>
<values> SIPG, NIPG, OOB, IIP </values>
<suggestion> SIPG </suggestion>
<comment> Experimental settings are enabled by the appropriate CMake flag.
</comment>
</parameter>
<parameter name="upwinding">
<definition> Upwinding method for skeleton terms.
**semiUpwind:** Apply upwinding to conductivity factor (only).
**fullUpwind:** Apply upwinding on numeric flux and conductivity.
</definition>
<values> none, semiUpwind, fullUpwind </values>
<suggestion> none </suggestion>
</parameter>
<parameter name="weights">
<definition> Apply harmonic weighting to skeleton term contributions.
</definition>
<values> true, false </values>
<suggestion> true </suggestion>
</parameter>
</category>
</dorie>
......@@ -283,9 +283,6 @@ public:
}
private:
/// Controls how to deal with jumps occurring in the reconstructed flux
enum class CheckLevel { none, warn, error };
std::shared_ptr<spdlog::logger> _log;
GV _gv;
FEMU _femu;
......@@ -297,25 +294,6 @@ private:
/// The jump check level
const CheckLevel _check_lvl;
const RF _check_tol;
/// Return a check level for a given string
CheckLevel to_check_lvl (const std::string check_lvl_str)
{
if (to_lower(check_lvl_str) == "none") {
return CheckLevel::none;
}
else if (check_lvl_str == "warn") {
return CheckLevel::warn;
}
else if (check_lvl_str == "error") {
return CheckLevel::error;
}
else {
_log->error("Unknown check for jumps in reconstructed flux: {}",
check_lvl_str);
DUNE_THROW(IOError, "Unknown check for flux reconstruction jumps");
}
}
};
} // namespace Dorie
......
#ifndef DUNE_DORIE_GRID_FUNCTION_INTEGRATOR_HH
#define DUNE_DORIE_GRID_FUNCTION_INTEGRATOR_HH
namespace Dune {
namespace Dorie {
/*-------------------------------------------------------------------------*//**
* @brief Integrate grid function
*
* @author Santiago Ospina De Los Ríos
* @date 2019
*
* @param[in] gf The grid function
* @param[in] int_order The integration order
*
* @tparam GF The grid function type
*
* @return The integral over the whole domain
*/
template<class GF>
typename GF::Traits::RangeType
gf_integrator(const GF& gf, unsigned int int_order)
{
using Traits = typename GF::Traits;
using Range = typename Traits::RangeType;
Range integral(0.);
const auto& gv = gf.getGridView();
// Loop over the grid
for (const auto& cell : Dune::elements(gv,Dune::Partitions::interior))
{
// Get geometry
auto geo = cell.geometry();
// Integrate on the cell
for (const auto& point : quadratureRule(geo, int_order))
{
// Evaluate function
Range y;
const auto& x = point.position();
gf.evaluate(cell,x,y);
auto factor = point.weight() * geo.integrationElement(x);
// Accumulate
integral += y*factor;
}
}
// Accumulate over processors
return gv.comm().sum(integral);
}
/*-------------------------------------------------------------------------*//**
* @brief Integrate grid function on the boundary of the domain
* @details If the grid function is a vector valued function, it will be
* integrated in the normal direction of the boundary.
*
* @author Santiago Ospina De Los Ríos
* @date 2019
*
* @param[in] gf The grid function
* @param[in] int_order The integration order
*
* @tparam GF The grid function type
*
* @return The integral over the whole domain boundary
*/
template<class GF>
Dune::FieldVector<typename GF::Traits::RangeFieldType,1>
boundary_gf_integrator(const GF& gf, unsigned int int_order)
{
using Traits = typename GF::Traits;
using Range = typename Traits::RangeType;
constexpr int dim_range = Traits::dimRange;
constexpr int dim_domain = Traits::dimDomain;
static_assert(dim_range == 1 || dim_range == dim_domain);
Dune::FieldVector<typename GF::Traits::RangeFieldType,1> integral(0.);
const auto& gv = gf.getGridView();
// Loop over the grid
for (const auto& cell : Dune::elements(gv,Dune::Partitions::interior))
{
// Loop over the intersections
for (const auto& face : Dune::intersections(gv,cell))
{
// Only integrate on boundary
if (face.boundary() and not face.neighbor())
{
// Get geometry
auto geo_f = face.geometry();
auto geo_in_i = face.geometryInInside();
// Integrate on the face
for (const auto& point : quadratureRule(geo_f, int_order))
{
// Relative positions
const auto& position_f = point.position();
const auto position_i = geo_in_i.global(position_f);
// Face normal vector
const auto normal = face.unitOuterNormal(position_f);
// Evaluate function
Range y;
gf.evaluate(cell,position_i,y);
const auto factor = point.weight() * geo_f.integrationElement(position_f);
// Accumulate
if constexpr (dim_range == dim_domain)
integral += y * normal * factor;
else
integral += y * factor;
}
}
}
}
// Accumulate over processors
return gv.comm().sum(integral);
}
}
}
#endif
// -*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_DORIE_LOCAL_OPERATOR_DG_HH
#define DUNE_DORIE_LOCAL_OPERATOR_DG_HH
/**
@defgroup LocalOperators Local operators
@{
Local operators are in many senses the heart of dorie; in them, we transform
finite element basis into a vector of residuals and a mass matrix taking into
account the specific equation and method we want to solve. These
transformations follow the structure of dune-pdelab. In particular,
some care has to be taken regarding the different frames of reference when
working with coordinates and gradients. It can get particularly messy
when working with intersections because one has to express the same
position/gradient with respect to the inside, outside, and face entities,
and sometimes with respect to the macro-grid (e.g. global coordinates).
Therefore, in dorie we try to use the following suffixes convention :
| Convention | Meaning |
|------------|---------|
| `<description>` + `_i` | `<description>` with respect to the **inside** entity. |
| `<description>` + `_o` | `<description>` with respect to the **outside** entity. |
| `<description>` + `_f` | `<description>` with respect to the **face** entity (or intersection). |
| `<description>` + `_g` | `<description>` with respect to the **global** grid. |
| `<description>` + `_n` | `<description>` in **normal** direction. |
Notice that when working with gradients, the suffix for normal direction can
be combined with the first three suffixes (`_n_i`, `_n_o`, and `_n_f`). When
it's not specified, is understood that normal directions are referenced with
respect to the face (`_n` \f$\equiv\f$. `_n_f`).
@todo Update Dune::Dorie::Operator::RichardsDGSpatialOperator to the
local operator convention.
@}
**/
namespace Dune {
namespace Dorie {
namespace Operator {
/**
* @brief Interior penalty type for dg methods
* @ingroup LocalOperators
*/
struct DGMethod
{
enum Type { NIPG, //!< Non-symmetric interior penalty
SIPG, //!< Symmetric interior penalty (default)
OOB, //!< Oden, Babuska, Baumann (no penalty term)
IIP //!< Incomplete interior penalty (no symmetry)
};
};
/**
* @brief Upwinding type
* @ingroup LocalOperators
*/
struct DGUpwinding
{
enum Type {
FullUpwind, //!< Upwind flux and conductiviy
SemiUpwind, //!< Only upwind conductivity factor
None //!< No upwinding
};
};
/**
* @brief Gradient weighting
* @ingroup LocalOperators
*/
struct DGWeights
{
enum Type {
WeightsOn, //!< harmonic weighting of flux contributions
WeightsOff //!< no weighting of flux contributions
};
};
/**
* @brief Operator internal BC handling
* @ingroup LocalOperators
*/
struct BCType
{
enum Type { Neumann, //!< Fixed flux at boundary
Dirichlet, //!< Fixed matric head at boundary
Outflow, //!<
None //!< No BC specified
};
};
/**
* @brief Read experimental operator settings from a parameter file
* @ingroup LocalOperators
*
* @param[in] inifile The inifile
*
* @return Tuple of DGMethod, DGUpwinding, and DGWeights to be inserted into
* the local operator constructor.
*/
inline auto read_experimental_operator_settings(
const Dune::ParameterTree &inifile,
std::shared_ptr<spdlog::logger> log)
-> std::tuple<Operator::DGMethod::Type,
Operator::DGUpwinding::Type,
Operator::DGWeights::Type>
{
log->debug("Reading experimental DG operator settings:");
auto& dg_ini = inifile.sub("numerics.experimental");
DGMethod::Type method;
DGUpwinding::Type upwinding;
DGWeights::Type weights;
const auto method_str = dg_ini.get<std::string>("method");
log->debug(" DG method: {}", method_str);
if (method_str == "SIPG")
method = DGMethod::SIPG;
else if (method_str == "NIPG")
method = DGMethod::NIPG;
else if (method_str == "OOB")
method = DGMethod::OOB;
else if (method_str == "IIP")
method = DGMethod::IIP;
else{
log->error("Unknown DG method: {}", method_str);
DUNE_THROW(Dune::IOError, "Unknown DG method");
}
const auto upwinding_str
= dg_ini.get<std::string>("upwinding");
log->debug(" Upwinding scheme: {}", upwinding_str);
if (upwinding_str == "none")
upwinding = DGUpwinding::None;
else if (upwinding_str == "semiUpwind")
upwinding = DGUpwinding::SemiUpwind;
else if (upwinding_str == "fullUpwind")
upwinding = DGUpwinding::FullUpwind;
else {
log->error("Unknown upwinding scheme: {}", upwinding_str);
DUNE_THROW(Dune::IOError, "Unknown upwinding scheme");
}
const auto weights_str = dg_ini.get<bool>("weights");
log->debug(" Harmonic weights at interfaces: {}", weights_str);
if (weights_str)
weights = DGWeights::WeightsOn;
else
weights = DGWeights::WeightsOff;
return std::make_tuple(method, upwinding, weights);
}
} // namespace Operator
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_LOCAL_OPERATOR_DG_HH
\ No newline at end of file
#ifndef DUNE_DORIE_ONE_STEP_METHOD_LINEAR_HH
#define DUNE_DORIE_ONE_STEP_METHOD_LINEAR_HH
#include <dune/pdelab/stationary/linearproblem.hh>
#include <dune/pdelab/gridoperator/onestep.hh>
#include <dune/pdelab/instationary/onestep.hh>
namespace Dune{
namespace Dorie{
template<class T, class IGOS, class LS, class TrlV, class TstV = TrlV, class TC = Dune::PDELab::SimpleTimeController<T> >
class OneStepMethodLinear;
// explicit case
template<class T, class GO0, class GO1, class LS, class TrlV, class TstV, class TC>
class OneStepMethodLinear<
T,
Dune::PDELab::OneStepGridOperator<GO0,GO1,false>,
LS,
TrlV,
TstV,
TC>
: public Dune::PDELab::ExplicitOneStepMethod<
T,
Dune::PDELab::OneStepGridOperator<GO0,GO1,false>,
LS,
TrlV,
TstV>
{
using Base = Dune::PDELab::ExplicitOneStepMethod<
T,
Dune::PDELab::OneStepGridOperator<GO0,GO1,false>,
LS,
TrlV,
TstV>;
public:
using Base::Base;
};
// Implicit
template<class T, class GO0, class GO1, class LS, class TrlV, class TstV, class TC>
class OneStepMethodLinear<
T,
Dune::PDELab::OneStepGridOperator<GO0,GO1,true>,
LS,
TrlV,
TstV,
TC>
: public Dune::PDELab::StationaryLinearProblemSolver<
Dune::PDELab::OneStepGridOperator<GO0,GO1,true>,
LS,
TrlV/*Not completeley sure about this one*/>
, public Dune::PDELab::OneStepMethod<
T,
Dune::PDELab::OneStepGridOperator<GO0,GO1,true>,
Dune::PDELab::StationaryLinearProblemSolver<
Dune::PDELab::OneStepGridOperator<GO0,GO1,true>,
LS,
TrlV/*Not completeley sure about this one*/>,
TrlV,
TstV/*,TC*/>
{
using IGO = Dune::PDELab::OneStepGridOperator<GO0,GO1,true>;
using Base1 = Dune::PDELab::StationaryLinearProblemSolver<
IGO,
LS,
TrlV/*Not completeley sure about this one*/>;
using Base2 = Dune::PDELab::OneStepMethod<T,IGO,Base1,TrlV,TstV/*,TC*/>;
public:
OneStepMethodLinear(
const Dune::PDELab::TimeSteppingParameterInterface<T>& method,
IGO& igo,
LS& ls
) : Base1(igo,ls,1e-6,1e-99,0)
, Base2(method,igo,*static_cast<Base1*>(this))
{}
using Base2::apply;
using Base2::result;
};
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_ONE_STEP_METHOD_LINEAR_HH
\ No newline at end of file
#ifndef DUNE_DORIE_STATISTICAL_MOMENTS_HH
#define DUNE_DORIE_STATISTICAL_MOMENTS_HH
namespace Dune {
namespace Dorie {
template<class GF>
std::vector<typename GF::Traits::DomainType>
statistical_moments(
const GF& gf,
unsigned int int_order,
unsigned int max_order = 2)
{
using Traits = typename GF::Traits;
using Domain = typename Traits::DomainType;
constexpr unsigned int dim = Traits::dimDomain;
assert(max_order>0);
std::vector<Domain> moment(max_order);
const auto& gv = gf.getGridView();
// Loop over the grid
for (const auto& cell : Dune::elements(gv,Dune::Partitions::interior))
{
auto geo = cell.geometry();
for (const auto& point : quadratureRule(geo, int_order))
{
const auto& x = point.position();
auto x_gl = geo.global(x);
typename Traits::RangeType y;
gf.evaluate(cell,x,y);
auto factor = point.weight() * geo.integrationElement(x);
Domain weight(1.);
for (unsigned int i = 0; i < max_order; ++i) {
for (unsigned int j = 0; j < dim; ++j)
{
weight[j] *= x_gl[j];
moment[i][j] += weight[j]*y*factor;
}
}
}
}
// Accumulate over processors
for (unsigned int i = 0; i < max_order; ++i)
moment[i] = gv.comm().sum(moment[i]);
return moment;
}
}
}
#endif
#ifndef DUNE_DORIE_TENSOR_TO_GLYPH_HH
#define DUNE_DORIE_TENSOR_TO_GLYPH_HH
#include <dune/common/fmatrix.hh>
namespace Dune{
namespace Dorie{
template<class RF, int dim>
Dune::FieldMatrix<RF,3,3> tensor_to_glyph(const Dune::FieldMatrix<RF,dim,dim>& y)
{
Dune::FieldMatrix<RF,3,3> y_glyph;
y_glyph[0][0] = y[0][0];
y_glyph[0][1] = y[0][1];
y_glyph[1][0] = y[1][0];
y_glyph[1][1] = y[1][1];
if constexpr (dim == 3){
y_glyph[0][2] = 0.;
y_glyph[1][2] = 0.;
y_glyph[2][2] = 0.;
y_glyph[1][2] = 0.;
y_glyph[0][2] = 0.;
}
return y_glyph;
}
}
}
#endif
\ No newline at end of file
......@@ -31,6 +31,10 @@ enum class OutputPolicy {None,EndOfRichardsStep,EndOfTransportStep};
*/
enum class AdaptivityPolicy {None,WaterFlux,SoluteFlux};
/// Controls how to deal with run-time checks
enum class CheckLevel { none, warn, error };
/// 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
......
......@@ -14,6 +14,7 @@
#include <dune/common/fvector.hh>
#include <dune/grid/common/gridview.hh>
#include <dune/dorie/common/util.hh>
namespace Dune {
namespace Dorie {
......@@ -116,6 +117,23 @@ inline bool is_none (const std::string& in)
return false;
}
/// Return a check level for a given string
inline CheckLevel to_check_lvl (const std::string check_lvl_str)
{
if (to_lower(check_lvl_str) == "none") {
return CheckLevel::none;
}
else if (check_lvl_str == "warn") {
return CheckLevel::warn;
}
else if (check_lvl_str == "error") {
return CheckLevel::error;
}
else {
DUNE_THROW(IOError, "Unknown check level");
}
}
/*-------------------------------------------------------------------------*//**
* @brief Transform a Yaml node with sequences to map
* @details This function check all the nodes within 'in' and transform every
......
#ifndef DUNE_DORIE_ADAPTIVITY_HANDLER_HH
#define DUNE_DORIE_ADAPTIVITY_HANDLER_HH
#ifndef DUNE_DORIE_RICHARDS_ADAPTIVITY_HANDLER_HH
#define DUNE_DORIE_RICHARDS_ADAPTIVITY_HANDLER_HH
#include <iostream>
#include <string>
......@@ -22,24 +22,25 @@
namespace Dune{
namespace Dorie{
/// Adaptivity base class. Does nothing.
template<typename Traits, typename GFS, typename Param, typename Boundary, typename U, int order>
class AdaptivityHandlerBase
class RichardsAdaptivityHandler
{
private:
using Grid = typename Traits::Grid;
using GV = typename Traits::GV;
using RF = typename Traits::RF;
using Grid = typename Traits::Grid;
using GV = typename Traits::GV;
using RF = typename Traits::RF;
public:
AdaptivityHandlerBase () = default;
virtual ~AdaptivityHandlerBase () = default;
RichardsAdaptivityHandler () = default;
virtual ~RichardsAdaptivityHandler () = default;
/// Do nothing.
virtual void mark_grid (Grid& grid, GV& gv, GFS& gfs, const Param& param, const Boundary& boundary, const RF time, U& u) { }
/// Do nothing.
virtual void mark_grid (Grid& grid, GV& gv, GFS& gfs, const Param& param, const Boundary& boundary, const RF time, U& u) { }
/// Do nothing.
virtual void adapt_grid (Grid& grid, GFS& gfs, U& u) { }
/// Do nothing.
virtual void adapt_grid (Grid& grid, GFS& gfs, U& u) { }
};
/// Specialized SimulationBase class providing functions and members for adaptive grid refinement
......@@ -49,7 +50,7 @@ template<typename Traits,
typename Boundary,
typename U,
int order>
class AdaptivityHandler : public AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U,order>
class WaterFluxAdaptivityHandler : public RichardsAdaptivityHandler<Traits,GFS,Param,Boundary,U,order>
{
private:
using RF = typename Traits::RF;
......@@ -61,7 +62,7 @@ private:
/// Adaptivity GFS
using AGFS = typename AGFSHelper::Type;
/// Error estimator local operator
using ESTLOP = Dune::Dorie::Operator::FluxErrorEstimator<Traits,Param,Boundary>;
using ESTLOP = Dune::Dorie::Operator::WaterFluxErrorEstimator<Traits,Param,Boundary>;
/// Empty constraints container
using NoTrafo = Dune::PDELab::EmptyTransformation;
/// Solution vector type
......@@ -90,8 +91,8 @@ public:
* \param _inifile Parameter file parser
* \param grid Grid to adapt (reference is not saved)
*/
AdaptivityHandler (Dune::ParameterTree& _inifile, Grid& grid)
: AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U,order>(),
WaterFluxAdaptivityHandler (Dune::ParameterTree& _inifile, Grid& grid)
: RichardsAdaptivityHandler<Traits,GFS,Param,Boundary,U,order>(),
inifile(_inifile),
maxLevel(inifile.get<int>("adaptivity.maxLevel")),
minLevel(inifile.get<int>("adaptivity.minLevel")),
......@@ -220,7 +221,7 @@ public:
* define which create() function is enabled at compile time.
*/
template<typename Traits, typename GFS, typename Param, typename Boundary, typename U, int order>
struct AdaptivityHandlerFactory
struct RichardsAdaptivityHandlerFactory
{
private:
......@@ -246,7 +247,7 @@ private:
using Grid = typename Traits::Grid;
using AHB = AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U,order>;
using AHB = RichardsAdaptivityHandler<Traits,GFS,Param,Boundary,U,order>;
Dune::ParameterTree& inifile;
Grid& grid;
......@@ -257,17 +258,17 @@ public:
/*---------------------------------------------------------------------*//**
* @brief Create the factory, taking the parameters for building an
* AdaptivityHandler
* RichardsAdaptivityHandler
*
* @param _inifile Parameter file parser
* @param _grid Grid to adapt (reference is not saved)
*/
AdaptivityHandlerFactory (Dune::ParameterTree& _inifile, Grid& _grid) :
RichardsAdaptivityHandlerFactory (Dune::ParameterTree& _inifile, Grid& _grid) :
inifile(_inifile),
grid(_grid)
{ }
/// Create a dummy AdaptivityHandler
/// Create a dummy RichardsAdaptivityHandler
/** Assert that adaptive grid (UG) is used.
*/
template<bool active = enabled>
......@@ -276,15 +277,15 @@ public:
return std::make_unique<AHB>();
}
/// Create a true AdaptivityHandler
/// Create a true RichardsAdaptivityHandler
template<bool active = enabled>
std::enable_if_t<active,std::unique_ptr<AHB>> create ()
{
std::unique_ptr<AHB> p;
if (inifile.hasSub("adaptivity") &&
inifile.get<std::string>("adaptivity.policy") != "none")
p = std::make_unique<AdaptivityHandler<Traits,GFS,Param,Boundary,U,order>>(inifile,grid);
else
inifile.get<std::string>("adaptivity.policy") == "waterFlux")
p = std::make_unique<WaterFluxAdaptivityHandler<Traits,GFS,Param,Boundary,U,order>>(inifile,grid);
else
p = std::make_unique<AHB>();
return p;
}
......@@ -293,5 +294,4 @@ public:
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_ADAPTIVITY_HANDLER_HH
#endif // DUNE_DORIE_RICHARDS_ADAPTIVITY_HANDLER_HH
......@@ -28,7 +28,7 @@ namespace Dune {
* Skeleton integral is mostly the same as in the DG operator.
*/
template<typename Traits, typename Parameter, typename Boundary>
class FluxErrorEstimator
class WaterFluxErrorEstimator
: public Dune::PDELab::LocalOperatorDefaultFlags
{
......@@ -61,16 +61,16 @@ namespace Dune {
RF penalty_factor;
RF time;
FluxErrorEstimator (const Dune::ParameterTree& config,
WaterFluxErrorEstimator (const Dune::ParameterTree& config,
const Parameter& param_, const Boundary& boundary_,
RichardsDGMethod::Type method_ = RichardsDGMethod::SIPG,
DGMethod::Type method_ = DGMethod::SIPG,
int intorderadd_ = 2, int quadrature_factor_ = 2)
: param(param_), boundary(boundary_),
intorderadd(intorderadd_), quadrature_factor(quadrature_factor_),
penalty_factor(config.get<RF>("numerics.penaltyFactor")),
time(0.0)
{
if (method_ == RichardsDGMethod::OOB)
if (method_ == DGMethod::OOB)
penalty_factor = 0.0;
}
......
This diff is collapsed.
#include <dune/pdelab/function/scalarscaled.hh>
#include <dune/dorie/model/richards/richards.hh>
#include <dune/dorie/common/grid_function_integrator.hh>
#include <dune/dorie/common/statistical_moments.hh>
#include <limits>
#include <cmath>
namespace Dune{
namespace Dorie{
......@@ -94,6 +100,16 @@ RichardsSimulation<Traits>::RichardsSimulation (
inifile.get<std::string>("output.fileName"),
inifile.get<std::string>("output.outputPath"),
"./");
// Configure statistics file
std::ofstream stats_file;
stats_file.open(inifile.get<std::string>("output.fileName") + ".dat", std::ofstream::trunc);
stats_file << "# Statistics of the Richards Equation generated by DORiE.\n";
stats_file << "#\n";
stats_file << "# run_time\t time\t mass\t boundary_flux\t";
stats_file << "mean_cw_x\t mean_cw_y\t mean_cw_z\t";
stats_file << "variance_cw_x\t variance_cw_y\t variance_cw_z\n";
stats_file.close();
}
// --- Adaptivity Setup --- //
......@@ -229,6 +245,8 @@ void RichardsSimulation<Traits>::step()
// invalidate cache for water flux reconstruction
cache_water_flux_gf_rt.reset();
time_steps++;
if (this->output_policy() == OutputPolicy::EndOfRichardsStep)
write_data();
}
......@@ -266,15 +284,102 @@ void RichardsSimulation<Traits>::post_adapt_grid()
template<typename Traits>
void RichardsSimulation<Traits>::write_data () const
{
if (vtkwriter)
const auto time = controller->getTime();
// create adapters
auto head = get_matric_head();
auto wflux = std::make_shared<const GFWaterFlux>(gfs, u, fparam);
auto cond = std::make_shared<const GFConductivity>(gv, fparam);
auto wc = std::make_shared<const GFWaterContent>(head, gv, fparam);
auto sat = std::make_shared<const GFSaturation>(head, gv, fparam);
// Compute mass of WC
RF mass = gf_integrator(*wc,2*order+2);
// Compute statistical moments of WC
using PDF = Dune::PDELab::ScalarScaledGridFunctionAdapter<const GFWaterContent>;
PDF pdf(1./mass,*wc);
auto moments = statistical_moments(pdf,order*4,2);
auto moments_squared = moments[0];
for (int i = 0; i < dim; ++i) moments_squared[i]*=moments_squared[i];
const auto variance = moments[1] - moments_squared;
Dune::FieldVector<double,3> moments0_3d(0.), variance_3d(0.);
for (int i = 0; i < dim; ++i)
{
// create adapters
auto head = get_matric_head();
auto wflux = std::make_shared<const GFWaterFlux>(gfs, u, fparam);
auto cond = std::make_shared<const GFConductivity>(gv, fparam);
auto wc = std::make_shared<const GFWaterContent>(head, gv, fparam);
auto sat = std::make_shared<const GFSaturation>(head, gv, fparam);
moments0_3d[i] = moments[0][i];
variance_3d[i] = variance[i];
}
RF boundary_flux = std::numeric_limits<RF>::quiet_NaN();;
if constexpr (enable_rt_engine) {
if (enable_fluxrc) {
auto wfluxr = get_water_flux_reconstructed();
boundary_flux = boundary_gf_integrator(*wfluxr,2*flux_order+2);
}
}
if (time_steps==0) cache_time = time;
if (time_steps==0) cache_mass = mass;
if (time_steps==0) cache_boundary_flux = boundary_flux;
// hard check on mass
if (std::isnan(mass)) {
this->_log->error("Mass of water is NaN: PDE solver has diverged.");
DUNE_THROW(MathError,"PDE solver has diverged!");
}
// check relative change of mass
if (Dune::FloatCmp::ne(mass,0.))
{
// change in time
auto dt = time - cache_time;
// change of mass by PDE solution
auto mass_dt = (mass - cache_mass);
// change of mass by boundary conditions
mass_dt += 0.5*dt*(boundary_flux + cache_boundary_flux);
// relative mass change
auto rel_mass_dt = mass_dt/mass;
auto mass_tol = inifile.get<double>("conservationOfMass.checkTol");
auto check_mass_str = inifile.get<std::string>("conservationOfMass.check");
auto check_mass_lvl = Dorie::to_check_lvl(check_mass_str);
const bool passed = rel_mass_dt < mass_tol;
if (passed) {
this->_log->trace(" Conservation of mass does not hold for"
" {}% of mass",
rel_mass_dt);
}
else if (not passed && check_mass_lvl == CheckLevel::warn) {
this->_log->warn("Conservation of mass above tolerance:"
" {}% of mass is not conserved.",
rel_mass_dt);
}
else if (not passed && check_mass_lvl == CheckLevel::error) {
this->_log->error("Conservation of mass above tolerance:"
" {}% of mass is not conserved.",
rel_mass_dt);
DUNE_THROW(MathError, "Conservation of mass does not obey tolerance");
}
}
// update cache
cache_time = time;
cache_mass = mass;
cache_boundary_flux = boundary_flux;
// Configure statistics file
std::ofstream stats_file;
stats_file.open(inifile.get<std::string>("output.fileName") + ".dat", std::ofstream::app);
stats_file << timer.elapsed() << "\t" << time << "\t" << mass << "\t";
stats_file << boundary_flux << "\t" << moments0_3d << "\t" << variance_3d;
stats_file << std::endl;
stats_file.close();
if (vtkwriter)
{
if (inifile.get<bool>("output.vertexData")) {
vtkwriter->addVertexData(head,"head");
vtkwriter->addVertexData(wflux,"flux");
......@@ -302,7 +407,6 @@ void RichardsSimulation<Traits>::write_data () const
}
try{
const auto time = controller->getTime();
this->_log->trace("Writing solution at time {:.2e}", time);
vtkwriter->write(time, output_type);
}
......
......@@ -15,6 +15,7 @@
#include <dune/dorie/common/simulation_base.hh>
#include <dune/dorie/common/util.hh>
#include <dune/dorie/common/utility.hh>
#include <dune/dorie/common/interpolator.hh>
#include <dune/dorie/common/time_controller.hh>
#include <dune/dorie/common/grid_creator.hh>
......@@ -163,9 +164,9 @@ struct RichardsSimulationTraits : public BaseTraits
/// Custom VTK output writer
using Writer = Dune::Dorie::GridFunctionVTKSequenceWriter<GV>;
/// Grid Adaptivity Handler base class
using AdaptivityHandler = Dune::Dorie::AdaptivityHandlerBase<BaseTraits,GFS,FlowParameters,FlowBoundary,U,order>;
using AdaptivityHandler = Dune::Dorie::RichardsAdaptivityHandler<BaseTraits,GFS,FlowParameters,FlowBoundary,U,order>;
/// Factory for creating the AdaptivityHandler
using AdaptivityHandlerFactory = Dune::Dorie::AdaptivityHandlerFactory<BaseTraits,GFS,FlowParameters,FlowBoundary,U,order>;
using AdaptivityHandlerFactory = Dune::Dorie::RichardsAdaptivityHandlerFactory<BaseTraits,GFS,FlowParameters,FlowBoundary,U,order>;
};
/*-------------------------------------------------------------------------*//**
......@@ -190,6 +191,8 @@ private:
static constexpr int order = Traits::order;
static constexpr int flux_order = Traits::flux_order;
static constexpr bool enable_rt_engine = Traits::enable_rt_engine;
using TimeField = typename Traits::TimeField;
using RF = typename Traits::RF;
using Grid = typename Traits::Grid;
using GV = typename Traits::GV;
......@@ -316,12 +319,18 @@ protected:
std::unique_ptr<Writer> vtkwriter;
std::unique_ptr<AdaptivityHandler> adaptivity;
mutable std::shared_ptr<GFFluxReconstruction> cache_water_flux_gf_rt;
Dune::Timer timer;
RF time_before;
RF dt_before;
const bool enable_fluxrc;
unsigned int time_steps;
mutable std::shared_ptr<GFFluxReconstruction> cache_water_flux_gf_rt;
mutable RF cache_mass;
mutable RF cache_boundary_flux;
mutable TimeField cache_time;
public:
/*------------------------------------------------------------------------*//**
......
This diff is collapsed.
add_library(dorie-transport STATIC
sim_yasp_2_0.cc
sim_yasp_2_1.cc
sim_yasp_2_2.cc
sim_ug_2_0.cc
sim_yasp_2_2.cc
sim_ug_2_1.cc
sim_yasp_2_3.cc
sim_ug_2_2.cc
sim_yasp_3_0.cc
sim_yasp_3_1.cc
sim_ug_2_2.cc
sim_ug_2_3.cc
sim_yasp_3_2.cc
sim_yasp_3_3.cc
sim_ug_2_3.cc
sim_ug_3_0.cc
sim_ug_3_1.cc
sim_ug_3_2.cc
......
......@@ -12,8 +12,8 @@ namespace Dorie{
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,1,1>>; *
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,2,1>>; *
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,3,1>>; *
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,1,1>>; +
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,2,1>>; +
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,1,1>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,2,1>>;
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,3,1>>; *
} // namespace Dorie
......
......@@ -12,8 +12,8 @@ namespace Dorie{
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,1,2>>; *
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,2,2>>; *
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,3,2>>; *
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,1,2>>; +
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,2,2>>; +
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,1,2>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,2,2>>;
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,3,2>>; *
} // namespace Dorie
......
......@@ -12,8 +12,8 @@ namespace Dorie{
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,1,3>>; *
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,2,3>>; *
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,3,3>>; *
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,1,3>>; +
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,2,3>>; +
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,1,3>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,2,3>>;
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,3,3>>; *
} // namespace Dorie
......
......@@ -9,9 +9,9 @@
namespace Dune{
namespace Dorie{
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,1,1>>; +
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,2,1>>; +
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,3,1>>; +
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,1,1>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,2,1>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,3,1>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -9,9 +9,9 @@
namespace Dune{
namespace Dorie{
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,1,2>>; +
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,2,2>>; +
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,3,2>>; +
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,1,2>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,2,2>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,3,2>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -9,9 +9,9 @@
namespace Dune{
namespace Dorie{
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,1,3>>; +
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,2,3>>; +
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,3,3>>; +
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,1,3>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,2,3>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,3,3>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -9,8 +9,8 @@
namespace Dune{
namespace Dorie{
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,1,1>>; +
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,2,1>>; +
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,1,1>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,2,1>>;
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,3,1>>; *
} // namespace Dorie
......
......@@ -9,8 +9,8 @@
namespace Dune{
namespace Dorie{
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,1,2>>; +
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,2,2>>; +
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,1,2>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,2,2>>;
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,3,2>>; *
} // namespace Dorie
......
......@@ -9,8 +9,8 @@
namespace Dune{
namespace Dorie{
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,1,3>>; +
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,2,3>>; +
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,1,3>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,2,3>>;
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,3,3>>; *
} // namespace Dorie
......
This diff is collapsed.
......@@ -101,10 +101,6 @@ private:
using WaterFlux = typename GFWaterFlux::Traits::RangeType;
using WaterContent = typename GFWaterContent::Traits::RangeType;
// Extract world dimension
enum {dim = GFWaterContent::Traits::dimDomain};
using Tensor = Dune::FieldMatrix<RF,dim,dim>;
public:
/**
......@@ -505,9 +501,9 @@ public:
public:
TransportFVTemporalOperator()
: _time(0.)
{ }
TransportFVTemporalOperator()
: _time(0.)
{ }
/*---------------------------------------------------------------------*//**
* @brief Volume integral depending on test and ansatz functions
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -11,4 +11,4 @@ spatial_resolution_front_we -1
spatial_resolution_back_sn -1
spatial_resolution_back_we -1
number_BC_change_times 1
0 neumann -1e-6 dirichlet 0
\ No newline at end of file
0 neumann -5.55e-6 dirichlet 0
\ No newline at end of file
......@@ -6,9 +6,11 @@ _evaluation = const
_model = transport
_quantity = solute
_tol = 1E-8
_diff = conv, diff | expand diff
_param = sand, silt, layered | expand param
__name = const_solute_{_diff}_{_param}
_param = sand, layered | expand param
_rorder = 1, 2 | expand rorder
_ycells = 160, 40 | expand rorder
_torder = 0, 1 | expand torder
__name = const_solute_{_param}_T{_torder}_R{_rorder}
[simulation]
mode = richards+transport
......@@ -16,11 +18,11 @@ mode = richards+transport
[grid]
gridType = rectangular
initialLevel = 0
cells = 1 160
cells = 1 {_ycells}
[grid.mapping]
file = "{_asset_path}/maps/cell_ids.h5"
volume = 0, 1, ode_layered_160 | expand param
volume = 0, ode_layered_{_ycells} | expand param
[richards]
output.policy = none
......@@ -33,7 +35,7 @@ time.end = 1E6
time.maxTimestep = 1E6
time.startTimestep = 1E4
numerics.FEorder = 1
numerics.FEorder = {_rorder}
[transport]
output.fileName = {__name}
......@@ -46,10 +48,10 @@ boundary.file = "{_asset_path}/bcs/solute_2d_const.dat"
parameters.file = "{_asset_path}/param/transport_param.yml"
parameters.effHydrDips = 0, 2E-9 | expand diff
time.end = 1E6
time.maxTimestep = 1E6
time.startTimestep = 1E4
numerics.timestepMethod = implicit_euler
numerics.FEorder = {_torder}
\ No newline at end of file