From 4c43ba827aa29c00d1f789ad677469ef77194470 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Thu, 12 Jul 2018 19:05:33 +0200
Subject: [PATCH] [periodic] Implement periodic boundaries for cctpfa
---
.../cellcentered/tpfa/darcyslaw.hh | 12 ++++---
.../cellcentered/tpfa/fvelementgeometry.hh | 34 ++++++-------------
.../cellcentered/tpfa/fvgridgeometry.hh | 16 ++++++---
3 files changed, 30 insertions(+), 32 deletions(-)
diff --git a/dumux/discretization/cellcentered/tpfa/darcyslaw.hh b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
index c80cc8c692..109ac32bc9 100644
--- a/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
+++ b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
@@ -191,10 +191,12 @@ class CCTpfaDarcysLaw<ScalarType, FVGridGeometry, /*isNetwork*/ false>
{
const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
const auto outsideK = outsideVolVars.permeability();
- const auto outsideTi = computeTpfaTransmissibility(scvf, outsideScv, outsideK, outsideVolVars.extrusionFactor());
+ const auto outsideTi = fvGeometry.fvGridGeometry().isPeriodic()
+ ? computeTpfaTransmissibility(fvGeometry.flipScvf(scvf.index()), outsideScv, outsideK, outsideVolVars.extrusionFactor())
+ : -1.0*computeTpfaTransmissibility(scvf, outsideScv, outsideK, outsideVolVars.extrusionFactor());
const auto alpha_outside = vtmv(scvf.unitOuterNormal(), outsideK, g)*outsideVolVars.extrusionFactor();
- flux += rho*tij/outsideTi*(alpha_inside - alpha_outside);
+ flux -= rho*tij/outsideTi*(alpha_inside - alpha_outside);
}
return flux;
@@ -241,9 +243,9 @@ class CCTpfaDarcysLaw<ScalarType, FVGridGeometry, /*isNetwork*/ false>
// refers to the scv of our element, so we use the scv method
const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
const auto& outsideVolVars = elemVolVars[outsideScvIdx];
- const Scalar tj = -1.0*computeTpfaTransmissibility(scvf, outsideScv,
- getPermeability_(problem, outsideVolVars, scvf.ipGlobal()),
- outsideVolVars.extrusionFactor());
+ const Scalar tj = fvGeometry.fvGridGeometry().isPeriodic()
+ ? computeTpfaTransmissibility(fvGeometry.flipScvf(scvf.index()), outsideScv, getPermeability_(problem, outsideVolVars, scvf.ipGlobal()), outsideVolVars.extrusionFactor())
+ : -1.0*computeTpfaTransmissibility(scvf, outsideScv, getPermeability_(problem, outsideVolVars, scvf.ipGlobal()), outsideVolVars.extrusionFactor());
// harmonic mean (check for division by zero!)
// TODO: This could lead to problems!? Is there a better way to do this?
diff --git a/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
index fb1cf2a7fd..c1e4b4388b 100644
--- a/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
+++ b/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
@@ -281,8 +281,10 @@ public:
std::vector<IndexType> handledNeighbors;
handledNeighbors.reserve(element.subEntities(1));
+
for (const auto& intersection : intersections(fvGridGeometry().gridView(), element))
{
+ // for inner intersections and periodic (according to grid interface) intersections make neighbor geometry
if (intersection.neighbor())
{
const auto outside = intersection.outside();
@@ -297,7 +299,8 @@ public:
}
}
- if (dim < dimWorld)
+ // build flip index set for network, surface, and periodic grids
+ if (dim < dimWorld || fvGridGeometry().isPeriodic())
{
flippedScvfIndices_.resize(scvfs_.size());
for (unsigned int localScvfIdx = 0; localScvfIdx < scvfs_.size(); ++localScvfIdx)
@@ -332,21 +335,6 @@ public:
}
}
}
-
- //! TODO Check if user added additional DOF dependencies, i.e. the residual of DOF globalI depends
- //! on additional DOFs not included in the discretization schemes' occupation pattern
- // const auto globalI = fvGridGeometry().elementMapper().index(element);
- // const auto& additionalDofDependencies = problem.getAdditionalDofDependencies(globalI);
- // if (!additionalDofDependencies.empty())
- // {
- // neighborScvs_.reserve(neighborScvs_.size() + additionalDofDependencies.size());
- // neighborScvIndices_.reserve(neighborScvIndices_.size() + additionalDofDependencies.size());
- // for (auto globalJ : additionalDofDependencies)
- // {
- // neighborScvs_.emplace_back(fvGridGeometry().element(globalJ).geometry(), globalJ);
- // neighborScvIndices_.emplace_back(globalJ);
- // }
- // }
}
//! Binding of an element preparing the geometries only inside the element
@@ -418,7 +406,6 @@ private:
continue;
const auto& scvfNeighborVolVarIndices = neighborVolVarIndices[scvfCounter];
-
if (intersection.neighbor() || intersection.boundary())
{
ScvfGridIndexStorage scvIndices;
@@ -429,7 +416,7 @@ private:
intersection.geometry(),
scvFaceIndices[scvfCounter],
scvIndices,
- intersection.boundary());
+ intersection.boundary() && !intersection.neighbor());
scvfIndices_.emplace_back(scvFaceIndices[scvfCounter]);
scvfCounter++;
@@ -466,14 +453,14 @@ private:
if (handledScvf[intersection.indexInInside()])
continue;
+ // this catches inner and periodic scvfs
const auto& scvfNeighborVolVarIndices = neighborVolVarIndices[scvfCounter];
-
- if (intersection.neighbor())
+ if (scvfNeighborVolVarIndices[0] < fvGridGeometry().gridView().size(0))
{
// only create subcontrol faces where the outside element is the bound element
if (dim == dimWorld)
{
- if (intersection.outside() == *elementPtr_)
+ if (scvfNeighborVolVarIndices[0] == fvGridGeometry().elementMapper().index(*elementPtr_))
{
ScvfGridIndexStorage scvIndices({eIdx, scvfNeighborVolVarIndices[0]});
neighborScvfs_.emplace_back(intersection,
@@ -509,7 +496,6 @@ private:
break;
}
}
-
}
// for surface and network grids mark that we handled this face
@@ -517,9 +503,11 @@ private:
handledScvf[intersection.indexInInside()] = true;
scvfCounter++;
}
+
+ // only increase counter for boundary intersections
+ // (exclude periodic boundaries according to dune grid interface, they have been handled in neighbor==true)
else if (intersection.boundary())
{
- // for surface and network grids mark that we handled this face
if (dim < dimWorld)
handledScvf[intersection.indexInInside()] = true;
scvfCounter++;
diff --git a/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh
index 71ea360571..a7d381f2c7 100644
--- a/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh
+++ b/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh
@@ -208,9 +208,13 @@ public:
for (const auto& intersection : intersections(this->gridView(), element))
{
- // inner sub control volume faces
+ // inner sub control volume faces (includes periodic boundaries)
if (intersection.neighbor())
{
+ // update the grid geometry if we have periodic boundaries
+ if (intersection.boundary())
+ this->setPeriodic();
+
if (dim == dimWorld)
{
const auto nIdx = this->elementMapper().index(intersection.outside());
@@ -257,8 +261,8 @@ public:
scvfIndicesOfScv_[eIdx] = scvfsIndexSet;
}
- // Make the flip index set for network and surface grids
- if (dim < dimWorld)
+ // Make the flip index set for network, surface, and periodic grids
+ if (dim < dimWorld || this->isPeriodic())
{
flipScvfIndices_.resize(scvfs_.size());
for (auto&& scvf : scvfs_)
@@ -467,9 +471,13 @@ public:
for (const auto& intersection : intersections(this->gridView(), element))
{
- // inner sub control volume faces
+ // inner sub control volume faces (includes periodic boundaries)
if (intersection.neighbor())
{
+ // update the grid geometry if we have periodic boundaries
+ if (intersection.boundary())
+ this->setPeriodic();
+
if (dim == dimWorld)
{
scvfsIndexSet.push_back(numScvf_++);
--
GitLab