...
 
Commits (2)
  • Lukas Riedel's avatar
    Add upwinding options to Richards FV local operator · 2c263469
    Lukas Riedel authored
    Run ODE tests now with upwinding options of FV and DG solver
    where a change of the residual tolerance limit was not necessary.
    
    Several adjustments were made to make applying local operator options
    easier:
    * Remove 'EXPERIMENTAL_DG_FEATURES' option. Settings in config file
      now apply by default.
    * Move previously experimental options to 'richards.numerics' section.
    * Rename settings 'method' to 'DGMethod' and 'weights' to 'DGWeights'.
    * Fix typo: Rename OOB to OBB scheme (Oden, Babuska, Baumann).
    2c263469
  • Lukas Riedel's avatar
    Merge branch '153-add-upwinding-option-to-richards-fv-solver' into 'master' · 373885e6
    Lukas Riedel authored
    Resolve "Add upwinding option to Richards FV solver"
    
    Closes #153
    
    See merge request !161
    373885e6
......@@ -46,6 +46,7 @@
* Define compile-time settings via CPP macros !144
* [Google Test](https://github.com/google/googletest) unit test framework
as Git Submodule !159
* Upwinding options for Richards finite volume local operator !161
### Changed
* `Simulation` is renamed `RichardsSimulation` and moved to
......
......@@ -14,12 +14,6 @@ endif()
# add extra flags to debug compiler flags
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall")
# option to change DG scheme via config file
option(EXPERIMENTAL_DG_FEATURES
"Enable experimental DG settings through the config file"
OFF
)
#
if(NOT (dune-common_DIR OR dune-common_ROOT OR
"${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*"))
......
......@@ -199,16 +199,6 @@ If you installed [Anaconda](https://conda.io/docs/user-guide/install/download.ht
path chosen as installation prefix when configuring HDF5.
### Experimental Features
The local operator implementing Richards equation's discretization supports
multiple scheme settings. Setting these via the config file is disabled by
default. You can enable this feature by reconfiguring DORiE with the CMake flag
`-DEXPERIMENTAL_DG_FEATURES=On`, and rebuilding it.
The configuration settings in the section `[dg.experimental]` will then override
the default settings.
### Recommended Third-Party Software
The following software packages are cross-platform, so you should be able to find a release that fits your operating system:
......
......@@ -231,13 +231,6 @@ 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="FEorder">
<definition> Polynomial order of the DG method used. Setting this value
to 0 (zero) selects the finite volume (FV) solver (only compatible to
......@@ -245,6 +238,52 @@ adding an empty line, make text **bold** or ``monospaced``.
</definition>
<suggestion> 1 </suggestion>
<values> 0, 1, 2, 3 </values>
<comment> Select '0' for the finite volume (FV) solver </comment>
</parameter>
<parameter name="upwinding">
<definition> Upwinding method for skeleton terms. Upwinding typically
increases numeric stability while reducing accuracy.
**semiUpwind:** Apply upwinding to conductivity factor (only).
**fullUpwind:** Apply upwinding to conductivity.
**Not recommended for DG**.
</definition>
<values> none, semiUpwind, fullUpwind </values>
<suggestion> none </suggestion>
<comment> Choose upwinding type: 'none', 'semiUpwind', 'fullUpwind'
</comment>
</parameter>
<parameter name="DGMethod">
<definition> DG discretization method for skeleton terms.
**SIPG:** Symmetric Interior Penalty
**NIPG:** Non-Symmetric Interior Penalty
**OBB:** Oden, Babuska, Baumann: no penalty term
**IIP:** Incomplete Interior Penalty: no symmetry term
</definition>
<values> SIPG, NIPG, OBB, IIP </values>
<suggestion> SIPG </suggestion>
</parameter>
<parameter name="DGWeights">
<definition> Apply harmonic weighting to skeleton term contributions
in DG.
</definition>
<values> true, false </values>
<suggestion> true </suggestion>
</parameter>
<parameter name="penaltyFactor">
<definition> Penalty factor to be used in the Discontinuous Galerkin scheme
</definition>
<values> float </values>
<suggestion> 10 </suggestion>
</parameter>
</category>
......@@ -328,41 +367,4 @@ adding an empty line, make text **bold** or ``monospaced``.
</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>
......@@ -22,8 +22,3 @@ endforeach()
# Global target
add_custom_target("dorie" DEPENDS richards transport)
# enable setting operator scheme from config file
if(EXPERIMENTAL_DG_FEATURES)
target_compile_definitions("dorie" PUBLIC -DEXPERIMENTAL_DG_FEATURES)
endif()
......@@ -15,7 +15,8 @@
#include <dune/pdelab/localoperator/defaultimp.hh>
#include <dune/pdelab/finiteelement/localbasiscache.hh>
#include <dune/dorie/model/richards/local_operator_DG.hh>
#include <dune/dorie/common/typedefs.hh>
#include <dune/dorie/model/richards/local_operator_cfg.hh>
namespace Dune {
namespace Dorie {
......@@ -63,14 +64,14 @@ namespace Dune {
FluxErrorEstimator (const Dune::ParameterTree& config,
const Parameter& param_, const Boundary& boundary_,
RichardsDGMethod::Type method_ = RichardsDGMethod::SIPG,
RichardsDGMethod method_ = RichardsDGMethod::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_ == RichardsDGMethod::OBB)
penalty_factor = 0.0;
}
......@@ -310,4 +311,4 @@ namespace Dune {
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_ERROR_INDICATOR_HH
\ No newline at end of file
#endif // DUNE_DORIE_ERROR_INDICATOR_HH
......@@ -2,10 +2,8 @@
#ifndef DUNE_DORIE_RICHARDS_OPERATOR_HH
#define DUNE_DORIE_RICHARDS_OPERATOR_HH
#include <string>
#include <tuple>
#include <map>
#include <vector>
#include <memory>
#include <dune/common/parametertree.hh>
#include <dune/common/exceptions.hh>
......@@ -22,40 +20,12 @@
#include <dune/pdelab/finiteelement/localbasiscache.hh>
#include <dune/dorie/common/typedefs.hh>
#include <dune/dorie/model/richards/local_operator_cfg.hh>
namespace Dune {
namespace Dorie {
namespace Operator {
//! Interior penalty type
struct RichardsDGMethod
{
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)
};
};
//! Upwinding type
struct RichardsDGUpwinding
{
enum Type {
fullUpwind, //!< Upwind flux and conductiviy
semiUpwind, //!< Only upwind conductivity factor
none //!< No upwinding
};
};
//! Gradient weighting
struct RichardsDGWeights
{
enum Type {
weightsOn, //!< harmonic weighting of flux contributions
weightsOff //!< no weighting of flux contributions
};
};
//! Operator internal BC handling
struct BCType
{
......@@ -65,62 +35,6 @@ struct BCType
};
};
/// Read experimental operator settings from a parameter file
/* \return Tuple of RichardsDGMethod, RichardsDGUpwinding, and
* RichardsDGWeights to be inserted into the local operator constructor.
*/
inline auto read_experimental_operator_settings(
const Dune::ParameterTree &inifile)
-> std::tuple<Operator::RichardsDGMethod::Type,
Operator::RichardsDGUpwinding::Type,
Operator::RichardsDGWeights::Type>
{
const auto log = Dorie::get_logger(log_richards);
log->debug("Reading experimental DG operator settings:");
RichardsDGMethod::Type method;
RichardsDGUpwinding::Type upwinding;
RichardsDGWeights::Type weights;
const auto method_str = inifile.get<std::string>("dg.experimental.method");
log->debug(" DG method: {}", method_str);
if (method_str == "SIPG")
method = RichardsDGMethod::SIPG;
else if (method_str == "NIPG")
method = RichardsDGMethod::NIPG;
else if (method_str == "OOB")
method = RichardsDGMethod::OOB;
else if (method_str == "IIP")
method = RichardsDGMethod::IIP;
else {
log->error("Unknown DG method: {}", method_str);
DUNE_THROW(Dune::IOError, "Unknown DG method");
}
const auto upwinding_str
= inifile.get<std::string>("dg.experimental.upwinding");
log->debug(" Upwinding scheme: {}", upwinding_str);
if (upwinding_str == "none")
upwinding = RichardsDGUpwinding::none;
else if (upwinding_str == "semiUpwind")
upwinding = RichardsDGUpwinding::semiUpwind;
else if (upwinding_str == "fullUpwind")
upwinding = RichardsDGUpwinding::fullUpwind;
else {
log->error("Unknown upwinding scheme: {}", upwinding_str);
DUNE_THROW(Dune::IOError, "Unknown upwinding scheme");
}
const auto weights_str = inifile.get<bool>("dg.experimental.weights");
log->debug(" Harmonic weights at interfaces: {}", weights_str);
if (weights_str)
weights = RichardsDGWeights::weightsOn;
else
weights = RichardsDGWeights::weightsOff;
return std::make_tuple(method, upwinding, weights);
}
/// A local operator for solving the Richards equation using upwinding.
/** Written by Ole Klein, IWR, Heidelberg and modified by the
* DORiE developers, IUP, Heidelberg \n
......@@ -167,9 +81,9 @@ protected:
std::shared_ptr<const Parameter> param; //!< class providing equation parameters
std::shared_ptr<const Boundary> boundary; //!< class providing boundary conditions
std::shared_ptr<const SourceTerm> sourceTerm; //!< class providing source term information
const RichardsDGMethod::Type method; //!< Interior penalty type
const RichardsDGUpwinding::Type upwinding; //!< DG method weights switch
const RichardsDGWeights::Type weights; //!< Gradient weighting
const RichardsDGMethod method; //!< Interior penalty type
const RichardsUpwinding upwinding; //!< DG method weights switch
const RichardsDGWeights weights; //!< Gradient weighting
const int intorderadd; //!< integration order addend
const int quadrature_factor; //!< factor to FEM integration order
......@@ -198,9 +112,9 @@ public:
* \param param_ Object providing equation parameters
* \param boundary_ Object prividing boundary conditions
* \param sourceTerm_ Object providing source term information
* \param method_ DG method type \see RichardsDGMethod::Type
* \param upwinding_ Upwinding type \see RichardsDGUpwinding::Type
* \param weights_ DG weigths switch \see RichardsDGWeights::Type
* \param method_ DG method type \see RichardsDGMethod
* \param upwinding_ Upwinding type \see RichardsUpwinding
* \param weights_ DG weigths switch \see RichardsDGWeights
* \param intorderadd_ Addend to integration order
* \param quadrature_factor_ Factor to FEM integration order
*/
......@@ -209,9 +123,9 @@ public:
std::shared_ptr<const Parameter> param_,
std::shared_ptr<const Boundary> boundary_,
std::shared_ptr<const SourceTerm> sourceTerm_,
const RichardsDGMethod::Type method_ = RichardsDGMethod::SIPG,
const RichardsDGUpwinding::Type upwinding_ = RichardsDGUpwinding::none,
const RichardsDGWeights::Type weights_ = RichardsDGWeights::weightsOn,
const RichardsDGMethod method_ = RichardsDGMethod::SIPG,
const RichardsUpwinding upwinding_ = RichardsUpwinding::none,
const RichardsDGWeights weights_ = RichardsDGWeights::weightsOn,
int intorderadd_ = 2,
int quadrature_factor_ = 2
) :
......@@ -231,12 +145,12 @@ public:
cache(20)
{
theta = 0.; // IIP
if (method == RichardsDGMethod::NIPG || method == RichardsDGMethod::OOB)
if (method == RichardsDGMethod::NIPG || method == RichardsDGMethod::OBB)
theta = 1.;
else if (method == RichardsDGMethod::SIPG)
theta = -1.;
if (method == RichardsDGMethod::OOB)
if (method == RichardsDGMethod::OBB)
penalty_factor = 0.;
}
......@@ -527,13 +441,13 @@ public:
RF omega_s = 0.5, omega_n = 0.5;
if (weights == RichardsDGWeights::weightsOn)
{
if (upwinding == RichardsDGUpwinding::none)
if (upwinding == RichardsUpwinding::none)
{
omega_s = cond_n / ((cond_s) + (cond_n));
omega_n = cond_s / ((cond_s) + (cond_n));
}
else if (upwinding == RichardsDGUpwinding::semiUpwind
|| upwinding == RichardsDGUpwinding::fullUpwind)
else if (upwinding == RichardsUpwinding::semiUpwind
|| upwinding == RichardsUpwinding::fullUpwind)
{
omega_s = satCond_n / (satCond_s + satCond_n);
omega_n = satCond_s / (satCond_s + satCond_n);
......@@ -560,7 +474,7 @@ public:
numFlux_s *= satCond_s;
numFlux_n *= satCond_n;
if (upwinding == RichardsDGUpwinding::none)
if (upwinding == RichardsUpwinding::none)
{
numFlux_s *= relCond_s;
numFlux_n *= relCond_n;
......@@ -569,8 +483,8 @@ public:
// Upwind in relative conductivity
RF relCond_upwind=0.0, satCond_upwind=0.0;
if (upwinding == RichardsDGUpwinding::semiUpwind
|| upwinding == RichardsDGUpwinding::fullUpwind)
if (upwinding == RichardsUpwinding::semiUpwind
|| upwinding == RichardsUpwinding::fullUpwind)
{
RF cond_upwind;
if (numFlux > 0)
......@@ -593,21 +507,21 @@ public:
// diffusion term
// consistency term
// + penalty term
if (upwinding == RichardsDGUpwinding::none)
if (upwinding == RichardsUpwinding::none)
{
for (unsigned int i = 0; i < lfsv_s.size(); i++)
r_s.accumulate(lfsv_s, i, numFlux * phiv_s[i] * factor);
for (unsigned int i = 0; i < lfsv_n.size(); i++)
r_n.accumulate(lfsv_n, i, - numFlux * phiv_n[i] * factor);
}
else if (upwinding == RichardsDGUpwinding::semiUpwind)
else if (upwinding == RichardsUpwinding::semiUpwind)
{
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, relCond_upwind * numFlux * phiv_s[i] * factor);
for (unsigned int i = 0; i<lfsv_n.size(); i++)
r_n.accumulate(lfsv_n,i, - relCond_upwind * numFlux * phiv_n[i] * factor);
}
else if (upwinding == RichardsDGUpwinding::fullUpwind)
else if (upwinding == RichardsUpwinding::fullUpwind)
{
// Upwind numerical flux
double numFlux_upwind = (numFlux > 0) ? numFlux_s : numFlux_n;
......@@ -628,21 +542,21 @@ public:
// discrete gradient sign: Negative for lifting
double dg_sign = (lopcase == LOPCase::RTVolume) ? -1. : 1.;
if (upwinding == RichardsDGUpwinding::none)
if (upwinding == RichardsUpwinding::none)
{
for (unsigned int i = 0; i < lfsv_s.size(); i++)
r_s.accumulate(lfsv_s, i, dg_sign * theta * omega_s * (gradphiv_s[i] * normal) * relCond_s * satCond_s * jump * factor);
for (unsigned int i = 0; i < lfsv_n.size(); i++)
r_n.accumulate(lfsv_n, i, dg_sign * theta * omega_n * (gradphiv_n[i] * normal) * relCond_n * satCond_n * jump * factor);
}
else if (upwinding == RichardsDGUpwinding::semiUpwind)
else if (upwinding == RichardsUpwinding::semiUpwind)
{
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, dg_sign * theta * omega_s * (gradphiv_s[i] * normal) * relCond_upwind * satCond_s * jump * factor);
for (unsigned int i = 0; i<lfsv_n.size(); i++)
r_n.accumulate(lfsv_n,i, dg_sign * theta * omega_n * (gradphiv_n[i] * normal) * relCond_upwind * satCond_n * jump * factor);
}
else if (upwinding == RichardsDGUpwinding::fullUpwind)
else if (upwinding == RichardsUpwinding::fullUpwind)
{
for (unsigned int i = 0; i < lfsv_s.size(); i++)
r_s.accumulate(lfsv_s, i, dg_sign * theta * (gradphiv_s[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
......@@ -813,15 +727,15 @@ public:
// compute conductivity factor
RF relCond;
if (upwinding == RichardsDGUpwinding::semiUpwind
|| upwinding == RichardsDGUpwinding::fullUpwind)
if (upwinding == RichardsUpwinding::semiUpwind
|| upwinding == RichardsUpwinding::fullUpwind)
{
if (numFlux > 0)
relCond = relCond_s;
else
relCond = relCond_n;
}
else // (upwinding == RichardsDGUpwinding::none)
else // (upwinding == RichardsUpwinding::none)
{
// harmonic average of conductivity factors
relCond = 2 * relCond_s * relCond_n / (relCond_s + relCond_n);
......
......@@ -10,6 +10,7 @@
#include <dune/pdelab/localoperator/defaultimp.hh>
#include <dune/dorie/common/typedefs.hh>
#include <dune/dorie/model/richards/local_operator_cfg.hh>
namespace Dune {
namespace Dorie {
......@@ -60,9 +61,11 @@ public:
*/
RichardsFVSpatialOperator(
const std::shared_ptr<const Parameter>& param,
const std::shared_ptr<const Boundary>& boundary
const std::shared_ptr<const Boundary>& boundary,
const RichardsUpwinding upwinding = RichardsUpwinding::none
) : _param(param)
, _boundary(boundary)
, _upwinding(upwinding)
, _time(0.)
{}
......@@ -126,10 +129,21 @@ public:
const auto center_position_i_g = geo_i.center();
const auto center_position_o_g = geo_o.center();
// Distance between the two entity centers
const auto center_position_diff = center_position_i_g
- center_position_o_g;
const auto distance = center_position_diff.two_norm();
// Inside/outside solute value
const RangeU u_i = x_i(lfsu_i,0);
const RangeU u_o = x_o(lfsu_o,0);
// Finite difference of u between the two entities
RangeU dudn = (u_o-u_i) / distance;
// Update gradient with gravity vector
dudn -= _param->gravity() * normal_f;
// bind parameterization and retrieve functions
_param->bind(entity_i);
const auto saturation_f_i = _param->saturation_f();
......@@ -143,28 +157,49 @@ public:
const RangeU cond_i = conductivity_f_i(saturation_f_i(u_i));
const RangeU cond_o = conductivity_f_o(saturation_f_o(u_o));
// Harmonic average for interface conductivity
const RangeU cond_f = 2.0/( 1.0/(cond_i + 1E-30)
+ 1.0/(cond_o + 1E-30) );
// Distance between the two entity centers
const auto center_position_diff = center_position_i_g - center_position_o_g;
const auto distance = center_position_diff.two_norm();
// Finite difference of u between the two entities
RangeU dudn = (u_o-u_i)/distance;
// Update gradient with gravity vector
dudn -= _param->gravity()*normal_f;
// Interface conductivity factor (Upwinding applies)
RangeU cond_f(0.0);
if (_upwinding == RichardsUpwinding::none)
{
// Harmonic average for interface conductivity
cond_f = 2.0/( 1.0/(cond_i + 1E-30)
+ 1.0/(cond_o + 1E-30) );
}
else if (_upwinding == RichardsUpwinding::semiUpwind)
{
// Harmonic average of saturated conductivities
const RangeU cond_sat_i = conductivity_f_i(1.0);
const RangeU cond_sat_o = conductivity_f_o(1.0);
// Upwinding conductivity factor
RangeU cond_upwind_factor(0.0);
if (dudn < 0.0) {
// inward flux
cond_upwind_factor = cond_o / cond_sat_o;
}
else {
// outward flux
cond_upwind_factor = cond_i / cond_sat_i;
}
cond_f = 2.0 * cond_upwind_factor / ( 1.0/(cond_sat_i + 1E-30)
+ 1.0/(cond_sat_o + 1E-30) );
}
else // RichardsFVUpwinding::fullUpwind
{
// No average: Use upwind conductivity
if (dudn < 0.0) {
// inward flux
cond_f = cond_o;
}
else {
// outward flux
cond_f = cond_i;
}
}
// Water flux in normal direction w.r.t the intersection
auto water_flux_n = - cond_f*dudn;
// // Upwinding
// if (water_flux_n > 0.)
// water_flux_n = - cond_i*dudn;
// else
// water_flux_n = - cond_o*dudn;
const auto water_flux_n = - cond_f * dudn;
// Symmetric contribution to residual on inside element
r_i.accumulate(lfsv_i, 0, water_flux_n*volume_f );
......@@ -290,6 +325,9 @@ private:
const std::shared_ptr<const Parameter> _param;
/// Class providing boundary conditions
const std::shared_ptr<const Boundary> _boundary;
/// Upwinding setting
const RichardsUpwinding _upwinding;
/// Operator internal time
double _time;
};
......
#ifndef DUNE_DORIE_RICHARDS_LOP_CFG_HH
#define DUNE_DORIE_RICHARDS_LOP_CFG_HH
#include <tuple>
#include <string>
#include <dune/common/exceptions.hh>
#include <dune/common/parametertree.hh>
#include <dune/dorie/common/logging.hh>
namespace Dune {
namespace Dorie {
namespace Operator {
/**
* \addtogroup RichardsModel
* \{
*/
//! Interior penalty type
enum class RichardsDGMethod
{
NIPG, //!< Non-symmetric interior penalty
SIPG, //!< Symmetric interior penalty (default)
OBB, //!< Oden, Babuska, Baumann (no penalty term)
IIP //!< Incomplete interior penalty (no symmetry term)
};
//! Gradient weighting
enum class RichardsDGWeights
{
weightsOn, //!< harmonic weighting of flux contributions
weightsOff //!< no weighting of flux contributions
};
//! Upwinding type
enum class RichardsUpwinding
{
fullUpwind, //!< Upwind conductivity
semiUpwind, //!< Upwind conductivity factor
none //!< No upwinding
};
/// Read operator settings from a parameter file
/* \return Tuple of RichardsUpwinding, RichardsDGMethod, and
* RichardsDGWeights to be inserted into the local operator constructor.
*/
inline auto read_operator_settings(
const Dune::ParameterTree& inifile)
-> std::tuple<Operator::RichardsUpwinding,
Operator::RichardsDGMethod,
Operator::RichardsDGWeights>
{
const auto log = Dorie::get_logger(log_richards);
log->debug("Reading local operator settings:");
RichardsUpwinding upwinding;
RichardsDGMethod method = RichardsDGMethod::SIPG;
RichardsDGWeights weights = RichardsDGWeights::weightsOn;
// Upwinding
const auto upwinding_str
= inifile.get<std::string>("numerics.upwinding");
log->debug(" Upwinding scheme: {}", upwinding_str);
if (upwinding_str == "none")
upwinding = RichardsUpwinding::none;
else if (upwinding_str == "semiUpwind")
upwinding = RichardsUpwinding::semiUpwind;
else if (upwinding_str == "fullUpwind")
upwinding = RichardsUpwinding::fullUpwind;
else {
log->error("Unknown upwinding scheme: {}", upwinding_str);
DUNE_THROW(Dune::IOError, "Unknown upwinding scheme");
}
// Return here if the local operator is FV only
if (inifile.get<int>("numerics.FEorder") == 0)
{
log->debug(" Ignoring settings 'numerics.DGMethod' and "
"'numerics.DGWeights' for finite volume solver.");
return std::make_tuple(upwinding, method, weights);
}
// DG Method
const auto method_str = inifile.get<std::string>("numerics.DGMethod");
log->debug(" DG method: {}", method_str);
if (method_str == "SIPG")
method = RichardsDGMethod::SIPG;
else if (method_str == "NIPG")
method = RichardsDGMethod::NIPG;
else if (method_str == "OBB")
method = RichardsDGMethod::OBB;
else if (method_str == "IIP")
method = RichardsDGMethod::IIP;
else {
log->error("Unknown DG method: {}", method_str);
DUNE_THROW(Dune::IOError, "Unknown DG method");
}
// Harmonic weighting of fluxes
const auto weights_str = inifile.get<bool>("numerics.DGWeights");
log->debug(" Harmonic weights for interface fluxes (DG): {}",
weights_str);
if (weights_str)
weights = RichardsDGWeights::weightsOn;
else
weights = RichardsDGWeights::weightsOff;
return std::make_tuple(upwinding, method, weights);
}
/**
* \}
* // group RichardsModel
*/
} // namespace Operator
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_RICHARDS_LOP_CFG_HH
......@@ -55,30 +55,27 @@ ModelRichards<Traits>::ModelRichards (
fsource = std::make_shared<const FlowSource>(inifile);
finitial = FlowInitialFactory::create(inifile, gv, this->_log);
// Read local operator settings
namespace OP = Dune::Dorie::Operator;
const auto settings = OP::read_operator_settings(inifile);
// --- Local Operators ---
if constexpr (order>0)
{
if constexpr (order>0)
{
this->_log->debug("Setting up local grid operators: DG method");
#ifdef EXPERIMENTAL_DG_FEATURES
// read experimental settings from inifile
namespace OP = Dune::Dorie::Operator;
const auto settings = OP::read_experimental_operator_settings(inifile);
const auto method = std::get<OP::RichardsDGMethod::Type>(settings);
const auto upwinding = std::get<OP::RichardsDGUpwinding::Type>(settings);
const auto weights = std::get<OP::RichardsDGWeights::Type>(settings);
const auto upwinding = std::get<OP::RichardsUpwinding>(settings);
const auto method = std::get<OP::RichardsDGMethod>(settings);
const auto weights = std::get<OP::RichardsDGWeights>(settings);
slop = std::make_unique<SLOP>(inifile, fparam, fboundary, fsource,
method, upwinding, weights);
#else
slop = std::make_unique<SLOP>(inifile, fparam, fboundary, fsource);
#endif // EXPERIMENTAL_DG_FEATURES
tlop = std::make_unique<TLOP>(inifile, fparam);
}
else {
this->_log->debug("Setting up local grid operators: FV method");
slop = std::make_unique<SLOP>(fparam, fboundary);
const auto upwinding = std::get<OP::RichardsUpwinding>(settings);
slop = std::make_unique<SLOP>(fparam, fboundary, upwinding);
tlop = std::make_unique<TLOP>(fparam);
}
......
......@@ -33,6 +33,7 @@
#include <dune/dorie/model/richards/flow_source.hh>
#include <dune/dorie/model/richards/local_operator_DG.hh>
#include <dune/dorie/model/richards/local_operator_FV.hh>
#include <dune/dorie/model/richards/local_operator_cfg.hh>
namespace Dune{
namespace Dorie{
......
......@@ -5,7 +5,10 @@ _asset_path = "${CMAKE_CURRENT_LIST_DIR}"
_evaluation = ode
_order = 0, 1, 2, 3 | expand prec
__name = ode_homogeneous_sand_k{_order}
_upwind = none, semiUpwind, fullUpwind | expand
__name = ode_homogeneous_sand_k{_order}_{_upwind}
{_order} > 0 and {_upwind} == fullUpwind | exclude
[grid]
gridType = rectangular
......@@ -34,6 +37,7 @@ NewtonParameters.Reduction = 1E-10
numerics.penaltyFactor = 10
numerics.FEorder = 0, 1, 2, 3 | expand prec
numerics.upwinding = {_upwind}
[_ode]
......@@ -41,4 +45,4 @@ headLower = 0.
flux = -5.55e-6
head_abstol = 3E-4, 3E-6, 2E-8, 6E-8 | expand prec
flux_abstol = 1E100, 5E-11, 1E-9, 3E-12 | expand prec
flux_rt_abstol = 4E-14
\ No newline at end of file
flux_rt_abstol = 4E-14
......@@ -5,7 +5,10 @@ _asset_path = "${CMAKE_CURRENT_LIST_DIR}"
_evaluation = ode
_order = 0, 1, 2, 3 | expand prec
__name = ode_homogeneous_silt_k{_order}
_upwind = none, semiUpwind, fullUpwind | expand
__name = ode_homogeneous_silt_k{_order}_{_upwind}
{_order} > 0 and {_upwind} == fullUpwind | exclude
[grid]
gridType = rectangular
......@@ -38,6 +41,7 @@ NewtonParameters.Reduction = 1E-10
numerics.penaltyFactor = 10
numerics.FEorder = 0, 1, 2, 3 | expand prec
numerics.upwinding = {_upwind}
[_ode]
......
......@@ -5,7 +5,10 @@ _asset_path = "${CMAKE_CURRENT_LIST_DIR}"
_evaluation = ode
_order = 0, 1, 2, 3 | expand prec
__name = ode_layered_k{_order}
_upwind = none, semiUpwind, fullUpwind | expand
__name = ode_layered_k{_order}_{_upwind}
{_order} > 0 and {_upwind} != none | exclude
[grid]
gridType = rectangular
......@@ -38,6 +41,7 @@ NewtonParameters.Reduction = 1E-10
numerics.penaltyFactor = 10
numerics.FEorder = {_order}
numerics.upwinding = {_upwind}
[_ode]
......@@ -45,4 +49,4 @@ headLower = 0.
flux = -5.55e-6
head_abstol = 1E-3, 3E-5, 4E-8, 4E-6 | expand prec
flux_abstol = 1E100, 5E-10, 4E-9, 5E-11 | expand prec
flux_rt_abstol = 4E-14
\ No newline at end of file
flux_rt_abstol = 4E-14