From 3133e81eda938ccfc4b52cfd4c2f23ff8a02a9f1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= <dennis.glaeser@iws.uni-stuttgart.de>
Date: Mon, 12 Nov 2018 19:11:13 +0100
Subject: [PATCH] [facet][tpfa][darcyslaw] implement gravity handling for xi !=
1.0
---
.../facet/cellcentered/tpfa/darcyslaw.hh | 106 +++++++++++++-----
test/multidomain/facet/1p_1p/CMakeLists.txt | 1 +
2 files changed, 81 insertions(+), 26 deletions(-)
diff --git a/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh b/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh
index 0875ddf35f..921630cce9 100644
--- a/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh
+++ b/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh
@@ -90,6 +90,9 @@ public:
//! Export transmissibility storage type
using AdvectionTransmissibilityContainer = std::array<Scalar, 3>;
+ //! Export the type used for the gravity coefficients
+ using GravityCoefficients = std::array<Scalar, 2>;
+
//! update subject to a given problem
template< class Problem, class ElementVolumeVariables >
void updateAdvection(const Problem& problem,
@@ -98,26 +101,35 @@ 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
//! that this cache and the law implementation for non-coupled
//! models can be reused here on facets that do not lie on an
//! interior boundary, i.e. do not coincide with a facet element
- Scalar advectionTij() const { return tij_[insideTijIdx]; }
+ Scalar advectionTij() const
+ { return tij_[insideTijIdx]; }
//! returns the transmissibility associated with the inside cell
- Scalar advectionTijInside() const { return tij_[insideTijIdx]; }
+ Scalar advectionTijInside() const
+ { return tij_[insideTijIdx]; }
//! returns the transmissibility associated with the outside cell
- Scalar advectionTijOutside() const {return tij_[outsideTijIdx]; }
+ Scalar advectionTijOutside() const
+ { return tij_[outsideTijIdx]; }
//! returns the transmissibility associated with the outside cell
- Scalar advectionTijFacet() const {return tij_[facetTijIdx]; }
+ Scalar advectionTijFacet() const
+ { return tij_[facetTijIdx]; }
+
+ //! return the coefficients for the computation of gravity at the scvf
+ const GravityCoefficients& gravityCoefficients() const
+ { return g_; }
private:
std::array<Scalar, 3> tij_;
+ GravityCoefficients g_;
};
/*!
@@ -188,27 +200,33 @@ class CCTpfaFacetCouplingDarcysLawImpl<ScalarType, FVGridGeometry, /*isNetwork*/
static const Scalar gravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
if (gravity)
{
- // this is inconsistent for xi != 1
- static const Scalar xi = getParamFromGroup<Scalar>(problem.paramGroup(), "FacetCoupling.Xi", 1.0);
- if (xi != 1.0)
- DUNE_THROW(Dune::NotImplemented, "Gravitational acceleration for facet coupling and xi != 1.0");
-
- // compute alpha := n^T*K*g
+ // compute alpha := n^T*K*g and add to flux (use arithmetic mean for density)
const auto& g = problem.gravityAtPos(scvf.ipGlobal());
- const auto alpha_inside = vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g)*insideVolVars.extrusionFactor();
-
- // for the density, use arithmetic average
const auto rho = 0.5*(insideVolVars.density(phaseIdx) + facetVolVars.density(phaseIdx));
- flux += rho*scvf.area()*alpha_inside;
+ const auto rhoTimesArea = rho*scvf.area();
+ const auto alpha_inside = rhoTimesArea*insideVolVars.extrusionFactor()
+ *vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g);
- // maybe add K-weighted gravitational contribution
+ flux += alpha_inside;
if (!scvf.boundary())
{
const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
- const auto wFacet = computeFacetTransmissibility_(insideVolVars, facetVolVars, scvf);
- const auto alpha_outside = vtmv(scvf.unitOuterNormal(), outsideVolVars.permeability(), g)*outsideVolVars.extrusionFactor();
- flux += rho*fluxVarsCache.advectionTijInside()/wFacet*(alpha_inside - alpha_outside);
+ // add further gravitational contributions
+ if ( problem.interiorBoundaryTypes(element, scvf).hasOnlyNeumann() )
+ {
+ static const Scalar xi = getParamFromGroup<Scalar>(problem.paramGroup(), "FacetCoupling.Xi", 1.0);
+ const auto alpha_facet = rhoTimesArea*insideVolVars.extrusionFactor()
+ *vtmv(scvf.unitOuterNormal(), facetVolVars.permeability(), g);
+ const auto alpha_outside = rhoTimesArea*outsideVolVars.extrusionFactor()
+ *vtmv(scvf.unitOuterNormal(), outsideVolVars.permeability(), g);
+
+ flux -= fluxVarsCache.gravityCoefficients()[0]*(xi*alpha_inside - alpha_facet + (1.0 - xi)*alpha_outside);
+ flux += fluxVarsCache.gravityCoefficients()[1]*(xi*alpha_outside - alpha_facet + (1.0 - xi)*alpha_inside);
+ }
+
+ // add outside contribution
+ flux += fluxVarsCache.advectionTijOutside()*outsideVolVars.pressure(phaseIdx);
}
return flux;
@@ -226,6 +244,20 @@ class CCTpfaFacetCouplingDarcysLawImpl<ScalarType, FVGridGeometry, /*isNetwork*/
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))
@@ -260,8 +292,9 @@ class CCTpfaFacetCouplingDarcysLawImpl<ScalarType, FVGridGeometry, /*isNetwork*/
// where \f$\gamma$\f denotes the domain living on the facets and \f$\bar{\mathbf{u}}$\f are
// intermediate face unknowns in the matrix domain. Equivalently, flux continuity reads:
// \f$\mathbf{A} \bar{\mathbf{u}} = \mathbf{B} \mathbf{u} + \mathbf{M} \mathbf{u}_\gamma\f$.
- // Combining the two, we can eliminate the intermediate unknowns and compute the transmissibilities.
- if (!scvf.boundary() && xi != 1.0)
+ // Combining the two, we can eliminate the intermediate unknowns and compute the transmissibilities
+ // that allow the description of the fluxes as functions of the cell and Dirichlet pressures only.
+ if (!scvf.boundary())
{
const auto outsideScvIdx = scvf.outsideScvIdx();
const auto& outsideVolVars = elemVolVars[outsideScvIdx];
@@ -270,11 +303,32 @@ class CCTpfaFacetCouplingDarcysLawImpl<ScalarType, FVGridGeometry, /*isNetwork*/
outsideVolVars.permeability(),
outsideVolVars.extrusionFactor());
- const Scalar factor = wIn * wFacet / ( wIn * wOut * ( 2.0 * xi - 1.0 ) + wFacet * ( xi * ( wIn + wOut ) + wFacet ) );
-
- tij[Cache::insideTijIdx] = factor * ( wOut * xi + wFacet );
- tij[Cache::outsideTijIdx] = factor * ( wOut * ( 1.0 - xi ) );
- tij[Cache::facetTijIdx] = factor * ( - wOut - wFacet );
+ if (xi != 1.0)
+ {
+ // The gravity coefficients are the first row of the inverse of the A matrix in the local eq system
+ // multiplied with wIn. Note that we never compute the inverse but use an optimized implementation below.
+ // The A matrix has the following coefficients:
+ // A = | xi*wIn + wFacet, (xi - 1.0)*wOut | -> AInv = 1/detA | xi*wOut + wFacet, -(xi - 1.0)*wOut |
+ // | wIn*(xi - 1.0) , xi*wOut + wFacet | | -wIn*(xi - 1.0) , xi*wIn + wFacet |
+ const Scalar xiMinusOne = (xi - 1.0);
+ const Scalar a01 = xiMinusOne*wOut;
+ const Scalar a11 = xi*wOut + wFacet;
+ const Scalar detA = (xi*wIn + wFacet)*a11 - xiMinusOne*wIn*a01;
+ g[0] = wIn*a11/detA; g[1] = -wIn*a01/detA;
+
+ // optimized implementation: factorization obtained using sympy
+ const Scalar factor = wIn * wFacet / ( wIn * wOut * ( 2.0 * xi - 1.0 ) + wFacet * ( xi * ( wIn + wOut ) + wFacet ) );
+ tij[Cache::insideTijIdx] = factor * ( wOut * xi + wFacet );
+ tij[Cache::outsideTijIdx] = factor * ( wOut * ( 1.0 - xi ) );
+ tij[Cache::facetTijIdx] = factor * ( - wOut - wFacet );
+ }
+ else
+ {
+ g[0] = wIn/(wIn+wFacet); g[1] = 0.0;
+ tij[Cache::insideTijIdx] = wFacet*g[0];
+ tij[Cache::facetTijIdx] = -tij[Cache::insideTijIdx];
+ tij[Cache::outsideTijIdx] = 0.0;
+ }
}
else
{
diff --git a/test/multidomain/facet/1p_1p/CMakeLists.txt b/test/multidomain/facet/1p_1p/CMakeLists.txt
index f15f388ce1..fc3443e51e 100644
--- a/test/multidomain/facet/1p_1p/CMakeLists.txt
+++ b/test/multidomain/facet/1p_1p/CMakeLists.txt
@@ -1,3 +1,4 @@
add_subdirectory(analytical)
+add_subdirectory(gravity)
add_subdirectory(linearprofile)
add_subdirectory(threedomain)
--
GitLab