Commit 99db5435 authored by Lukas Riedel's avatar Lukas Riedel 📝

Merge branch 'bugix/conductivity-in-dirichlet-bc' into feature/local-operator-schemes

parents 6cdbdb13 237faaa9
...@@ -543,10 +543,22 @@ public: ...@@ -543,10 +543,22 @@ public:
// compute numerical flux estimation // compute numerical flux estimation
const RF numFlux = satCond_s * ( - (gradu_s * normal) + penalty * jump); 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 // switches between Neumann and Dirichlet bc
// this choice is kept constant for one time step // this choice is kept constant for one time step
normal_flux = boundary.j(ig.intersection(),it.position(),time); normal_flux = boundary.j(ig.intersection(),it.position(),time);
if(numFlux < normal_flux) { if(relCond_upwind * numFlux < normal_flux) {
bcType = BCType::dirichlet; bcType = BCType::dirichlet;
bc_type_cache.emplace(it_id,BoundaryCondition::Dirichlet); bc_type_cache.emplace(it_id,BoundaryCondition::Dirichlet);
} }
...@@ -626,7 +638,7 @@ public: ...@@ -626,7 +638,7 @@ public:
// consistency term // consistency term
// + penalty term // + penalty term
for (unsigned int i = 0; i<lfsv_s.size(); i++) 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) // update residual (symmetry term)
// (non-)symmetric IP 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