From ed887c6a9a4bdaf4928ace80264a5e9651472648 Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser" <dennis.glaeser@iws.uni-stuttgart.de>
Date: Sat, 9 Feb 2019 14:48:19 +0100
Subject: [PATCH] [facet][box][darcyslaw] use tpfa flux also on interior
Dirichlet BCS
---
dumux/multidomain/facet/box/darcyslaw.hh | 91 ++++++-------------
.../facet/box/fluxvariablescache.hh | 10 +-
2 files changed, 34 insertions(+), 67 deletions(-)
diff --git a/dumux/multidomain/facet/box/darcyslaw.hh b/dumux/multidomain/facet/box/darcyslaw.hh
index 29de16a088..c204533701 100644
--- a/dumux/multidomain/facet/box/darcyslaw.hh
+++ b/dumux/multidomain/facet/box/darcyslaw.hh
@@ -83,82 +83,49 @@ public:
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
const auto& insideVolVars = elemVolVars[insideScv];
- // evaluate user-defined interior boundary types
- const auto bcTypes = problem.interiorBoundaryTypes(element, scvf);
-
- static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
+ // evaluate the flux in a tpfa manner
+ const auto& supportPointShapeValues = fluxVarCache.shapeValuesAtTpfaSupportPoint();
- // on interior Neumann boundaries, evaluate the flux using the facet permeability
- if (bcTypes.hasOnlyNeumann())
+ Scalar rho = 0.0;
+ Scalar supportPressure = 0.0;
+ for (const auto& scv : scvs(fvGeometry))
{
- const auto& supportPointShapeValues = fluxVarCache.shapeValuesAtTpfaSupportPoint();
+ const auto& volVars = elemVolVars[scv];
+ rho += volVars.density(phaseIdx)*shapeValues[scv.indexInElement()][0];
+ supportPressure += volVars.pressure(phaseIdx)*supportPointShapeValues[scv.indexInElement()][0];
+ }
- Scalar rho = 0.0;
- Scalar supportPressure = 0.0;
- for (const auto& scv : scvs(fvGeometry))
- {
- const auto& volVars = elemVolVars[scv];
- rho += volVars.density(phaseIdx)*shapeValues[scv.indexInElement()][0];
- supportPressure += volVars.pressure(phaseIdx)*supportPointShapeValues[scv.indexInElement()][0];
- }
+ // the transmissibility on the matrix side
+ const auto dm = (scvf.ipGlobal() - fluxVarCache.tpfaSupportPoint()).two_norm();
+ const auto tm = 1.0/dm*vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), scvf.unitOuterNormal());
- // compute tpfa flux such that flux and pressure continuity holds
- const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf);
+ // compute flux depending on the user's choice of boundary types
+ const auto bcTypes = problem.interiorBoundaryTypes(element, scvf);
+ const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf);
+ // compute tpfa flux such that flux continuity holds with flux in fracture
+ Scalar flux;
+ if (bcTypes.hasOnlyNeumann())
+ {
// On surface grids, use sqrt of aperture as distance measur
using std::sqrt;
- const auto df = dim == dimWorld ? 0.5*facetVolVars.extrusionFactor()
- : 0.5*sqrt(facetVolVars.extrusionFactor());
- const auto dm = (scvf.ipGlobal() - fluxVarCache.tpfaSupportPoint()).two_norm();
-
- const auto tm = 1.0/dm*vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), scvf.unitOuterNormal());
+ const auto df = dim == dimWorld ? 0.5*facetVolVars.extrusionFactor() : 0.5*sqrt(facetVolVars.extrusionFactor());
const auto tf = 1.0/df*vtmv(scvf.unitOuterNormal(), facetVolVars.permeability(), scvf.unitOuterNormal());
- auto flux = tm*tf/(tm+tf)*(supportPressure - facetVolVars.pressure(phaseIdx))
- *scvf.area()*insideVolVars.extrusionFactor();
-
- if (enableGravity)
- flux -= rho*scvf.area()*insideVolVars.extrusionFactor()
- *vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), problem.gravityAtPos(scvf.center()));
-
- return flux;
+ flux = tm*tf/(tm+tf)*(supportPressure - facetVolVars.pressure(phaseIdx))*scvf.area()*insideVolVars.extrusionFactor();
}
-
- // on interior Dirichlet boundaries use the facet pressure and evaluate flux
else if (bcTypes.hasOnlyDirichlet())
- {
- // create vector with nodal pressures
- std::vector<Scalar> pressures(element.subEntities(dim));
- for (const auto& scv : scvs(fvGeometry))
- pressures[scv.localDofIndex()] = elemVolVars[scv].pressure(phaseIdx);
-
- // substitute with facet pressures for those scvs touching this facet
- for (const auto& scvfJ : scvfs(fvGeometry))
- if (scvfJ.interiorBoundary() && scvfJ.facetIndexInElement() == scvf.facetIndexInElement())
- pressures[ fvGeometry.scv(scvfJ.insideScvIdx()).localDofIndex() ]
- = problem.couplingManager().getLowDimVolVars(element, scvfJ).pressure(phaseIdx);
-
- // evaluate gradP - rho*g at integration point
- Scalar rho(0.0);
- Dune::FieldVector<Scalar, dimWorld> gradP(0.0);
- for (const auto& scv : scvs(fvGeometry))
- {
- rho += elemVolVars[scv].density(phaseIdx)*shapeValues[scv.indexInElement()][0];
- gradP.axpy(pressures[scv.localDofIndex()], fluxVarCache.gradN(scv.indexInElement()));
- }
-
- if (enableGravity)
- gradP.axpy(-rho, problem.gravityAtPos(scvf.center()));
-
- // apply matrix permeability and return the flux
- return -1.0*scvf.area()
- *insideVolVars.extrusionFactor()
- *vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), gradP);
- }
-
+ flux = tm*(supportPressure - facetVolVars.pressure(phaseIdx))*scvf.area()*insideVolVars.extrusionFactor();
// mixed boundary types are not supported
else
DUNE_THROW(Dune::NotImplemented, "Mixed boundary types are not supported");
+
+ static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
+ if (enableGravity)
+ flux -= rho*scvf.area()*insideVolVars.extrusionFactor()
+ *vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), problem.gravityAtPos(scvf.center()));
+
+ return flux;
}
// compute transmissibilities ti for analytical jacobians
diff --git a/dumux/multidomain/facet/box/fluxvariablescache.hh b/dumux/multidomain/facet/box/fluxvariablescache.hh
index 41593ec3ef..43cc613541 100644
--- a/dumux/multidomain/facet/box/fluxvariablescache.hh
+++ b/dumux/multidomain/facet/box/fluxvariablescache.hh
@@ -69,9 +69,9 @@ public:
// on interior boundaries with Neumann BCs, prepare the shape values at a point
// inside the element whose orthogonal projection is the integration point on scvf
- if (scvf.interiorBoundary() && problem.interiorBoundaryTypes(element, scvf).hasOnlyNeumann())
+ if (scvf.interiorBoundary())
{
- isInteriorNeumannCache_ = true;
+ isInteriorBoundaryCache_ = true;
const auto& geometry = element.geometry();
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
@@ -96,14 +96,14 @@ public:
//! returns the integration point inside the element for interior boundaries
const GlobalPosition& tpfaSupportPoint() const
- { assert(isInteriorNeumannCache_); return ipGlobalInside_; }
+ { assert(isInteriorBoundaryCache_); return ipGlobalInside_; }
//! returns the shape values at ip inside the element for interior boundaries
const std::vector<ShapeValue>& shapeValuesAtTpfaSupportPoint() const
- { assert(isInteriorNeumannCache_); return shapeValuesInside_; }
+ { assert(isInteriorBoundaryCache_); return shapeValuesInside_; }
private:
- bool isInteriorNeumannCache_{false};
+ bool isInteriorBoundaryCache_{false};
GlobalPosition ipGlobalInside_;
std::vector<ShapeValue> shapeValuesInside_;
};
--
GitLab