From aca0fe325c3006ee09c32354270f378c8b1543d3 Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser" <dennis.glaeser@iws.uni-stuttgart.de>
Date: Wed, 30 Sep 2020 14:11:33 +0200
Subject: [PATCH] [facet][tpfa] support gravity on surface grids
---
.../facet/cellcentered/tpfa/darcyslaw.hh | 72 +++++++++++++++----
1 file changed, 60 insertions(+), 12 deletions(-)
diff --git a/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh b/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh
index a1a032e730..937c76b9a3 100644
--- a/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh
+++ b/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh
@@ -367,6 +367,9 @@ public:
//! Export transmissibility storage type
using AdvectionTransmissibilityContainer = std::array<Scalar, 2>;
+ //! Export the type used for the gravity coefficients
+ using GravityCoefficients = std::array<Scalar, 1>;
+
//! update subject to a given problem
template< class Problem, class ElementVolumeVariables >
void updateAdvection(const Problem& problem,
@@ -375,7 +378,7 @@ public:
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace &scvf)
{
- tij_ = AdvectionType::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf);
+ tij_ = AdvectionType::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, g_);
}
//! We use the same name as in the TpfaDarcysLawCache so
@@ -393,8 +396,13 @@ public:
Scalar advectionTijFacet() const
{ return tij_[facetTijIdx]; }
+ //! return the coefficients for the computation of gravity at the scvf
+ const GravityCoefficients& gravityCoefficients() const
+ { return g_; }
+
private:
AdvectionTransmissibilityContainer tij_;
+ GravityCoefficients g_;
};
/*!
@@ -436,24 +444,49 @@ class CCTpfaFacetCouplingDarcysLawImpl<ScalarType, GridGeometry, /*isNetwork*/tr
int phaseIdx,
const ElementFluxVarsCache& elemFluxVarsCache)
{
- static const Scalar gravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
- if (gravity)
- DUNE_THROW(Dune::NotImplemented, "Gravity for darcys law with facet coupling on surface grids");
-
if (!problem.couplingManager().isOnInteriorBoundary(element, scvf))
return TpfaDarcysLaw::flux(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarsCache);
+ // On surface grids only xi = 1.0 can be used, as the coupling condition
+ // for xi != 1.0 does not generalize for surface grids where there can be
+ // seveal neighbor meeting at a branching point.
+ static const Scalar xi = getParamFromGroup<Scalar>(problem.paramGroup(), "FacetCoupling.Xi", 1.0);
+ if (Dune::FloatCmp::ne(xi, 1.0, 1e-6))
+ DUNE_THROW(Dune::InvalidStateException, "Xi != 1.0 cannot be used on surface grids");
+
// Obtain inside and fracture pressures
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+ const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf);
const auto pInside = insideVolVars.pressure(phaseIdx);
- const auto pFacet = problem.couplingManager().getLowDimVolVars(element, scvf).pressure(phaseIdx);
+ const auto pFacet = facetVolVars.pressure(phaseIdx);
- // return flux
+ // compute and return flux
const auto& fluxVarsCache = elemFluxVarsCache[scvf];
- if (scvf.boundary())
- return fluxVarsCache.advectionTijInside()*pInside + fluxVarsCache.advectionTijFacet()*pFacet;
- else
- return fluxVarsCache.advectionTijInside()*pInside + fluxVarsCache.advectionTijFacet()*pFacet;
+ Scalar flux = fluxVarsCache.advectionTijInside()*pInside + fluxVarsCache.advectionTijFacet()*pFacet;
+
+ static const Scalar gravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
+ if (gravity)
+ {
+ // compute alpha := n^T*K*g and add to flux (use arithmetic mean for density)
+ const auto& g = problem.spatialParams().gravity(scvf.ipGlobal());
+ const auto rho = 0.5*(insideVolVars.density(phaseIdx) + facetVolVars.density(phaseIdx));
+ const auto rhoTimesArea = rho*Extrusion::area(scvf);
+ const auto alpha_inside = rhoTimesArea*insideVolVars.extrusionFactor()
+ *vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g);
+
+ flux += alpha_inside;
+
+ // maybe add further gravitational contributions
+ if ( !scvf.boundary() && problem.interiorBoundaryTypes(element, scvf).hasOnlyNeumann() )
+ {
+ const auto alpha_facet = rhoTimesArea*insideVolVars.extrusionFactor()
+ *vtmv(scvf.unitOuterNormal(), facetVolVars.permeability(), g);
+
+ flux -= fluxVarsCache.gravityCoefficients()[0]*(alpha_inside - alpha_facet);
+ }
+ }
+
+ return flux;
}
// The flux variables cache has to be bound to an element prior to flux calculations
@@ -464,6 +497,20 @@ class CCTpfaFacetCouplingDarcysLawImpl<ScalarType, GridGeometry, /*isNetwork*/tr
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf)
+ {
+ typename Cache::GravityCoefficients g;
+ return calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, g);
+ }
+
+ // This overload additionally receives a container in which the coefficients required
+ // for the computation of the gravitational acceleration ar the scvf are stored
+ template< class Problem, class ElementVolumeVariables >
+ static TijContainer calculateTransmissibility(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ typename Cache::GravityCoefficients& g)
{
TijContainer tij;
if (!problem.couplingManager().isCoupled(element, scvf))
@@ -504,7 +551,8 @@ class CCTpfaFacetCouplingDarcysLawImpl<ScalarType, GridGeometry, /*isNetwork*/tr
*vtmv(scvf.unitOuterNormal(), facetVolVars.permeability(), scvf.unitOuterNormal());
// TODO: check for division by zero??
- tij[Cache::insideTijIdx] = wFacet*wIn/(wIn+wFacet);
+ g[0] = wIn/(wIn+wFacet);
+ tij[Cache::insideTijIdx] = wFacet*g[0];
tij[Cache::facetTijIdx] = -tij[Cache::insideTijIdx];
}
else if (iBcTypes.hasOnlyDirichlet())
--
GitLab