Commit f723a504 authored by Lukas Riedel's avatar Lukas Riedel

guess evaporation switch by evaluating numeric flux from Dirichlet BC

parent 870121ae
......@@ -436,14 +436,51 @@ public:
// 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;
bc_type_cache.emplace(it_id,BoundaryCondition::Neumann);
}
else {
normal_flux = boundary.j(ig.intersection(),it.position(),time);
// 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++)
jac.mv(gradphi_s[i][0],tgradphi_s[i]);
// compute gradient of u
Vector gradu_s(0.);
for (unsigned int i = 0; i<lfsu_s.size(); i++)
gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[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);
if(numFlux < normal_flux) {
bcType = BCType::dirichlet;
bc_type_cache.emplace(it_id,BoundaryCondition::Dirichlet);
std::cout << "- Falling back to Dirichlet BC." << std::endl;
}
else {
bcType = BCType::vonNeumann;
bc_type_cache.emplace(it_id,BoundaryCondition::Neumann);
}
}
......
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