Skip to content

GitLab

D
dorie
Commit 052a90b3 authored by Lukas Riedel's avatar Lukas Riedel
Browse files
Options

Merge branch 'cleanup/dg-operator' into cleanup/dg-operator-check-residuals

parents c2ff5c7a 969acd6b
Hide whitespace changes
Inline Side-by-side
Showing with 47 additions and 38 deletions
dune/dorie/solver/operator_DG.hh
View file @ 052a90b3
......@@ -162,8 +162,12 @@ public:
// loop over quadrature points
for (const auto& it : quadratureRule(gt,intorder))
{
// evaluate position in element local and global coordinates
const auto local = it.position();
const auto global = gt.global(local);
// evaluate basis functions
auto& phi = cache[order].evaluateFunction(it.position(),lfsu.finiteElement().localBasis());
auto& phi = cache[order].evaluateFunction(local,lfsu.finiteElement().localBasis());
// evaluate u
RF u = 0.;
......@@ -171,10 +175,10 @@ public:
u += x(lfsu,i) * phi[i];
// evaluate gradient of basis functions (we assume Galerkin method lfsu = lfsv)
const auto& js = cache[order].evaluateJacobian(it.position(),lfsu.finiteElement().localBasis());
const auto& js = cache[order].evaluateJacobian(local,lfsu.finiteElement().localBasis());
// transform gradients to real element
const auto jac = gt.jacobianInverseTransposed(it.position());
const auto jac = gt.jacobianInverseTransposed(local);
std::vector<Vector> gradphi(lfsu.size());
for (unsigned int i = 0; i<lfsu.size(); i++)
jac.mv(js[i][0],gradphi[i]);
......@@ -186,20 +190,17 @@ public:
gradu -= param.gravity();
const Domain xGlobal = eg.entity().geometry().global(it.position());
const RF head = param.headMillerToRef(u,xGlobal);
// compute saturation from matrix head
const RF saturation = param.saturation(head,xGlobal);
const RF head = param.headMillerToRef(u,global);
const RF saturation = param.saturation(head,global);
// compute conductivity (with regard to saturation)
const RF relCond = param.condFactor(saturation,xGlobal);
const RF satCond = param.K(xGlobal);
const RF cond = param.condRefToMiller(relCond*satCond, xGlobal);
const RF relCond = param.condFactor(saturation,global);
const RF satCond = param.K(global);
const RF cond = param.condRefToMiller(relCond*satCond, global);
// integration factor
const RF factor = it.weight() * gt.integrationElement(it.position());
const RF factor = it.weight() * gt.integrationElement(local);
// update residual
for (unsigned int i = 0; i<lfsv.size(); i++)
......@@ -254,6 +255,17 @@ public:
const auto local_s = ig.geometryInInside().global(it.position());
const auto local_n = ig.geometryInOutside().global(it.position());
// face normal vector
const Domain normal = ig.unitOuterNormal(it.position());
// position of quadrature point in global coordinates
// add offset for parameter queries
const RF eps = 1e-9;
auto global_s = ig.inside().geometry().global(local_s);
auto global_n = ig.outside().geometry().global(local_n);
global_s = global_s.axpy(-eps,normal);
global_n = global_n.axpy(eps,normal);
// evaluate basis functions
const auto& phi_s = cache[order_s].evaluateFunction(local_s,lfsu_s.finiteElement().localBasis());
const auto& phi_n = cache[order_n].evaluateFunction(local_n,lfsu_n.finiteElement().localBasis());
......@@ -291,18 +303,12 @@ public:
gradu_s -= param.gravity();
gradu_n -= param.gravity();
// face normal vector
const Domain normal = ig.unitOuterNormal(it.position());
// compute jump in solution
const RF jump = u_s - u_n;
const Domain xGlobal_s = ig.inside().geometry().global(local_s);
const Domain xGlobal_n = ig.outside().geometry().global(local_n);
// conductivity at saturation
const RF satCond_s = param.condRefToMiller(param.K(xGlobal_s), xGlobal_s);
const RF satCond_n = param.condRefToMiller(param.K(xGlobal_n), xGlobal_n);
const RF satCond_s = param.condRefToMiller(param.K(global_s), global_s);
const RF satCond_n = param.condRefToMiller(param.K(global_n), global_n);
// compute weights
RF omega_s = 0.5, omega_n = 0.5;
......@@ -324,15 +330,15 @@ public:
RF head_upwind;
if (numFlux > 0)
head_upwind = param.headMillerToRef(u_s, xGlobal_s);
head_upwind = param.headMillerToRef(u_s, global_s);
else
head_upwind = param.headMillerToRef(u_n, xGlobal_n);
head_upwind = param.headMillerToRef(u_n, global_n);
// compute saturation from matrix head
const RF saturation_upwind = param.saturation(head_upwind, xGlobal_s);
const RF saturation_upwind = param.saturation(head_upwind, global_s);
// compute relative conductivity
const RF relCond_upwind = param.condFactor(saturation_upwind, xGlobal_s);
const RF relCond_upwind = param.condFactor(saturation_upwind, global_s);
// update residual (flux term)
// diffusion term
......@@ -612,15 +618,17 @@ public:
// loop over quadrature points
for (const auto& it : quadratureRule(gt,intorder))
{
// evaluate position in element local coordinates
const auto local = it.position();
// evaluate values of shape functions (we assume Galerkin method lfsu = lfsv)
std::vector<Scalar> phi(lfsv.size());
lfsv.finiteElement().localBasis().evaluateFunction(it.position(),phi);
lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
// scaled weight factor
const RF factor = it.weight() * eg.geometry().integrationElement(it.position());
const RF factor = it.weight() * gt.integrationElement(local);
const RF source = sourceTerm.q(eg.entity(),it.position(),time) / eg.geometry().volume();
const RF source = sourceTerm.q(eg.entity(),local,time) / gt.volume();
// update residual
for (unsigned int i = 0; i<lfsv.size(); i++)
......@@ -733,26 +741,27 @@ template<typename Traits, typename Parameter, typename FEM, bool adjoint>
// loop over quadrature points
for (const auto& it : quadratureRule(gt,intorder))
{
// evaluate position in element local and global coordinates
const auto local = it.position();
const auto global = gt.global(local);
// evaluate basis functions
std::vector<Scalar> phi(lfsu.size());
lfsu.finiteElement().localBasis().evaluateFunction(it.position(),phi);
lfsu.finiteElement().localBasis().evaluateFunction(local,phi);
// evaluate u
RF u = 0.;
for (unsigned int i = 0; i<lfsu.size(); i++)
u += x(lfsu,i) * phi[i];
u += x(lfsu,i) * phi[i];
// integration factor
const RF factor = it.weight() * eg.geometry().integrationElement(it.position());
const Domain xGlobal = eg.entity().geometry().global(it.position());
// calculate water content from matric head
const RF head = param.headMillerToRef(u,global);
const RF water_content = param.waterContent(head,global);
const RF head = param.headMillerToRef(u,xGlobal);
// integration factor
const RF factor = it.weight() * gt.integrationElement(local);
// update residual
const RF water_content = param.waterContent(head,xGlobal);
//std::cout << "water cont " << water_content << std::endl;
for (unsigned int i = 0; i<lfsv.size(); i++)
r.accumulate(lfsv,i, water_content * phi[i] * factor);
......
dune/dorie/solver/util_interpolator.hh
View file @ 052a90b3
......@@ -34,7 +34,7 @@ namespace Dune {
const InterpolationMethod method;
explicit InterpolatorBase(InterpolationMethod method_, const ParameterTree& config_, const Dune::MPIHelper& helper_, const int verbose_):
config(config_), helper(helper_), verbose(verbose_), eps(1e-8), initialized(false), method(method_)
config(config_), helper(helper_), verbose(verbose_), eps(1e-9), initialized(false), method(method_)
{}
RF evaluate(const Array& field, const Domain& pos)
......
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