Commit a43a81c4 authored by Dennis Gläser's avatar Dennis Gläser
Browse files

Merge branch 'feature/pass-fluxcache-to-neumann' into 'master'

Feature/pass fluxcache to neumann

Closes #660

See merge request !1626
parents 86f9a7ab a386f36b
......@@ -28,6 +28,7 @@
#include <dune/geometry/type.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/deprecated.hh>
#include <dumux/discretization/method.hh>
namespace Dumux {
......@@ -163,10 +164,14 @@ public:
// get the fvGeometry and elementVolVars needed for the bc and source interfaces
auto fvGeometry = localView(*fvGridGeometry_);
fvGeometry.bindElement(element);
fvGeometry.bind(element);
auto elemVolVars = localView(gridVariables_->curGridVolVars());
elemVolVars.bindElement(element, fvGeometry, sol);
elemVolVars.bind(element, fvGeometry, sol);
// elemFluxVarsCache for neumann interface
auto elemFluxVarsCache = localView(gridVariables_->gridFluxVarsCache());
elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
//! Check if we have to refine around a source term
if (refineAtSource_)
......@@ -208,7 +213,7 @@ public:
// we are on a pure Neumann boundary
else if(refineAtFluxBC_)
{
const auto fluxes = problem_->neumann(element, fvGeometry, elemVolVars, scvf);
const auto fluxes = Deprecated::neumann(*problem_, element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
if (fluxes.infinity_norm() > eps_)
{
indicatorVector_[eIdx] = true;
......@@ -242,7 +247,7 @@ public:
//! check if scvf is on Neumann boundary
if (scvf.boundary() && bcTypes[scvf.insideScvIdx()].hasNeumann())
{
const auto fluxes = problem_->neumann(element, fvGeometry, elemVolVars, scvf);
const auto fluxes = Deprecated::neumann(*problem_, element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
if (fluxes.infinity_norm() > eps_)
{
indicatorVector_[eIdx] = true;
......
......@@ -28,6 +28,7 @@
#include <dune/geometry/type.hh>
#include <dune/istl/matrix.hh>
#include <dumux/common/deprecated.hh>
#include <dumux/common/properties.hh>
#include <dumux/assembly/fvlocalresidual.hh>
......@@ -111,7 +112,7 @@ public:
// are enforced strongly by replacing the residual entry afterwards.
if (bcTypes.hasNeumann())
{
auto neumannFluxes = problem.neumann(element, fvGeometry, elemVolVars, scvf);
auto neumannFluxes = Deprecated::neumann(problem, element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
// multiply neumann fluxes with the area and the extrusion factor
neumannFluxes *= scvf.area()*elemVolVars[scv].extrusionFactor();
......
......@@ -25,6 +25,7 @@
#ifndef DUMUX_CC_LOCAL_RESIDUAL_HH
#define DUMUX_CC_LOCAL_RESIDUAL_HH
#include <dumux/common/deprecated.hh>
#include <dumux/common/reservedblockvector.hh>
#include <dumux/common/properties.hh>
#include <dumux/assembly/fvlocalresidual.hh>
......@@ -96,7 +97,7 @@ public:
// Neumann and Robin ("solution dependent Neumann") boundary conditions
else if (bcTypes.hasNeumann() && !bcTypes.hasDirichlet())
{
auto neumannFluxes = problem.neumann(element, fvGeometry, elemVolVars, scvf);
auto neumannFluxes = Deprecated::neumann(problem, element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
// multiply neumann fluxes with the area and the extrusion factor
const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
......
......@@ -76,6 +76,49 @@ auto effectiveThermalConductivity(const VV& volVars,
{
return ETC::effectiveThermalConductivity(volVars);
}
// support old interface of the neumann() function on problems
template<class E, class FVEG, class EVV, class EFVC>
class HasNewNeumannIF
{
using SCVF = typename FVEG::SubControlVolumeFace;
public:
template<class P>
auto operator()(P&& p) -> decltype(p.neumann(std::declval<const E&>(),
std::declval<const FVEG&>(),
std::declval<const EVV&>(),
std::declval<const EFVC&>(),
std::declval<const SCVF&>()))
{}
};
template<class P, class E, class FVEG, class EVV, class EFVC,
typename std::enable_if_t<!decltype(isValid(HasNewNeumannIF<E, FVEG, EVV, EFVC>()).template check<P>())::value, int> = 0>
auto DUNE_DEPRECATED_MSG("Use new neumann() interface (see common/fvproblem.hh) that additionally receives the element flux variables cache in your problem!. Will be removed after 3.1 release")
neumann(const P& problem,
const E& element,
const FVEG& fvGeometry,
const EVV& elemVolVars,
const EFVC& elemFluxVarsCache,
const typename FVEG::SubControlVolumeFace& scvf)
{
return problem.neumann(element, fvGeometry, elemVolVars, scvf);
}
template<class P, class E, class FVEG, class EVV, class EFVC,
typename std::enable_if_t<decltype(isValid(HasNewNeumannIF<E, FVEG, EVV, EFVC>()).template check<P>())::value, int> = 0>
auto neumann(const P& problem,
const E& element,
const FVEG& fvGeometry,
const EVV& elemVolVars,
const EFVC& elemFluxVarsCache,
const typename FVEG::SubControlVolumeFace& scvf)
{
return problem.neumann(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
}
////////////////////////////////////////////////////////
////////////////////////////////////////////////////////
......
......@@ -57,17 +57,17 @@ class FVProblem
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
enum {
dim = GridView::dimension
};
enum { dim = GridView::dimension };
using PointSource = GetPropType<TypeTag, Properties::PointSource>;
using PointSourceHelper = GetPropType<TypeTag, Properties::PointSourceHelper>;
using PointSourceMap = std::map<std::pair<std::size_t, std::size_t>,
std::vector<PointSource> >;
using PointSourceMap = std::map< std::pair<std::size_t, std::size_t>,
std::vector<PointSource> >;
using GridVolumeVariables = GetPropType<TypeTag, Properties::GridVolumeVariables>;
using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
......@@ -88,7 +88,6 @@ public:
using PrimaryVariables = FVProblem::PrimaryVariables;
using NumEqVector = FVProblem::NumEqVector;
using BoundaryTypes = FVProblem::BoundaryTypes;
};
/*!
......@@ -275,14 +274,16 @@ public:
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param elemVolVars All volume variables for the element
* \param elemFluxVarsCache Flux variables caches for all faces in stencil
* \param scvf The sub control volume face
*
* Negative values mean influx.
* E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
* E.g. for the mass balance that would be the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
*/
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
// forward it to the interface with only the global position
......
......@@ -27,6 +27,7 @@
#include <dune/geometry/type.hh>
#include <dumux/common/deprecated.hh>
#include <dumux/common/properties.hh>
#include <dumux/assembly/fvlocalresidual.hh>
......@@ -109,7 +110,7 @@ public:
// Neumann and Robin ("solution dependent Neumann") boundary conditions
if (bcTypes.hasNeumann() && !bcTypes.hasDirichlet())
{
auto neumannFluxes = problem.neumann(element, fvGeometry, elemVolVars, scvf);
auto neumannFluxes = Deprecated::neumann(problem, element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
// multiply neumann fluxes with the area and the extrusion factor
neumannFluxes *= scvf.area()*elemVolVars[scv].extrusionFactor();
......
......@@ -29,6 +29,7 @@
#include <dune/geometry/referenceelements.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/deprecated.hh>
#include <dumux/io/velocityoutput.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/elementsolution.hh>
......@@ -310,7 +311,7 @@ public:
if (bcTypes.hasNeumann())
{
// check if we have Neumann no flow, we can just use 0
const auto neumannFlux = problem_.neumann(element, fvGeometry, elemVolVars, scvf);
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;
......
......@@ -115,7 +115,10 @@ class DamBreakProblem : public ShallowWaterProblem<TypeTag>
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using NeumannFluxes = GetPropType<TypeTag, Properties::NumEqVector>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GridView = GetPropType<TypeTag, Properties::GridView>;
......@@ -220,19 +223,20 @@ public:
* \param element
* \param fvGeometry
* \param elemVolVars
* \param elemFluxVarsCache
* \param scvf
*/
NeumannFluxes neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
NeumannFluxes values(0.0);
//we need the Riemann invariants to compute the values depending of the boundary type
//since we use a weak imposition we do not have a dirichlet value. We impose fluxes
//based on q,h, etc. computed with the Riemann invariants
// we need the Riemann invariants to compute the values depending of the boundary type
// since we use a weak imposition we do not have a dirichlet value. We impose fluxes
// based on q,h, etc. computed with the Riemann invariants
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
const auto& insideVolVars = elemVolVars[insideScv];
const auto& nxy = scvf.unitOuterNormal();
......
......@@ -72,7 +72,10 @@ class ElasticProblem : public GeomechanicsFVProblem<TypeTag>
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using FVElementGeometry = typename FVGridGeometry::LocalView;
......@@ -127,6 +130,7 @@ public:
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolvars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
GradU gradU = exactGradient(scvf.ipGlobal());
......
......@@ -112,15 +112,17 @@ public:
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param elemVolVars All volume variables for the element
* \param elemFluxVarsCache Flux variables caches for all faces in stencil
* \param scvf The sub control volume face
*
* Negative values mean influx.
* E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
*/
template<class ElementVolumeVariables>
template<class ElementVolumeVariables, class ElementFluxVarsCache>
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFluxVarsCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
......
......@@ -108,15 +108,17 @@ public:
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param elemVolVars All volume variables for the element
* \param elemFluxVarsCache Flux variables caches for all faces in stencil
* \param scvf The sub control volume face
*
* Negative values mean influx.
* E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
*/
template<class ElementVolumeVariables>
template<class ElementVolumeVariables, class ElementFluxVarsCache>
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFluxVarsCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
......
......@@ -199,14 +199,16 @@ public:
* \param element The element for which the Neumann boundary condition is set
* \param fvGeometry The fvGeometry
* \param elemVolVars The element volume variables
* \param elemFluxVarsCache Flux variables caches for all faces in stencil
* \param scvf The boundary sub control volume face
*
* For this method, the \a values variable stores primary variables.
*/
template<class ElementVolumeVariables>
template<class ElementVolumeVariables, class ElementFluxVarsCache>
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFluxVarsCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
......
......@@ -205,14 +205,16 @@ public:
* \param element The element for which the Neumann boundary condition is set
* \param fvGeometry The fvGeometry
* \param elemVolVars The element volume variables
* \param elemFluxVarsCache Flux variables caches for all faces in stencil
* \param scvf The boundary sub control volume face
*
* For this method, the \a values variable stores primary variables.
*/
template<class ElementVolumeVariables>
template<class ElementVolumeVariables, class ElementFluxVarsCache>
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFluxVarsCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
......
......@@ -99,7 +99,9 @@ class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
......@@ -246,12 +248,13 @@ public:
* \param element The element for which the Neumann boundary condition is set
* \param fvGeometry The fvGeometry
* \param elemVolVars The element volume variables
* \param elemFluxVarsCache Flux variables caches for all faces in stencil
* \param scvf The boundary sub control volume face
*
*/
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
......
......@@ -192,14 +192,16 @@ public:
* \param element The element for which the Neumann boundary condition is set
* \param fvGeometry The fvGeometry
* \param elemVolVars The element volume variables
* \param elemFluxVarsCache Flux variables caches for all faces in stencil
* \param scvf The boundary sub control volume face
*
* For this method, the \a values variable stores primary variables.
*/
template<class ElementVolumeVariables>
template<class ElementVolumeVariables, class ElementFluxVarsCache>
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFluxVarsCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
......
......@@ -174,14 +174,16 @@ public:
* \param element The element for which the Neumann boundary condition is set
* \param fvGeometry The fvGeometry
* \param elemVolVars The element volume variables
* \param elemFluxVarsCache Flux variables caches for all faces in stencil
* \param scvf The boundary sub control volume face
*
* For this method, the \a values variable stores primary variables.
*/
template<class ElementVolumeVariables>
template<class ElementVolumeVariables, class ElementFluxVarsCache>
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFluxVarsCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
......
......@@ -83,12 +83,15 @@ class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
using VolumeVariables = typename GridVariables::GridVolumeVariables::VolumeVariables;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
......@@ -203,6 +206,7 @@ public:
* \param element The element for which the Neumann boundary condition is set
* \param fvGeometry The fvGeometry
* \param elemVolVars The element volume variables
* \param elemFluxVarsCache Flux variables caches for all faces in stencil
* \param scvf The boundary sub control volume face
*
* For this method, the \a values variable stores primary variables.
......@@ -210,6 +214,7 @@ public:
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
......
......@@ -151,14 +151,17 @@ double computeGlobalBoundaryMass(const Problem& problem, const SolutionVector& s
for (const auto& element : elements(gg.gridView()))
{
auto fvGeometry = localView(gg);
fvGeometry.bindElement(element);
fvGeometry.bind(element);
auto elemVolVars = localView(gridVars.curGridVolVars());
elemVolVars.bindElement(element, fvGeometry, sol);
elemVolVars.bind(element, fvGeometry, sol);
auto elemFluxVarsCache = localView(gridVars.gridFluxVarsCache());
elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
for (auto&& scvf : scvfs(fvGeometry))
if (scvf.boundary())
mass += problem.neumann(element, fvGeometry, elemVolVars, scvf)[transportEqIdx]
mass += problem.neumann(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf)[transportEqIdx]
* scvf.area() * elemVolVars[scvf.insideScvIdx()].extrusionFactor()
* Problem::FluidSystem::molarMass(transportCompIdx)
* dt;
......
......@@ -221,10 +221,11 @@ public:
* in normal direction of each component. Negative values mean
* influx.
*/
template<class ElementVolumeVariables>
template<class ElementVolumeVariables, class ElementFluxVarsCache>
NeumannFluxes neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolvars,
const ElementFluxVarsCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
NeumannFluxes values(0.0);
......
......@@ -205,10 +205,11 @@ public:
* in normal direction of each component. Negative values mean
* influx.
*/
template<class ElementVolumeVariables>
template<class ElementVolumeVariables, class ElementFluxVarsCache>
NeumannFluxes neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolvars,
const ElementFluxVarsCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
NeumannFluxes values(0.0);
......
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