Commit 3133e81e authored by Dennis Gläser's avatar Dennis Gläser Committed by Timo Koch
Browse files

[facet][tpfa][darcyslaw] implement gravity handling for xi != 1.0

parent 529f0571
......@@ -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
{
......
add_subdirectory(analytical)
add_subdirectory(gravity)
add_subdirectory(linearprofile)
add_subdirectory(threedomain)
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