Commit ac23c2d0 authored by Lukas Riedel's avatar Lukas Riedel

Use relative conductivity factor in Dirichlet flux term

Multiply relative conductivity to numeric flux now.
Add it also to the evaporation estimation.
parent 290fdd5a
......@@ -474,10 +474,22 @@ public:
// compute numerical flux estimation
const RF numFlux = satCond_s * ( - (gradu_s * normal) + penalty * jump);
RF head_upwind;
if (numFlux > 0)
head_upwind = param.headMillerToRef(u_s, p_global_s);
else
head_upwind = g; // Miller here?!
// compute saturation from matrix head
const RF saturation_upwind = param.saturation(head_upwind, p_global_s);
// compute relative conductivity
const RF relCond_upwind = param.condFactor(saturation_upwind, p_global_s);
// switches between Neumann and Dirichlet bc
// this choice is kept constant for one time step
normal_flux = boundary.j(ig.intersection(),it.position(),time);
if(numFlux < normal_flux) {
if(relCond_upwind * numFlux < normal_flux) {
bcType = BCType::dirichlet;
bc_type_cache.emplace(it_id,BoundaryCondition::Dirichlet);
}
......@@ -557,7 +569,7 @@ public:
// consistency term
// + penalty term
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, numFlux * phi_s[i] * factor);
r_s.accumulate(lfsv_s,i, relCond_upwind * numFlux * phi_s[i] * factor);
// update residual (symmetry term)
// (non-)symmetric IP term
......
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