Commit 7ae6d3ef authored by Lukas Riedel's avatar Lukas Riedel

Merge branch 'debug/improve-evaporation-BC' into 'master'

Improve evaporation BC switch

I improved the evaporation BC by changing the condition for which the applied BC switches from Neumann to Dirichlet BC: For any quadrature point with evaporation BC, the numeric flux is calculated once per time step. If this numeric flux is smaller than the flux applied by the Neumann BC, the Dirichlet BC is applied. This smoothes the transition between the two BCs.

Tested improvements:
* The head value given in the evaporation BC is not a mere cutoff value but the target boundary condition if the Neumann flux cannot be sustained
* Evaporation from homogeneous surfaces is very stable and reaches high time step values
* The new condition allows for a smooth transition back from Dirichlet to Neumann if the soil regains water saturation (e.g., by a rising water table)

* Evaporation from a Miller regime is still not very stable. This is due to the fact that the Dirichlet BC itself is not very stable for a Miller regime.
* Higher FEM orders can reduce the numeric stability
* Numeric flux is not a physical flux and depends on the `dg.penaltyFactor`.

* The code for the switch check is copypasta from the Dirichlet BC.

See merge request !8
parents 870121ae 4a88ce1a
......@@ -99,7 +99,8 @@ adding an empty line, make text **bold** or ``monospaced``.
<parameter name="cells">
<definition> Initial number of cells in each dimension (x, then y, then z) if ``gridType``
is set to ``rectangular``. </definition>
is set to ``rectangular``. This represents the coarsest level of the grid,
- i.e. refinement level 0. </definition>
<values> int &times; int (&times; int) </values>
<suggestion> 100 100 </suggestion>
......@@ -175,13 +176,13 @@ adding an empty line, make text **bold** or ``monospaced``.
<definition> Minimum time step that is allowed before DORiE stops running.
with an error, in seconds. </definition>
<values> float </values>
<suggestion> 10 </suggestion>
<suggestion> 0.1 </suggestion>
<parameter name="startTimestep">
<definition> Value of the first time step in seconds. </definition>
<values> float </values>
<suggestion> 100 </suggestion>
<suggestion> 10 </suggestion>
<parameter name="maxTimestep">
......@@ -311,9 +312,10 @@ adding an empty line, make text **bold** or ``monospaced``.
<definition> The maximum refinement level kept in the grid. This is a useful
tool to prevent over-refinement. If this value is high, the grid can
be refined to an arbitrary degree, leading to an evenly distributed error
across the grid. </definition>
across the grid. Make sure you avoid refinement levels which imply grid cell sizes
beyond the Richards continuum scale. </definition>
<values> int </values>
<suggestion> 15 </suggestion>
<suggestion> 10 </suggestion>
<parameter name="minLevel">
......@@ -42,7 +42,7 @@ Evaporation
.. object:: evaporation <flux> <head>
Set a volumetric flux [m/s] at the boundary (Neumann BC). If this flux cannot be sustained, set a fixed matric head at the boundary instead (Dirichlet BC).
Set a volumetric flux [m/s] (Neumann BC) and a fallback matric head [m] (Dirichlet BC) at the boundary. The Neumann BC is invoked as long as the water content inside the soil can sustain the flux. Before every time step, the numeric flux caused by the given Dirichlet BC is estimated. If it is lower than the given Neumann flux, the Dirichlet BC is applied.
Datafile Structure
......@@ -433,18 +433,52 @@ public:
else if (BoundaryCondition::isEvaporation(bc))
// estimate the dirichlet flux
// evaluate gradient of basis functions (we assume Galerkin method lfsu = lfsv)
const auto& gradphi_s = cache[order_s].evaluateJacobian(local_s,lfsu_s.finiteElement().localBasis());
// transform gradients to real element
const auto jac = ig.inside().geometry().jacobianInverseTransposed(local_s);
std::vector<Vector> tgradphi_s(lfsu_s.size());
for (unsigned int i = 0; i<lfsu_s.size(); i++)[i][0],tgradphi_s[i]);
// compute gradient of u
Vector gradu_s(0.);
for (unsigned int i = 0; i<lfsu_s.size(); i++)
// update with gravity vector
gradu_s -= param.gravity();
// face normal vector
const Domain normal = ig.unitOuterNormal(it.position());
// jump relative to Dirichlet value
const RF g = boundary.g(ig.intersection(),it.position(),time);
const RF jump = u_s - g;
const Domain xGlobal_s = ig.inside().geometry().global(local_s);
// conductivity at saturation
const RF satCond_s = param.condRefToMiller(param.K(xGlobal_s), xGlobal_s);
const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);
// compute numerical flux estimation
const RF numFlux = satCond_s * ( - (gradu_s * normal) + penalty * jump);
// switches between Neumann and Dirichlet bc
// this choice is kept constant for one time step
RF u_0 = boundary.g(ig.intersection(),it.position(),time);
if(u_s > u_0) {
normal_flux = boundary.j(ig.intersection(),it.position(),time);
bcType = BCType::vonNeumann;
else {
normal_flux = boundary.j(ig.intersection(),it.position(),time);
if(numFlux < normal_flux) {
bcType = BCType::dirichlet;
else {
bcType = BCType::vonNeumann;
else if (BoundaryCondition::isOther(bc))
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment