Commit c15f923b authored by Santiago Ospina's avatar Santiago Ospina
Browse files

[Flux Reconstruction] Fix error on lifting

* DOF of flux reconstruction and lifting were combined in a wrong way.
* Discrete gradient sign was removed because it was wrong. The sign has been always part of the dG scheme.
parent 9d5a88f8
......@@ -483,6 +483,8 @@ public:
}
if constexpr (extended) {
rl.assign(lfsw_s_cache.size(),0.0);
Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaSkeleton>::
alpha_skeleton(lop,ig,
lfsw_s_cache.localFunctionSpace(),pl,lfsv_volume,
......@@ -577,6 +579,7 @@ public:
}
if constexpr (extended) {
rl.assign(lfsw_s_cache.size(),0.0);
Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaBoundary>::
alpha_boundary(lop,ig,lfsw_s_cache.localFunctionSpace(),pl,lfsv_volume,rl_view);
......
......@@ -484,9 +484,6 @@ public:
// (we assume non-Galerkin method lfsu != lfsv)
lfsv_s.finiteElement().localBasis().evaluateFunction(sub_p_local_s,phiv_s);
lfsv_n.finiteElement().localBasis().evaluateFunction(sub_p_local_n,phiv_n);
// discrete gradient sign.
dg_sign = -1.;
} else if constexpr (lopcase == LOPCase::RTVolume) {
// (we assume lfsv = RTk)
lfsv_s.finiteElement().localBasis().evaluateFunction(p_local_s,gradphiv_s);
......@@ -740,7 +737,6 @@ public:
if constexpr (lopcase == LOPCase::RTSkeleton) {
// (we assume non-Galerkin method lfsu != lfsv)
lfsv_s.finiteElement().localBasis().evaluateFunction(it.position(),phiv_s);
dg_sign = -1.;
} else if constexpr (lopcase == LOPCase::RTVolume) {
// (we assume lfsv = RTk)
lfsv_s.finiteElement().localBasis().evaluateFunction(p_local_s,gradphiv_s);
......
......@@ -241,11 +241,11 @@ void RichardsSimulation<Traits>::step()
try {
// compute lifting
slop->penalty_term(false);
slop->consistency_term(false);
slop->symmetry_term(false);
slop->diffusion_term(false);
lifting_gf->update(*u);
slop->penalty_term(true);
slop->consistency_term(true);
slop->symmetry_term(true);
slop->diffusion_term(true);
} catch (Dune::Exception &e) {
this->_log->debug(" Local lifting failed");
......
......@@ -285,6 +285,10 @@ public:
acc += deviation;
acc_square += deviation * deviation;
if (flux_rec_gf->check_jumps(1e-10))
DUNE_THROW(Dune::MathError,
"Jumps in flux reconstruction are not admissible!");
std::cout << "wc_old: " << wc_old << std::endl;
std::cout << "wc_new: " << wc_new << std::endl;
std::cout << "integrated flux: " << flux_int << std::endl;
......
Supports Markdown
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