Commit 1391a864 authored by Dennis Gläser's avatar Dennis Gläser Committed by Timo Koch
Browse files

[porousmedium] free volume variables from type tag

Contributions by Timo Koch, Kilian Weishaupt.
parent dd01374b
......@@ -88,7 +88,7 @@ class StaggeredLocalAssembler<TypeTag,
using FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
using FaceSolution = typename GET_PROP_TYPE(TypeTag, StaggeredFaceSolution);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
......@@ -97,7 +97,6 @@ class StaggeredLocalAssembler<TypeTag,
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
static constexpr bool enableGridFluxVarsCache = GET_PROP_VALUE(TypeTag, EnableGridFluxVariablesCache);
static constexpr auto faceOffset = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
public:
......@@ -319,37 +318,40 @@ protected:
const auto& connectivityMap = assembler.fvGridGeometry().connectivityMap();
for(const auto& globalJ : connectivityMap(cellCenterIdx, cellCenterIdx, cellCenterGlobalI))
{
// get the volVars of the element with respect to which we are going to build the derivative
auto&& scvJ = fvGeometry.scv(globalJ);
const auto elementJ = fvGeometry.fvGridGeometry().element(globalJ);
auto& curVolVars = getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scvJ);
VolumeVariables origVolVars(curVolVars);
for(auto pvIdx : priVarIndices_(cellCenterIdx))
{
CellCenterPrimaryVariables priVars(curSol[cellCenterIdx][globalJ]);
const Scalar eps = numericEpsilon(priVars[pvIdx], cellCenterIdx, cellCenterIdx);
priVars[pvIdx] += eps;
auto elemSol = elementSolution<FVElementGeometry>(std::move(priVars));
curVolVars.update(elemSol, problem, elementJ, scvJ);
const auto deflectedResidual = localResidual.evalCellCenter(problem, element, fvGeometry, prevElemVolVars, curElemVolVars,
prevElemFaceVars, curElemFaceVars,
elemBcTypes, elemFluxVarsCache);
auto partialDeriv = (deflectedResidual - ccResidual);
partialDeriv /= eps;
// update the global jacobian matrix with the current partial derivatives
updateGlobalJacobian_(matrix[cellCenterIdx][cellCenterIdx], cellCenterGlobalI, globalJ, pvIdx, partialDeriv);
// restore the original volVars
curVolVars = origVolVars;
}
}
for(const auto& globalJ : connectivityMap(cellCenterIdx, cellCenterIdx, cellCenterGlobalI))
{
// get the volVars of the element with respect to which we are going to build the derivative
auto&& scvJ = fvGeometry.scv(globalJ);
const auto elementJ = fvGeometry.fvGridGeometry().element(globalJ);
auto& curVolVars = getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scvJ);
VolumeVariables origVolVars(curVolVars);
for(auto pvIdx : priVarIndices_(cellCenterIdx))
{
using PrimaryVariables = typename VolumeVariables::PrimaryVariables;
const auto& cellCenterPriVars = curSol[cellCenterIdx][globalJ];
PrimaryVariables priVars = makePriVarsFromCellCenterPriVars<PrimaryVariables>(cellCenterPriVars);
constexpr auto offset = PrimaryVariables::dimension - CellCenterPrimaryVariables::dimension;
const Scalar eps = numericEpsilon(priVars[pvIdx + offset], cellCenterIdx, cellCenterIdx);
priVars[pvIdx + offset] += eps;
auto elemSol = elementSolution<FVElementGeometry>(std::move(priVars));
curVolVars.update(elemSol, problem, elementJ, scvJ);
const auto deflectedResidual = localResidual.evalCellCenter(problem, element, fvGeometry, prevElemVolVars, curElemVolVars,
prevElemFaceVars, curElemFaceVars,
elemBcTypes, elemFluxVarsCache);
auto partialDeriv = (deflectedResidual - ccResidual);
partialDeriv /= eps;
// update the global jacobian matrix with the current partial derivatives
updateGlobalJacobian_(matrix[cellCenterIdx][cellCenterIdx], cellCenterGlobalI, globalJ, pvIdx, partialDeriv);
// restore the original volVars
curVolVars = origVolVars;
}
}
}
/*!
......@@ -369,17 +371,17 @@ protected:
JacobianMatrix& matrix,
const NumCellCenterEqVector& ccResidual)
{
const auto& problem = assembler.problem();
auto& localResidual = assembler.localResidual();
auto& gridVariables = assembler.gridVariables();
static constexpr auto cellCenterIdx = FVGridGeometry::cellCenterIdx();
static constexpr auto faceIdx = FVGridGeometry::faceIdx();
const auto& problem = assembler.problem();
auto& localResidual = assembler.localResidual();
auto& gridVariables = assembler.gridVariables();
static constexpr auto cellCenterIdx = FVGridGeometry::cellCenterIdx();
static constexpr auto faceIdx = FVGridGeometry::faceIdx();
// build derivatives with for cell center dofs w.r.t. cell center dofs
const auto cellCenterGlobalI = assembler.fvGridGeometry().elementMapper().index(element);
// build derivatives with for cell center dofs w.r.t. cell center dofs
const auto cellCenterGlobalI = assembler.fvGridGeometry().elementMapper().index(element);
for(const auto& scvfJ : scvfs(fvGeometry))
{
for(const auto& scvfJ : scvfs(fvGeometry))
{
const auto globalJ = scvfJ.dofIndex();
// get the faceVars of the face with respect to which we are going to build the derivative
......@@ -403,7 +405,7 @@ protected:
partialDeriv /= eps;
// update the global jacobian matrix with the current partial derivatives
updateGlobalJacobian_(matrix[cellCenterIdx][faceIdx], cellCenterGlobalI, globalJ, pvIdx - faceOffset, partialDeriv);
updateGlobalJacobian_(matrix[cellCenterIdx][faceIdx], cellCenterGlobalI, globalJ, pvIdx, partialDeriv);
// restore the original faceVars
faceVars = origFaceVars;
......@@ -451,17 +453,20 @@ protected:
for(auto pvIdx : priVarIndices_(cellCenterIdx))
{
CellCenterPrimaryVariables priVars(curSol[cellCenterIdx][globalJ]);
using PrimaryVariables = typename VolumeVariables::PrimaryVariables;
const auto& cellCenterPriVars = curSol[cellCenterIdx][globalJ];
PrimaryVariables priVars = makePriVarsFromCellCenterPriVars<PrimaryVariables>(cellCenterPriVars);
const Scalar eps = numericEpsilon(priVars[pvIdx], faceIdx, cellCenterIdx);
priVars[pvIdx] += eps;
constexpr auto offset = PrimaryVariables::dimension - CellCenterPrimaryVariables::dimension;
const Scalar eps = numericEpsilon(priVars[pvIdx + offset], faceIdx, cellCenterIdx);
priVars[pvIdx + offset] += eps;
auto elemSol = elementSolution<FVElementGeometry>(std::move(priVars));
curVolVars.update(elemSol, problem, elementJ, scvJ);
const auto deflectedResidual = localResidual.evalFace(problem, element, fvGeometry, scvf,
prevElemVolVars, curElemVolVars,
prevElemFaceVars, curElemFaceVars,
elemBcTypes, elemFluxVarsCache);
prevElemVolVars, curElemVolVars,
prevElemFaceVars, curElemFaceVars,
elemBcTypes, elemFluxVarsCache);
auto partialDeriv = (deflectedResidual - cachedResidual[scvf.localFaceIdx()]);
partialDeriv /= eps;
......@@ -529,7 +534,7 @@ protected:
partialDeriv /= eps;
// update the global jacobian matrix with the current partial derivatives
updateGlobalJacobian_(matrix[faceIdx][faceIdx], faceGlobalI, globalJ, pvIdx - faceOffset, partialDeriv);
updateGlobalJacobian_(matrix[faceIdx][faceIdx], faceGlobalI, globalJ, pvIdx, partialDeriv);
// restore the original faceVars
faceVars = origFaceVars;
......@@ -601,7 +606,6 @@ protected:
static auto priVarIndices_(typename FVGridGeometry::DofTypeIndices::CellCenterIdx)
{
constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
return Dune::range(0, numEqCellCenter);
#else
......@@ -613,12 +617,11 @@ protected:
//! Specialization for face dofs.
static auto priVarIndices_(typename FVGridGeometry::DofTypeIndices::FaceIdx)
{
constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
constexpr auto numEqFace = GET_PROP_VALUE(TypeTag, NumEqFace);
#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
return Dune::range(numEqCellCenter, numEqCellCenter + numEqFace);
return Dune::range(0, numEqFace);
#else
return IntRange(numEqCellCenter, numEqCellCenter + numEqFace);
return IntRange(0, numEqFace);
#endif
}
......
......@@ -42,7 +42,6 @@ NEW_PROP_TAG(ModelParameterGroup); //!< Property which defines the group that
NEW_PROP_TAG(ModelDefaultParameters); //!< Property which defines the group that is queried for parameters by default
NEW_PROP_TAG(GridCreator); //!< Property which provides a GridCreator (manages grids)
NEW_PROP_TAG(Grid); //!< The DUNE grid type
NEW_PROP_TAG(Indices); //!< Enumerations for the numeric model
NEW_PROP_TAG(PrimaryVariables); //!< A vector of primary variables
NEW_PROP_TAG(NumEqVector); //!< A vector of size number equations that can be used for Neumann fluxes, sources, residuals, ...
NEW_PROP_TAG(GridView); //!< The type of the grid view according to the grid type
......
......@@ -135,8 +135,12 @@ public:
const SubControlVolume& scv,
const PrimaryVariables& initSol) const
{
constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
constexpr auto numEq = GET_PROP_TYPE(TypeTag, ModelTraits)::numEq();
constexpr auto offset = numEq - numEqCellCenter;
for(auto&& i : priVarIndices_(cellCenterIdx))
sol[cellCenterIdx][scv.dofIndex()][i] = initSol[i];
sol[cellCenterIdx][scv.dofIndex()][i] = initSol[i + offset];
}
//! Applys the initial face solution
......@@ -162,7 +166,6 @@ protected:
static auto priVarIndices_(typename FVGridGeometry::DofTypeIndices::CellCenterIdx)
{
constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
return Dune::range(0, numEqCellCenter);
#else
......@@ -174,12 +177,11 @@ protected:
//! Specialization for face dofs.
static auto priVarIndices_(typename FVGridGeometry::DofTypeIndices::FaceIdx)
{
constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
constexpr auto numEq = GET_PROP_TYPE(TypeTag, ModelTraits)::numEq();
constexpr auto numEqFace = GET_PROP_VALUE(TypeTag, NumEqFace);
#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
return Dune::range(numEqCellCenter, numEq);
return Dune::range(0, numEqFace);
#else
return IntRange(numEqCellCenter, numEq);
return IntRange(0, numEqFace);
#endif
}
......
......@@ -54,16 +54,17 @@ class FicksLawImplementation<TypeTag, DiscretizationMethod::box>
using FluxVarCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using BalanceEqOpts = typename GET_PROP_TYPE(TypeTag, BalanceEqOpts);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
using Element = typename GridView::template Codim<0>::Entity;
using IndexType = typename GridView::IndexSet::IndexType;
using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
using Indices = typename ModelTraits::Indices;
enum { dim = GridView::dimension} ;
enum { dimWorld = GridView::dimensionworld} ;
enum
{
numPhases = GET_PROP_TYPE(TypeTag, ModelTraits)::numPhases(),
numComponents = GET_PROP_TYPE(TypeTag, ModelTraits)::numComponents()
numPhases = ModelTraits::numPhases(),
numComponents = ModelTraits::numComponents()
};
using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
......
......@@ -57,7 +57,7 @@ class MaxwellStefansLawImplementation<TypeTag, DiscretizationMethod::box >
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, GridFluxVariablesCache)::LocalView;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using Element = typename GridView::template Codim<0>::Entity;
enum { dim = GridView::dimension} ;
......
......@@ -58,12 +58,14 @@ class FicksLawImplementation<TypeTag, DiscretizationMethod::cctpfa>
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using BalanceEqOpts = typename GET_PROP_TYPE(TypeTag, BalanceEqOpts);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
using Indices = typename ModelTraits::Indices;
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
static const int numPhases = GET_PROP_TYPE(TypeTag, ModelTraits)::numPhases();
static const int numComponents = GET_PROP_TYPE(TypeTag, ModelTraits)::numComponents();
static const int numPhases = ModelTraits::numPhases();
static const int numComponents = ModelTraits::numComponents();
using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
......
......@@ -58,7 +58,7 @@ class MaxwellStefansLawImplementation<TypeTag, DiscretizationMethod::cctpfa >
using Element = typename GridView::template Codim<0>::Entity;
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, GridFluxVariablesCache)::LocalView;
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
static const int dim = GridView::dimension;
......
......@@ -30,6 +30,26 @@
namespace Dumux {
/*!
* \ingroup StaggeredDiscretization
* \brief Helper function to creater a PrimaryVariables object from CellCenterPrimaryVariables
* \tparam PrimaryVariables The type of the desired primary variables object
* \tparam CellCenterPrimaryVariables The type of the cell center (input) primary variables object
* \param CellCenterPrimaryVariables The cell center (input) primary variables object
*/
template<class PrimaryVariables, class CellCenterPrimaryVariables>
PrimaryVariables makePriVarsFromCellCenterPriVars(const CellCenterPrimaryVariables& cellCenterPriVars)
{
static_assert(int(PrimaryVariables::dimension) > int(CellCenterPrimaryVariables::dimension),
"PrimaryVariables' size must be greater than the one of CellCenterPrimaryVariables");
PrimaryVariables priVars(0.0);
constexpr auto offset = PrimaryVariables::dimension - CellCenterPrimaryVariables::dimension;
for (std::size_t i = 0; i < cellCenterPriVars.size(); ++i)
priVars[i + offset] = cellCenterPriVars[i];
return priVars;
}
template<class PrimaryVariables>
using StaggeredElementSolution = Dune::BlockVector<PrimaryVariables>;
......
......@@ -47,7 +47,6 @@ public:
resetEq(eqIdx);
}
/*!
* \brief Reset the boundary types for one equation.
*/
......
......@@ -54,7 +54,7 @@ class FicksLawImplementation<TypeTag, DiscretizationMethod::staggered >
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
static constexpr int dim = GridView::dimension;
......@@ -106,7 +106,7 @@ public:
if(scvf.boundary())
{
const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
if(bcTypes.isOutflow(eqIdx) && eqIdx != pressureIdx)
if(bcTypes.isOutflow(eqIdx))
return flux;
}
......
......@@ -51,7 +51,7 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethod::staggered >
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
using Element = typename GridView::template Codim<0>::Entity;
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
enum { energyBalanceIdx = Indices::energyBalanceIdx };
......
......@@ -54,7 +54,7 @@ class MaxwellStefansLawImplementation<TypeTag, DiscretizationMethod::staggered >
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
static const int dim = GridView::dimension;
......
......@@ -97,10 +97,7 @@ public:
//! Boundary types at a single degree of freedom
SET_PROP(StaggeredFreeFlowModel, BoundaryTypes)
{
private:
static constexpr auto size = GET_PROP_VALUE(TypeTag, NumEqCellCenter) + GET_PROP_VALUE(TypeTag, NumEqFace);
public:
using type = StaggeredFreeFlowBoundaryTypes<size>;
using type = StaggeredFreeFlowBoundaryTypes<GET_PROP_TYPE(TypeTag, ModelTraits)::numEq()>;
};
//! The velocity output
......
......@@ -49,6 +49,7 @@ class StaggeredGridVolumeVariables<Traits, /*cachingEnabled*/true>
{
using ThisType = StaggeredGridVolumeVariables<Traits, true>;
using Problem = typename Traits::Problem;
using PrimaryVariables = typename Traits::VolumeVariables::PrimaryVariables;
public:
//! export the type of the indices TODO: get them out of the volvars
......@@ -81,7 +82,10 @@ public:
for (auto&& scv : scvs(fvGeometry))
{
CellCenterPrimaryVariables priVars = sol[scv.dofIndex()];
// construct a privars object from the cell center solution vector
const auto& cellCenterPriVars = sol[scv.dofIndex()];
PrimaryVariables priVars = makePriVarsFromCellCenterPriVars<PrimaryVariables>(cellCenterPriVars);
auto elemSol = elementSolution<typename FVGridGeometry::LocalView>(std::move(priVars));
volumeVariables_[scv.dofIndex()].update(elemSol, problem(), element, scv);
}
......@@ -97,14 +101,14 @@ public:
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideScv = fvGeometry.scv(insideScvIdx);
CellCenterPrimaryVariables boundaryPriVars(0.0);
for(int eqIdx = 0; eqIdx < CellCenterPrimaryVariables::dimension; ++eqIdx)
PrimaryVariables boundaryPriVars(0.0);
constexpr auto offset = PrimaryVariables::dimension - CellCenterPrimaryVariables::dimension;
for(int eqIdx = offset; eqIdx < PrimaryVariables::dimension; ++eqIdx)
{
if(bcTypes.isDirichlet(eqIdx) || bcTypes.isDirichletCell(eqIdx))
boundaryPriVars[eqIdx] = problem().dirichlet(element, scvf)[eqIdx];
else if(bcTypes.isNeumann(eqIdx) || bcTypes.isOutflow(eqIdx) || bcTypes.isSymmetry())
boundaryPriVars[eqIdx] = sol[scvf.insideScvIdx()][eqIdx];
boundaryPriVars[eqIdx] = sol[scvf.insideScvIdx()][eqIdx - offset];
//TODO: this assumes a zero-gradient for e.g. the pressure on the boundary
// could be made more general by allowing a non-zero-gradient, provided in problem file
else
......
......@@ -81,7 +81,7 @@ private:
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices); // TODO extract indices from volumevariables
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices; // TODO extract indices from volumevariables
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
static constexpr auto enableCache = GET_PROP_VALUE(TypeTag, EnableGridVolumeVariablesCache);
......@@ -135,6 +135,7 @@ SET_TYPE_PROP(StaggeredModel,
CellCenterPrimaryVariables,
Dune::FieldVector<typename GET_PROP_TYPE(TypeTag, Scalar),
GET_PROP_VALUE(TypeTag, NumEqCellCenter)>);
//! The face primary variables
SET_TYPE_PROP(StaggeredModel,
FacePrimaryVariables,
......@@ -142,14 +143,7 @@ SET_TYPE_PROP(StaggeredModel,
GET_PROP_VALUE(TypeTag, NumEqFace)>);
//! Boundary types at a single degree of freedom
SET_PROP(StaggeredModel, BoundaryTypes)
{
private:
static constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
static constexpr auto numEqFace = GET_PROP_VALUE(TypeTag, NumEqFace);
public:
using type = BoundaryTypes<numEqCellCenter + numEqFace>;
};
SET_TYPE_PROP(StaggeredModel, BoundaryTypes, BoundaryTypes<GET_PROP_TYPE(TypeTag, ModelTraits)::numEq()>);
//! Set one or different base epsilons for the calculations of the localJacobian's derivatives
SET_PROP(StaggeredModel, BaseEpsilon)
......
......@@ -24,42 +24,33 @@
#ifndef DUMUX_NAVIERSTOKES_INDICES_HH
#define DUMUX_NAVIERSTOKES_INDICES_HH
#include <dumux/common/properties.hh>
namespace Dumux {
// \{
/*!
* \ingroup NavierStokesModel
* \brief The common indices for the isothermal Navier-Stokes model.
*
* \tparam dimension The dimension of the problem
* \tparam numEquations the number of model equations
* \tparam PVOffset The first index in a primary variable vector.
*/
template <int dimension, int numEquations, int PVOffset = 0>
template <int dimension>
struct NavierStokesIndices
{
static constexpr int dimXIdx = 0; //!< Index of the x-component of a vector of size dim
static constexpr int dimYIdx = 1; //!< Index of the y-component of a vector of size dim
static constexpr int dimZIdx = 2; //!< Index of the z-component of a vector of size dim
static constexpr int totalMassBalanceIdx = PVOffset + 0; //!< Index of the total mass balance equation
static constexpr int pressureIdx = totalMassBalanceIdx; //!< Index of the pressure in a solution vector
static constexpr int conti0EqIdx = dimension; //!< Index of the total mass balance equation
static constexpr int pressureIdx = conti0EqIdx; //!< Index of the pressure in a solution vector
static constexpr auto dim = dimension;
static constexpr auto numEq = numEquations;
static constexpr auto momentumBalanceOffset = numEquations - dim;
static constexpr int momentumBalanceIdx = PVOffset + momentumBalanceOffset; //!< Index of the momentum balance equation
static constexpr int momentumXBalanceIdx = momentumBalanceIdx; //!< Index of the momentum balance equation
static constexpr int momentumYBalanceIdx = momentumBalanceIdx + 1; //!< Index of the momentum balance equation
static constexpr int momentumZBalanceIdx = momentumBalanceIdx + 2; //!< Index of the momentum balance equation
static constexpr int momentumXBalanceIdx = 0; //!< Index of the momentum balance equation
static constexpr int momentumYBalanceIdx = 1; //!< Index of the momentum balance equation
static constexpr int momentumZBalanceIdx = 2; //!< Index of the momentum balance equation
static constexpr int velocityXIdx = momentumBalanceIdx; //!< Index of the velocity in a solution vector
static constexpr int velocityYIdx = momentumBalanceIdx + 1; //!< Index of the velocity in a solution vector
static constexpr int velocityZIdx = momentumBalanceIdx + 2; //!< Index of the velocity in a solution vector
static constexpr int velocityXIdx = 0; //!< Index of the velocity in a solution vector
static constexpr int velocityYIdx = 1; //!< Index of the velocity in a solution vector
static constexpr int velocityZIdx = 2; //!< Index of the velocity in a solution vector
/*!
* \brief Index of the velocity in a solution vector given a certain direction.
......@@ -68,11 +59,10 @@ struct NavierStokesIndices
*/
static constexpr int velocity(int dirIdx)
{
return dirIdx + momentumBalanceIdx;
return dirIdx;
}
};
// \}
} // end namespace
} // end namespace Dumux
#endif
......@@ -65,19 +65,21 @@
#include <dumux/discretization/methods.hh>
#include <dumux/discretization/fourierslaw.hh>
namespace Dumux
{
namespace Dumux {
/*!
* \ingroup NavierStokesModel
* \brief Traits for the Navier-Stokes model
*/
template<int dim>
template<int dimension>
struct NavierStokesModelTraits
{
//! The dimension of the model
static constexpr int dim() { return dimension; }
//! There are as many momentum balance equations as dimensions
//! and one mass balance equation.
static constexpr int numEq() { return dim+1; }
static constexpr int numEq() { return dimension+1; }
//! The number of phases is 1
static constexpr int numPhases() { return 1; }
......@@ -93,6 +95,9 @@ struct NavierStokesModelTraits
//! The model is isothermal
static constexpr bool enableEnergyBalance() { return false; }
//! the indices
using Indices = NavierStokesIndices<dim()>;
};
// \{
......@@ -154,16 +159,6 @@ SET_TYPE_PROP(NavierStokes, FluxVariables, NavierStokesFluxVariables<TypeTag>);
//! The flux variables cache class, by default the one for free flow
SET_TYPE_PROP(NavierStokes, FluxVariablesCache, FreeFlowFluxVariablesCache<TypeTag>);
//! The indices required by the isothermal single-phase model
SET_PROP(NavierStokes, Indices)
{
private:
static constexpr int numEq = GET_PROP_TYPE(TypeTag, ModelTraits)::numEq();
static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension;
public:
using type = NavierStokesIndices<dim, numEq>;
};
//! The specific vtk output fields
SET_PROP(NavierStokes, VtkOutputFields)
{
......@@ -187,17 +182,6 @@ public:
using type = NavierStokesNIModelTraits<IsothermalTraits>;
};
//! The indices required by the non-isothermal single-phase model
SET_PROP(NavierStokesNI, Indices)
{
private:
static constexpr int numEq = GET_PROP_TYPE(TypeTag, ModelTraits)::numEq();
static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension;
using IsothermalIndices = NavierStokesIndices<dim, numEq>;
public:
using type = NavierStokesNonIsothermalIndices<dim, numEq, IsothermalIndices>;
};