Commit 0f598cd0 authored by Dennis Gläser's avatar Dennis Gläser
Browse files

[mpfa] do not do neumann flux scaling

For non-zero neumann boundary conditions we did a scaling to recover
DgradU at the boundary from the prescribed flux. We don't do this anymore
and assume that when the UseTpfaBoundary property is set to false, that the
user provides -DgradU at the boundary instead of, e.g., mass fluxes.
parent 20f0c059
...@@ -66,6 +66,7 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa> ...@@ -66,6 +66,7 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using IndexType = typename GridView::IndexSet::IndexType; using IndexType = typename GridView::IndexSet::IndexType;
using Stencil = std::vector<IndexType>; using Stencil = std::vector<IndexType>;
static const bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
static const MpfaMethods method = GET_PROP_VALUE(TypeTag, MpfaMethod); static const MpfaMethods method = GET_PROP_VALUE(TypeTag, MpfaMethod);
public: public:
...@@ -108,7 +109,10 @@ public: ...@@ -108,7 +109,10 @@ public:
flux += tij[localIdx++]*h; flux += tij[localIdx++]*h;
} }
return flux + fluxVarsCache.advectionNeumannFlux(phaseIdx); if (useTpfaBoundary)
return flux;
else
return flux + fluxVarsCache.advectionNeumannFlux(phaseIdx);
} }
static Stencil stencil(const Problem& problem, static Stencil stencil(const Problem& problem,
......
...@@ -203,8 +203,7 @@ public: ...@@ -203,8 +203,7 @@ public:
T_ += D; T_ += D;
} }
template<typename UpwindFactorFunction> void assembleNeumannFluxes(const unsigned int eqIdx)
void assembleNeumannFluxes(const UpwindFactorFunction& upwindFactor, const unsigned int eqIdx)
{ {
if (!onBoundary() || GET_PROP_VALUE(TypeTag, UseTpfaBoundary)) if (!onBoundary() || GET_PROP_VALUE(TypeTag, UseTpfaBoundary))
return; return;
...@@ -221,11 +220,7 @@ public: ...@@ -221,11 +220,7 @@ public:
auto neumannFlux = problem_().neumann(element, this->fvGeometry_(), this->elemVolVars_(), globalScvf)[eqIdx]; auto neumannFlux = problem_().neumann(element, this->fvGeometry_(), this->elemVolVars_(), globalScvf)[eqIdx];
neumannFlux *= globalScvf.area(); neumannFlux *= globalScvf.area();
// recover -k*gradh // The flux is assumed to be prescribed in the form of -D*gradU
const auto& insideScv = fvGeometry_().scv(globalScvf.insideScvIdx());
const auto& volVars = elemVolVars_()[insideScv];
neumannFlux /= upwindFactor(volVars);
neumannFluxes_[fluxFaceIdx] = neumannFlux; neumannFluxes_[fluxFaceIdx] = neumannFlux;
} }
......
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