Skip to content

GitLab

Commit d98bfde1 authored by Timo Koch's avatar Timo Koch
Browse files

[pmflow] Separate velocity computation from velocity output class

parent a5fe8ab8
Expand all Hide whitespace changes
Inline Side-by-side
dumux/porousmediumflow/velocity.hh 0 → 100644
View file @ d98bfde1
This diff is collapsed.
dumux/porousmediumflow/velocityoutput.hh
View file @ d98bfde1
......@@ -25,6 +25,7 @@
#ifndef DUMUX_POROUSMEDIUMFLOW_VELOCITYOUTPUT_HH
#define DUMUX_POROUSMEDIUMFLOW_VELOCITYOUTPUT_HH
#include <memory>
#include <dune/common/float_cmp.hh>
#include <dune/geometry/referenceelements.hh>
......@@ -33,6 +34,7 @@
#include <dumux/io/velocityoutput.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/elementsolution.hh>
#include <dumux/porousmediumflow/velocity.hh>
namespace Dumux {
......@@ -67,6 +69,7 @@ class PorousMediumFlowVelocityOutput : public VelocityOutput<GridVariables>
using Problem = typename GridVolumeVariables::Problem;
using BoundaryTypes = typename Problem::Traits::BoundaryTypes;
using VelocityBackend = PorousMediumFlowVelocity<GridVariables, FluxVariables>;
public:
using VelocityVector = typename ParentType::VelocityVector;
......@@ -77,25 +80,12 @@ public:
* \param gridVariables The grid variables
*/
PorousMediumFlowVelocityOutput(const GridVariables& gridVariables)
: problem_(gridVariables.curGridVolVars().problem())
, fvGridGeometry_(gridVariables.fvGridGeometry())
, gridVariables_(gridVariables)
: gridVariables_(gridVariables) // TODO: can be removed after the deprecated calculateVelocity interface is removed
{
// check, if velocity output can be used (works only for cubes so far)
enableOutput_ = getParamFromGroup<bool>(problem_.paramGroup(), "Vtk.AddVelocity");
enableOutput_ = getParamFromGroup<bool>(gridVariables.curGridVolVars().problem().paramGroup(), "Vtk.AddVelocity");
if (enableOutput_)
{
// set the number of scvs the vertices are connected to
if (isBox && dim > 1)
{
// resize to the number of vertices of the grid
cellNum_.assign(fvGridGeometry_.gridView().size(dim), 0);
for (const auto& element : elements(fvGridGeometry_.gridView()))
for (unsigned int vIdx = 0; vIdx < element.subEntities(dim); ++vIdx)
++cellNum_[fvGridGeometry_.vertexMapper().subIndex(element, vIdx, dim)];
}
}
velocityBackend = std::make_unique<VelocityBackend>(gridVariables);
}
//! Returns whether or not velocity output is enabled.
......@@ -130,329 +120,14 @@ public:
const ElementFluxVarsCache& elemFluxVarsCache,
int phaseIdx) const override
{
using Velocity = typename VelocityVector::value_type;
if (!enableOutput_) return;
const auto geometry = element.geometry();
const Dune::GeometryType geomType = geometry.type();
// the upwind term to be used for the volume flux evaluation
auto upwindTerm = [phaseIdx](const auto& volVars) { return volVars.mobility(phaseIdx); };
if(isBox && dim == 1)
{
Velocity tmpVelocity(0.0);
tmpVelocity = (geometry.corner(1) - geometry.corner(0));
tmpVelocity /= tmpVelocity.two_norm();
for (auto&& scvf : scvfs(fvGeometry))
{
if (scvf.boundary())
continue;
// insantiate the flux variables
FluxVariables fluxVars;
fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
// get the volume flux divided by the area of the
// subcontrolvolume face in the reference element
Scalar localArea = scvfReferenceArea_(geomType, scvf.index());
Scalar flux = fluxVars.advectiveFlux(phaseIdx, upwindTerm) / localArea;
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
flux /= insideVolVars.extrusionFactor();
tmpVelocity *= flux;
const int eIdxGlobal = fvGridGeometry_.elementMapper().index(element);
velocity[eIdxGlobal] = tmpVelocity;
}
return;
}
// get the transposed Jacobian of the element mapping
const auto referenceElement = ReferenceElements::general(geomType);
const auto& localPos = referenceElement.position(0, 0);
const auto jacobianT2 = geometry.jacobianTransposed(localPos);
if(isBox)
{
using ScvVelocities = Dune::BlockVector<Velocity>;
ScvVelocities scvVelocities(fvGeometry.numScv());
scvVelocities = 0;
for (auto&& scvf : scvfs(fvGeometry))
{
if (scvf.boundary())
continue;
// local position of integration point
const auto localPosIP = geometry.local(scvf.ipGlobal());
// Transformation of the global normal vector to normal vector in the reference element
const auto jacobianT1 = geometry.jacobianTransposed(localPosIP);
const auto globalNormal = scvf.unitOuterNormal();
GlobalPosition localNormal(0);
jacobianT1.mv(globalNormal, localNormal);
localNormal /= localNormal.two_norm();
// instantiate the flux variables
FluxVariables fluxVars;
fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
// get the volume flux divided by the area of the
// subcontrolvolume face in the reference element
Scalar localArea = scvfReferenceArea_(geomType, scvf.index());
Scalar flux = fluxVars.advectiveFlux(phaseIdx, upwindTerm) / localArea;
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
flux /= insideVolVars.extrusionFactor();
// transform the volume flux into a velocity vector
Velocity tmpVelocity = localNormal;
tmpVelocity *= flux;
scvVelocities[scvf.insideScvIdx()] += tmpVelocity;
scvVelocities[scvf.outsideScvIdx()] += tmpVelocity;
}
// transform vertex velocities from local to global coordinates
for (auto&& scv : scvs(fvGeometry))
{
int vIdxGlobal = scv.dofIndex();
// calculate the subcontrolvolume velocity by the Piola transformation
Velocity scvVelocity(0);
jacobianT2.mtv(scvVelocities[scv.indexInElement()], scvVelocity);
scvVelocity /= geometry.integrationElement(localPos)*cellNum_[vIdxGlobal];
// add up the wetting phase subcontrolvolume velocities for each vertex
velocity[vIdxGlobal] += scvVelocity;
}
}
else
{
// For the number of scvfs per facet (mpfa) we simply obtain the number of
// corners of the first facet as prisms/pyramids are not supported here anyway
// -> for prisms/pyramids the number of scvfs would differ from facet to facet
static constexpr bool isMpfa = FVGridGeometry::discMethod == DiscretizationMethod::ccmpfa;
const int numScvfsPerFace = isMpfa ? element.template subEntity<1>(0).geometry().corners() : 1;
if (fvGeometry.numScvf() != element.subEntities(1)*numScvfsPerFace)
DUNE_THROW(Dune::NotImplemented, "Velocity output for non-conforming grids");
if (!geomType.isCube() && !geomType.isSimplex())
DUNE_THROW(Dune::NotImplemented, "Velocity output for other geometry types than cube and simplex");
// first we extract the corner indices for each scv for the CIV method
// for network grids there might be multiple intersection with the same geometryInInside
// we identify those by the indexInInside for now (assumes conforming grids at branching facets)
// here we keep track of them
std::vector<bool> handledScvf;
if (dim < dimWorld)
handledScvf.resize(element.subEntities(1), false);
// find the local face indices of the scvfs (for conforming meshes)
std::vector<unsigned int> scvfIndexInInside(fvGeometry.numScvf());
int localScvfIdx = 0;
for (const auto& intersection : intersections(fvGridGeometry_.gridView(), element))
{
if (dim < dimWorld)
if (handledScvf[intersection.indexInInside()])
continue;
if (intersection.neighbor() || intersection.boundary())
{
for (int i = 0; i < numScvfsPerFace; ++i)
scvfIndexInInside[localScvfIdx++] = intersection.indexInInside();
// for surface and network grids mark that we handled this face
if (dim < dimWorld)
handledScvf[intersection.indexInInside()] = true;
}
}
std::vector<Scalar> scvfFluxes(element.subEntities(1), 0.0);
localScvfIdx = 0;
for (auto&& scvf : scvfs(fvGeometry))
{
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
if (!scvf.boundary())
{
FluxVariables fluxVars;
fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
scvfFluxes[scvfIndexInInside[localScvfIdx]] += fluxVars.advectiveFlux(phaseIdx, upwindTerm)/insideVolVars.extrusionFactor();
}
else
{
auto bcTypes = problemBoundaryTypes_(element, scvf);
if (bcTypes.hasOnlyDirichlet())
{
FluxVariables fluxVars;
fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
scvfFluxes[scvfIndexInInside[localScvfIdx]] += fluxVars.advectiveFlux(phaseIdx, upwindTerm)/insideVolVars.extrusionFactor();
}
}
// increment scvf counter
localScvfIdx++;
}
// Correct boundary fluxes in case of Neumann conditions.
// In this general setting, it would be very difficult to
// calculate correct phase, i.e., volume, fluxes from arbitrary
// Neumann conditions. We approximate the Neumann flux by the
// flux on the opposite face. For extremely distorted grids this can
// lead to unexpected results (but then TPFA also leads to unexpected results).
localScvfIdx = 0;
for (auto&& scvf : scvfs(fvGeometry))
{
if (scvf.boundary())
{
auto bcTypes = problemBoundaryTypes_(element, scvf);
if (bcTypes.hasNeumann())
{
// check if we have Neumann no flow, we can just use 0
const auto neumannFlux = Deprecated::neumann(problem_, element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
using NumEqVector = std::decay_t<decltype(neumannFlux)>;
if (Dune::FloatCmp::eq<NumEqVector, Dune::FloatCmp::CmpStyle::absolute>(neumannFlux, NumEqVector(0.0), 1e-30))
scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0;
// cubes
else if (dim == 1 || geomType.isCube())
{
const auto fIdx = scvfIndexInInside[localScvfIdx];
const auto fIdxOpposite = fIdx%2 ? fIdx-1 : fIdx+1;
scvfFluxes[fIdx] = -scvfFluxes[fIdxOpposite];
}
// simplices
else if (geomType.isSimplex())
scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0;
}
}
// increment scvf counter
localScvfIdx++;
}
Velocity refVelocity;
// cubes: On the reference element simply average over opposite fluxes
// note that this is equal to a corner velocity interpolation method
if (dim == 1 || geomType.isCube())
{
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (scvfFluxes[2*i + 1] - scvfFluxes[2*i]);
}
// simplices: Raviart-Thomas-0 interpolation evaluated at the cell center
else if (geomType.isSimplex())
{
for (int dimIdx = 0; dimIdx < dim; dimIdx++)
{
refVelocity[dimIdx] = -scvfFluxes[dim - 1 - dimIdx];
for (int fIdx = 0; fIdx < dim + 1; fIdx++)
refVelocity[dimIdx] += scvfFluxes[fIdx]/(dim + 1);
}
}
// 3D prism and pyramids
else
DUNE_THROW(Dune::NotImplemented, "velocity output for cell-centered and prism/pyramid");
Velocity scvVelocity(0);
jacobianT2.mtv(refVelocity, scvVelocity);
scvVelocity /= geometry.integrationElement(localPos);
int eIdxGlobal = fvGridGeometry_.elementMapper().index(element);
velocity[eIdxGlobal] = scvVelocity;
} // cell-centered
}
private:
// The area of a subcontrolvolume face in a reference element.
// The 3d non-cube values have been calculated with quadrilateralArea3D
// of boxfvelementgeometry.hh.
static Scalar scvfReferenceArea_(Dune::GeometryType geomType, int fIdx)
{
if (dim == 1 || geomType.isCube())
{
return 1.0/(1 << (dim-1));
}
else if (geomType.isTriangle())
{
static const Scalar faceToArea[] = {0.372677996249965,
0.372677996249965,
0.235702260395516};
return faceToArea[fIdx];
}
else if (geomType.isTetrahedron())
{
static const Scalar faceToArea[] = {0.102062072615966,
0.102062072615966,
0.058925565098879,
0.102062072615966,
0.058925565098879,
0.058925565098879};
return faceToArea[fIdx];
}
else if (geomType.isPyramid())
{
static const Scalar faceToArea[] = {0.130437298687488,
0.130437298687488,
0.130437298687488,
0.130437298687488,
0.150923085635624,
0.1092906420717,
0.1092906420717,
0.0781735959970571};
return faceToArea[fIdx];
}
else if (geomType.isPrism())
{
static const Scalar faceToArea[] = {0.166666666666667,
0.166666666666667,
0.166666666666667,
0.186338998124982,
0.186338998124982,
0.117851130197758,
0.186338998124982,
0.186338998124982,
0.117851130197758};
return faceToArea[fIdx];
}
else {
DUNE_THROW(Dune::NotImplemented, "scvf area for unknown GeometryType");
}
if (enableOutput_)
velocityBackend->calculateVelocity(velocity, element, fvGeometry, elemVolVars, elemFluxVarsCache, phaseIdx);
}
private:
// The following SFINAE enable_if usage allows compilation, even if only a
//
// boundaryTypes(const Element&, const scv&)
//
// is provided in the problem file. In that case, the compiler cannot detect
// (without additional measures like "using...") the signature
//
// boundaryTypes(const Element&, const scvf&)
//
// in the problem base class. Therefore, calls to this method trigger a
// compiler error. However, that call is needed for calculating velocities
// if the cell-centered discretization is used. By proceeding as in the
// following lines, that call will only be compiled if cell-centered
// actually is used.
template <bool enable = isBox, typename std::enable_if_t<!enable, int> = 0>
BoundaryTypes problemBoundaryTypes_(const Element& element, const SubControlVolumeFace& scvf) const
{ return problem_.boundaryTypes(element, scvf); }
//! we should never call this method for box models
template <bool enable = isBox, typename std::enable_if_t<enable, int> = 0>
BoundaryTypes problemBoundaryTypes_(const Element& element, const SubControlVolumeFace& scvf) const
{ return BoundaryTypes(); }
const Problem& problem_;
const FVGridGeometry& fvGridGeometry_;
const GridVariables& gridVariables_;
const GridVariables& gridVariables_; // TODO: can be removed after the deprecated calculateVelocity interface is removed
bool enableOutput_;
std::vector<int> cellNum_;
std::unique_ptr<VelocityBackend> velocityBackend;
};
} // end namespace Dumux
......
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