...
 
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;
}
......
......@@ -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));
// Interface conductivity factor (Upwinding applies)
RangeU cond_f(0.0);
if (_upwinding == RichardsUpwinding::none)
{
// Harmonic average for interface conductivity
const RangeU cond_f = 2.0/( 1.0/(cond_i + 1E-30)
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;
}
// 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;
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)
{
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]
......
......@@ -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]
......