Commit 35c3b1f6 authored by Timo Koch's avatar Timo Koch
Browse files

[volvars] Reimplement local elementvolvars for 1pcc 1pbox

There are a restriction of the global vector obtained
by a localView call and then bound to an element.
parent 4056c771
......@@ -57,6 +57,7 @@ class DarcysLaw<TypeTag, typename std::enable_if<GET_PROP_VALUE(TypeTag, Discret
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Element = typename GridView::template Codim<0>::Entity;
......@@ -76,6 +77,7 @@ public:
static Scalar flux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const IndexType phaseIdx)
{
......@@ -83,7 +85,8 @@ public:
const auto& faceData = problem.model().fluxVarsCache()[scvf].faceData();
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
const auto extrusionFactor = problem.model().curVolVars(insideScv).extrusionFactor();
const auto& insideVolVars = elemVolVars[insideScv];
const auto extrusionFactor = insideVolVars.extrusionFactor();
const auto K = problem.spatialParams().intrinsicPermeability(insideScv);
// evaluate gradP - rho*g at integration point
......@@ -94,7 +97,7 @@ public:
DimVector gradI;
faceData.jacInvT.mv(faceData.shapeJacobian[scv.index()][0], gradI);
gradI *= problem.model().curVolVars(scv).pressure(phaseIdx);
gradI *= elemVolVars[scv].pressure(phaseIdx);
gradP += gradI;
}
if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
......@@ -103,9 +106,8 @@ public:
DimVector g(problem.gravityAtPos(scvf.center()));
// interpolate the density at the IP
const auto& insideVolVars = problem.model().curVolVars(insideScv);
const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
const auto& outsideVolVars = problem.model().curVolVars(outsideScv);
const auto& outsideVolVars = elemVolVars[outsideScv];
Scalar rho = 0.5*(insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx));
// turn gravity into a force
......
......@@ -18,10 +18,10 @@
*****************************************************************************/
/*!
* \file
* \brief Base class for the volume variables vector
* \brief The local volume variables class
*/
#ifndef DUMUX_DISCRETIZATION_BOX_VOLVARSVECTOR_HH
#define DUMUX_DISCRETIZATION_BOX_VOLVARSVECTOR_HH
#ifndef DUMUX_DISCRETIZATION_BOX_ELEMENT_VOLUMEVARIABLES_HH
#define DUMUX_DISCRETIZATION_BOX_ELEMENT_VOLUMEVARIABLES_HH
#include <dumux/implicit/properties.hh>
......@@ -32,20 +32,20 @@ namespace Dumux
* \ingroup ImplicitModel
* \brief Base class for the volume variables vector
*/
template<class TypeTag, bool useOldSol, bool enableVolVarsCache>
class BoxVolumeVariablesVector
template<class TypeTag, bool enableGlobalVolVarsCache>
class BoxElementVolumeVariables
{};
// specialization in case of storing the volume variables
template<class TypeTag, bool useOldSol>
class BoxVolumeVariablesVector<TypeTag, useOldSol,/*enableVolVarCaching*/true>
template<class TypeTag>
class BoxElementVolumeVariables<TypeTag,/*enableGlobalVolVarCache*/true>
{
friend BoxVolumeVariablesVector<TypeTag, !useOldSol, true>;
friend ImplicitModel<TypeTag>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using GlobalVolumeVariables = typename GET_PROP_TYPE(TypeTag, GlobalVolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using IndexType = typename GridView::IndexSet::IndexType;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
......@@ -55,86 +55,63 @@ class BoxVolumeVariablesVector<TypeTag, useOldSol,/*enableVolVarCaching*/true>
enum{ isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
BoxVolumeVariablesVector<TypeTag, useOldSol, true>& operator= (const BoxVolumeVariablesVector<TypeTag, useOldSol, true>& other) = default;
BoxVolumeVariablesVector<TypeTag, useOldSol, true>& operator= (const BoxVolumeVariablesVector<TypeTag, !useOldSol, true>& other)
{
volVars_ = other.volVars_;
eIdx_ = other.eIdx_;
return *this;
};
public:
void update(Problem& problem, const SolutionVector& sol)
{
problemPtr_ = &problem;
volVars_.resize(problem.gridView().size(0));
for (const auto& element : elements(problem.gridView()))
{
auto fvGeometry = localView(problem.model().globalFvGeometry());
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
(*this)[scv].update(sol[scv.dofIndex()],
problem,
element,
scv);
}
}
}
//! Constructor
BoxElementVolumeVariables(GlobalVolumeVariables& globalVolVars)
: globalVolVarsPtr_(&globalVolVars) {}
const VolumeVariables& operator [](IndexType scvIdx) const
{
return volVars_[eIdx_][scvIdx];
}
{ return globalVolVars().volVars(eIdx_, scvIdx); }
VolumeVariables& operator [](IndexType scvIdx)
{
return volVars_[eIdx_][scvIdx];
}
{ return globalVolVars().volVars(eIdx_, scvIdx); }
const VolumeVariables& operator [](const SubControlVolume& scv) const
{
return volVars_[eIdx_][scv.index()];
}
{ return globalVolVars().volVars(eIdx_, scv.index()); }
VolumeVariables& operator [](const SubControlVolume& scv)
{
return volVars_[eIdx_][scv.index()];
}
{ return globalVolVars().volVars(eIdx_, scv.index()); }
// For compatibility reasons with the case of not storing the vol vars.
// function to be called before assembling an element, preparing the vol vars within the stencil
void bind(const Element& element)
{ bindElement(element); }
void bind(const Element& element,
const FVElementGeometry& fvGeometry,
const SolutionVector& sol)
{
bindElement(element, fvGeometry, sol);
}
// function to prepare the vol vars within the element
void bindElement(const Element& element)
{ eIdx_ = problem_().elementMapper().index(element); }
void bindElement(const Element& element,
const FVElementGeometry& fvGeometry,
const SolutionVector& sol)
{
eIdx_ = globalVolVars().problem_().elementMapper().index(element);
}
private:
const Problem& problem_() const
{ return *problemPtr_; }
//! The global volume variables object we are a restriction of
const GlobalVolumeVariables& globalVolVars() const
{ return *globalVolVarsPtr_; }
const Problem* problemPtr_;
//! The global volume variables object we are a restriction of
GlobalVolumeVariables& globalVolVars()
{ return *globalVolVarsPtr_; }
private:
GlobalVolumeVariables* globalVolVarsPtr_;
IndexType eIdx_;
std::vector<std::vector<VolumeVariables>> volVars_;
};
// Specialization when the current volume variables are not stored
template<class TypeTag>
class BoxVolumeVariablesVector<TypeTag, /*isOldSol*/false, /*enableVolVarCaching*/false>
class BoxElementVolumeVariables<TypeTag, /*enableGlobalVolVarCache*/false>
{
// prev vol vars have to be a friend class in order for the assignment operator to work
friend BoxVolumeVariablesVector<TypeTag, true, false>;
friend ImplicitModel<TypeTag>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using GlobalVolumeVariables = typename GET_PROP_TYPE(TypeTag, GlobalVolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using IndexType = typename GridView::IndexSet::IndexType;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
......@@ -142,144 +119,53 @@ class BoxVolumeVariablesVector<TypeTag, /*isOldSol*/false, /*enableVolVarCaching
static const int dim = GridView::dimension;
using Element = typename GridView::template Codim<0>::Entity;
BoxVolumeVariablesVector& operator= (const BoxVolumeVariablesVector<TypeTag, false, false>& other) = default;
// operator curVolVars = prevVolVars
void operator= (const BoxVolumeVariablesVector<TypeTag, true, false>& other)
{}
public:
void update(Problem& problem, const SolutionVector& sol)
{
problemPtr_ = &problem;
}
//! Constructor
BoxElementVolumeVariables(const GlobalVolumeVariables& globalVolVars)
: globalVolVarsPtr_(&globalVolVars) {}
// specialization for box models, simply forwards to the bindElement method
void bind(const Element& element, const FVElementGeometry& fvGeometry)
void bind(const Element& element,
const FVElementGeometry& fvGeometry,
const SolutionVector& sol)
{
bindElement(element, fvGeometry);
bindElement(element, fvGeometry, sol);
}
// specialization for box models
void bindElement(const Element& element, const FVElementGeometry& fvGeometry)
void bindElement(const Element& element,
const FVElementGeometry& fvGeometry,
const SolutionVector& sol)
{
// resize volume variables to the required size
volVars_.resize(fvGeometry.numScv());
volumeVariables_.resize(fvGeometry.numScv());
for (auto&& scv : scvs(fvGeometry))
{
const auto& sol = problem_().model().curSol()[scv.dofIndex()];
// TODO: INTERFACE SOLVER
// problem_().model().boxInterfaceConditionSolver().updateScvVolVars(element, scv, sol);
volVars_[scv.index()].update(sol, problem_(), element, scv);
// globalVolVars().problem_().model().boxInterfaceConditionSolver().updateScvVolVars(element, scv, sol);
volumeVariables_[scv.index()].update(sol[scv.dofIndex()], globalVolVars().problem_(), element, scv);
}
}
const VolumeVariables& operator [](IndexType scvIdx) const
{
return volVars_[scvIdx];
}
{ return volumeVariables_[scvIdx]; }
VolumeVariables& operator [](IndexType scvIdx)
{
return volVars_[scvIdx];
}
{ return volumeVariables_[scvIdx]; }
const VolumeVariables& operator [](const SubControlVolume& scv) const
{
return volVars_[scv.index()];
}
{ return volumeVariables_[scv.index()]; }
VolumeVariables& operator [](const SubControlVolume& scv)
{
return volVars_[scv.index()];
}
private:
Problem& problem_() const
{ return *problemPtr_;}
Problem* problemPtr_;
std::vector<VolumeVariables> volVars_;
};
// Specialization when the previous volume variables are not stored
template<class TypeTag>
class BoxVolumeVariablesVector<TypeTag, /*isOldSol*/true, /*enableVolVarCaching*/false>
{
// current vol vars have to be a friend class in order for the assignment operator to work
friend BoxVolumeVariablesVector<TypeTag, false, false>;
friend ImplicitModel<TypeTag>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using IndexType = typename GridView::IndexSet::IndexType;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
static const int dim = GridView::dimension;
using Element = typename GridView::template Codim<0>::Entity;
BoxVolumeVariablesVector& operator= (const BoxVolumeVariablesVector<TypeTag, /*isOldSol*/false, /*enableVolVarCaching*/false>& other)
{
volVars_.clear();
problemPtr_ = other.problemPtr_;
return *this;
}
public:
{ return volumeVariables_[scv.index()]; }
void update(const Problem& problem, const SolutionVector& sol)
{
problemPtr_ = &problem;
}
// bind an element
void bindElement(const Element& element, const FVElementGeometry& fvGeometry)
{
// resize volume variables to the required size
volVars_.resize(fvGeometry.numScv());
for (auto&& scv : scvs(fvGeometry))
{
const auto& sol = problem_().model().prevSol()[scv.dofIndex()];
//TODO: INTERFACE SOLVER?
// problem_().model().boxInterfaceConditionSolver().updateScvVolVars(element, scv, sol);
volVars_[scv.index()].update(sol, problem_(), element, scv);
}
}
const VolumeVariables& operator [](IndexType scvIdx) const
{
return volVars_[scvIdx];
}
VolumeVariables& operator [](IndexType scvIdx)
{
return volVars_[scvIdx];
}
const VolumeVariables& operator [](const SubControlVolume& scv) const
{
return volVars_[scv.index()];
}
VolumeVariables& operator [](const SubControlVolume& scv)
{
return volVars_[scv.index()];
}
//! The global volume variables object we are a restriction of
const GlobalVolumeVariables& globalVolVars() const
{ return *globalVolVarsPtr_; }
private:
Problem& problem_() const
{ return *problemPtr_;}
Problem* problemPtr_;
std::vector<VolumeVariables> volVars_;
const GlobalVolumeVariables* globalVolVarsPtr_;
std::vector<VolumeVariables> volumeVariables_;
};
} // end namespace
......
......@@ -51,6 +51,7 @@ class BoxFluxVariablesCacheVector<TypeTag, true>
using Element = typename GridView::template Codim<0>::Entity;
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
public:
......@@ -76,12 +77,16 @@ public:
// Function is called by the BoxLocalJacobian prior to flux calculations on the element.
// We assume the FVGeometries to be bound at this point
void bind(const Element& element, const FVElementGeometry& fvGeometry)
void bind(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars)
{
bindElement(element, fvGeometry);
bindElement(element, fvGeometry, elemVolVars);
}
void bindElement(const Element& element, const FVElementGeometry& fvGeometry)
void bindElement(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars)
{ eIdx_ = problem_().elementMapper().index(element); }
// access operator
......@@ -114,6 +119,7 @@ class BoxFluxVariablesCacheVector<TypeTag, false>
using IndexType = typename GridView::IndexSet::IndexType;
using Element = typename GridView::template Codim<0>::Entity;
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
......@@ -126,12 +132,16 @@ public:
// Function is called by the BoxLocalJacobian prior to flux calculations on the element.
// We assume the FVGeometries to be bound at this point
void bind(const Element& element, const FVElementGeometry& fvGeometry)
void bind(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars)
{
bindElement(element, fvGeometry);
bindElement(element, fvGeometry, elemVolVars);
}
void bindElement(const Element& element, const FVElementGeometry& fvGeometry)
void bindElement(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars)
{
// temporary resizing of the cache
fluxVarsCache_.resize(fvGeometry.numScvf());
......
......@@ -83,13 +83,13 @@ public:
//! Get a sub control volume with a local scv index
const SubControlVolume& scv(IndexType scvIdx) const
{
return globalFvGeometry().scv(scvIdx, eIdx_);
return globalFvGeometry().scvs(eIdx_)[scvIdx];
}
//! Get a sub control volume face with a local scvf index
const SubControlVolumeFace& scvf(IndexType scvfIdx) const
{
return globalFvGeometry().scvf(scvfIdx, eIdx_);
return globalFvGeometry().scvfs(eIdx_)[scvfIdx];
}
//! iterator range for sub control volumes. Iterates over
......
......@@ -164,17 +164,17 @@ public:
static_cast<IndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
// compute the scvf normal unit outer normal
auto normal = geometryHelper.normal(elementGeometry, *scvfGeometry);
auto normal = geometryHelper.normal(elementGeometry, scvfGeometry);
const auto v = elementGeometry.corner(localScvIndices[1]) - elementGeometry.corner(localScvIndices[0]);
const auto s = v*normal;
if (std::signbit(s))
normal *= -1;
scvfs_.emplace_back(std::move(scvfGeometry),
normal,
scvfLocalIdx,
localScvIndices,
false);
scvfs_[eIdx].emplace_back(std::move(scvfGeometry),
normal,
scvfLocalIdx,
localScvIndices,
false);
// increment local counter
scvfLocalIdx++;
......@@ -199,11 +199,11 @@ public:
{static_cast<IndexType>(referenceElement.subEntity(intersection.indexInInside(), 1,
isScvfLocalIdx++, dim))};
scvfs_.emplace_back(std::move(scvfGeometry),
intersection.centerUnitOuterNormal(),
scvfLocalIdx,
localScvIndices,
true);
scvfs_[eIdx].emplace_back(std::move(scvfGeometry),
intersection.centerUnitOuterNormal(),
scvfLocalIdx,
localScvIndices,
true);
// increment local counter
scvfLocalIdx++;
......@@ -228,13 +228,12 @@ public:
{ return feCache_; }
private:
//! Get the local scvs for an element
const std::vector<SubControlVolume>& scvs(IndexType eIdx)
const std::vector<SubControlVolume>& scvs(IndexType eIdx) const
{ return scvs_[eIdx]; }
//! Get the local scvfs for an element
const std::vector<SubControlVolume>& scvfs(IndexType eIdx)
const std::vector<SubControlVolumeFace>& scvfs(IndexType eIdx) const
{ return scvfs_[eIdx]; }
const Problem& problem_() const
......
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief The global volume variables class
*/
#ifndef DUMUX_DISCRETIZATION_BOX_GLOBAL_VOLUMEVARIABLES_HH
#define DUMUX_DISCRETIZATION_BOX_GLOBAL_VOLUMEVARIABLES_HH
#include <dumux/implicit/properties.hh>
#include <dumux/discretization/box/elementvolumevariables.hh>
namespace Dumux
{
/*!
* \ingroup ImplicitModel
* \brief Base class for the volume variables vector
*/
template<class TypeTag, bool enableGlobalVolVarsCache>
class BoxGlobalVolumeVariables
{};
// specialization in case of storing the volume variables
template<class TypeTag>
class BoxGlobalVolumeVariables<TypeTag,/*enableGlobalVolVarCache*/true>
{
friend BoxElementVolumeVariables<TypeTag, true>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using IndexType = typename GridView::IndexSet::IndexType;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
static const int dim = GridView::dimension;
using Element = typename GridView::template Codim<0>::Entity;
enum{ isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
public:
void update(Problem& problem, const SolutionVector& sol)
{
problemPtr_ = &problem;
volumeVariables_.resize(problem.gridView().size(0));
for (const auto& element : elements(problem.gridView()))
{
auto eIdx = problem_().elementMapper().index(element);
auto fvGeometry = localView(problem.model().globalFvGeometry());
fvGeometry.bindElement(element);
volumeVariables_[eIdx].resize(fvGeometry.numScv());
for (auto&& scv : scvs(fvGeometry))
{
volumeVariables_[eIdx][scv.index()].update(sol[scv.dofIndex()],
problem,
element,
scv);
}
}
}
/*!
* \brief Return a local restriction of this global object
* The local object is only functional after calling its bind/bindElement method
* This is a free function that will be found by means of ADL
*/
friend inline ElementVolumeVariables localView(BoxGlobalVolumeVariables& global)