Skip to content
GitLab
Commit 0934d808 authored by Thomas Fetzer's avatar Thomas Fetzer
Browse files

[navierstokes] Improve indices 

* better distinguish between total mass balance and component mass balances
parent 856d107c
Hide whitespace changes
Inline Side-by-side
dumux/freeflow/navierstokes/indices.hh
View file @ 0934d808
......@@ -45,9 +45,8 @@ struct NavierStokesIndices
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 massBalanceIdx = PVOffset + 0; //!< Index of the mass balance equation
static constexpr int conti0EqIdx = massBalanceIdx; //!< Index of the mass balance equation
static constexpr int pressureIdx = massBalanceIdx; //!< Index of the pressure in a solution vector
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 auto dim = dimension;
static constexpr auto numEq = numEquations;
......
dumux/freeflow/navierstokes/staggered/fluxoverplane.hh
View file @ 0934d808
......@@ -417,7 +417,7 @@ private:
}();
CellCenterPrimaryVariables tmp(0.0);
tmp[Indices::conti0EqIdx] = cumulativeFlux / avgDensity;
tmp[Indices::totalMassBalanceIdx] = cumulativeFlux / avgDensity;
return tmp;
};
......@@ -438,7 +438,9 @@ private:
const auto& scvf,
const auto& elemFluxVarsCache)
{
const Scalar massFlux = localResidual_.computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache)[Indices::conti0EqIdx];
const Scalar totalMassFlux = localResidual_.computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars,
elemFaceVars, scvf, elemFluxVarsCache)
[Indices::totalMassBalanceIdx];
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
......@@ -446,7 +448,7 @@ private:
const auto avgDensity = 0.5*insideVolVars.density() + 0.5*outsideVolVars.density();
CellCenterPrimaryVariables tmp(0.0);
tmp[Indices::conti0EqIdx] = massFlux / avgDensity;
tmp[Indices::totalMassBalanceIdx] = totalMassFlux / avgDensity;
return tmp;
};
......
dumux/freeflow/navierstokes/staggered/fluxvariables.hh
View file @ 0934d808
......@@ -132,14 +132,14 @@ public:
// Check if we are on an outflow boundary.
const bool isOutflow = scvf.boundary() ?
problem.boundaryTypesAtPos(scvf.center()).isOutflow(Indices::massBalanceIdx)
problem.boundaryTypesAtPos(scvf.center()).isOutflow(Indices::totalMassBalanceIdx)
: false;
// Call the generic flux function.
const Scalar flux = advectiveFluxForCellCenter(elemVolVars, elemFaceVars, scvf, upwindTerm, isOutflow);
CellCenterPrimaryVariables result(0.0);
result[Indices::massBalanceIdx] = flux;
result[Indices::totalMassBalanceIdx] = flux;
return result;
}
......
dumux/freeflow/navierstokes/staggered/localresidual.hh
View file @ 0934d808
......@@ -73,7 +73,7 @@ class NavierStokesResidualImpl<TypeTag, DiscretizationMethod::staggered>
enum {
pressureIdx = Indices::pressureIdx,
massBalanceIdx = Indices::massBalanceIdx,
totalMassBalanceIdx = Indices::totalMassBalanceIdx,
momentumBalanceIdx = Indices::momentumBalanceIdx
};
......@@ -132,7 +132,7 @@ public:
const VolumeVariables& volVars) const
{
CellCenterPrimaryVariables storage;
storage[massBalanceIdx] = volVars.density();
storage[totalMassBalanceIdx] = volVars.density();
computeStorageForCellCenterNonIsothermal_(std::integral_constant<bool, ModelTraits::enableEnergyBalance() >(),
problem, scv, volVars, storage);
......@@ -251,7 +251,7 @@ protected:
const BoundaryTypes& bcTypes) const
{
// set a fixed pressure for cells adjacent to a wall
if(bcTypes.isDirichletCell(massBalanceIdx))
if(bcTypes.isDirichletCell(totalMassBalanceIdx))
{
const auto& insideVolVars = elemVolVars[insideScv];
residual[pressureIdx] = insideVolVars.pressure() - problem.dirichletAtPos(insideScv.center())[pressureIdx];
......@@ -297,7 +297,7 @@ protected:
// outflow condition for the momentum balance equation
if(bcTypes.isOutflow(momentumBalanceIdx))
{
if(bcTypes.isDirichlet(massBalanceIdx))
if(bcTypes.isDirichlet(totalMassBalanceIdx))
residual += computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elementFaceVars, elemFluxVarsCache);
else
DUNE_THROW(Dune::InvalidStateException, "Face at " << scvf.center() << " has an outflow BC for the momentum balance but no Dirichlet BC for the pressure!");
......
dumux/freeflow/navierstokesnc/indices.hh
View file @ 0934d808
......@@ -45,12 +45,13 @@ private:
using ParentType = NavierStokesIndices<dimension, numEquations, PVOffset>;
public:
static constexpr int phaseIdx = thePhaseIdx; //!< The phase index
static constexpr int mainCompIdx = phaseIdx; //!< The index of the main component
//! The index of the component whose mass balance will be replaced by the total one
static constexpr int replaceCompEqIdx = theReplaceCompEqIdx;
static constexpr int totalMassBalanceIdx = replaceCompEqIdx; //!< Index of the total mass balance equation
static constexpr int conti0EqIdx = PVOffset; //!< The base index of the transport equations
};
// \}
......
dumux/freeflow/rans/problem.hh
View file @ 0934d808
......@@ -73,11 +73,6 @@ class RANSProblem : public NavierStokesProblem<TypeTag>
using DimVector = Dune::FieldVector<Scalar, dim>;
using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
enum {
massBalanceIdx = Indices::massBalanceIdx,
momentumBalanceIdx = Indices::momentumBalanceIdx
};
using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
typename DofTypeIndices::CellCenterIdx cellCenterIdx;
typename DofTypeIndices::FaceIdx faceIdx;
......@@ -236,7 +231,7 @@ public:
for (auto&& scvf : scvfs(fvGeometry))
{
const int dofIdxFace = scvf.dofIndex();
const auto numericalSolutionFace = curSol[faceIdx][dofIdxFace][momentumBalanceIdx];
const auto numericalSolutionFace = curSol[faceIdx][dofIdxFace][Indices::momentumBalanceIdx];
velocityTemp[scvf.directionIndex()] += numericalSolutionFace;
}
for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
......
dumux/freeflow/rans/zeroeq/problem.hh
View file @ 0934d808
......@@ -69,11 +69,6 @@ class ZeroEqProblem : public RANSProblem<TypeTag>
using DimVector = Dune::FieldVector<Scalar, dim>;
using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
enum {
massBalanceIdx = Indices::massBalanceIdx,
momentumBalanceIdx = Indices::momentumBalanceIdx
};
using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
typename DofTypeIndices::CellCenterIdx cellCenterIdx;
typename DofTypeIndices::FaceIdx faceIdx;
......
test/freeflow/navierstokes/angelitestproblem.hh
View file @ 0934d808
......@@ -84,7 +84,7 @@ class AngeliTestProblem : public NavierStokesProblem<TypeTag>
dimWorld = GridView::dimensionworld
};
enum {
massBalanceIdx = Indices::massBalanceIdx,
totalMassBalanceIdx = Indices::totalMassBalanceIdx,
momentumBalanceIdx = Indices::momentumBalanceIdx,
pressureIdx = Indices::pressureIdx,
velocityXIdx = Indices::velocityXIdx,
......@@ -186,7 +186,7 @@ public:
values.setDirichlet(momentumBalanceIdx);
// set a fixed pressure in one cell
values.setDirichletCell(massBalanceIdx);
values.setDirichletCell(totalMassBalanceIdx);
return values;
}
......
test/freeflow/navierstokes/channeltestproblem.hh
View file @ 0934d808
......@@ -92,7 +92,7 @@ class ChannelTestProblem : public NavierStokesProblem<TypeTag>
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum { dimWorld = GridView::dimensionworld };
enum {
massBalanceIdx = Indices::massBalanceIdx,
totalMassBalanceIdx = Indices::totalMassBalanceIdx,
momentumBalanceIdx = Indices::momentumBalanceIdx,
pressureIdx = Indices::pressureIdx,
#if NONISOTHERMAL
......@@ -178,11 +178,11 @@ public:
// set a fixed pressure in one cell
if (isOutlet(globalPos))
{
values.setDirichlet(massBalanceIdx);
values.setDirichlet(totalMassBalanceIdx);
values.setOutflow(momentumBalanceIdx);
}
else
values.setOutflow(massBalanceIdx);
values.setOutflow(totalMassBalanceIdx);
return values;
}
......
test/freeflow/navierstokes/closedsystemtestproblem.hh
View file @ 0934d808
......@@ -77,7 +77,7 @@ class ClosedSystemTestProblem : public NavierStokesProblem<TypeTag>
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum { dimWorld = GridView::dimensionworld };
enum {
massBalanceIdx = Indices::massBalanceIdx,
totalMassBalanceIdx = Indices::totalMassBalanceIdx,
momentumBalanceIdx = Indices::momentumBalanceIdx,
pressureIdx = Indices::pressureIdx,
velocityXIdx = Indices::velocityXIdx,
......@@ -156,9 +156,9 @@ public:
// set a fixed pressure in one cell
if (isLowerLeftCell_(globalPos))
values.setDirichletCell(massBalanceIdx);
values.setDirichletCell(totalMassBalanceIdx);
else
values.setNeumann(massBalanceIdx);
values.setNeumann(totalMassBalanceIdx);
return values;
}
......
test/freeflow/navierstokes/doneatestproblem.hh
View file @ 0934d808
......@@ -87,7 +87,7 @@ class DoneaTestProblem : public NavierStokesProblem<TypeTag>
dimWorld = GridView::dimensionworld
};
enum {
massBalanceIdx = Indices::massBalanceIdx,
totalMassBalanceIdx = Indices::totalMassBalanceIdx,
momentumBalanceIdx = Indices::momentumBalanceIdx,
momentumXBalanceIdx = Indices::momentumXBalanceIdx,
momentumYBalanceIdx = Indices::momentumYBalanceIdx,
......@@ -205,7 +205,7 @@ public:
// set Dirichlet values for the velocity and pressure everywhere
values.setDirichlet(momentumBalanceIdx);
values.setDirichletCell(massBalanceIdx);
values.setDirichletCell(totalMassBalanceIdx);
return values;
}
......
test/freeflow/navierstokes/kovasznaytestproblem.hh
View file @ 0934d808
......@@ -82,7 +82,7 @@ class KovasznayTestProblem : public NavierStokesProblem<TypeTag>
dimWorld = GridView::dimensionworld
};
enum {
massBalanceIdx = Indices::massBalanceIdx,
totalMassBalanceIdx = Indices::totalMassBalanceIdx,
momentumBalanceIdx = Indices::momentumBalanceIdx,
pressureIdx = Indices::pressureIdx,
velocityXIdx = Indices::velocityXIdx,
......@@ -190,9 +190,9 @@ public:
// set a fixed pressure in one cell
if (isLowerLeftCell_(globalPos))
values.setDirichletCell(massBalanceIdx);
values.setDirichletCell(totalMassBalanceIdx);
else
values.setOutflow(massBalanceIdx);
values.setOutflow(totalMassBalanceIdx);
return values;
}
......
test/freeflow/navierstokesnc/channeltestproblem.hh
View file @ 0934d808
......@@ -54,15 +54,10 @@ NEW_PROP_TAG(FluidSystem);
SET_TYPE_PROP(ChannelNCTestTypeTag, FluidSystem,
FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)/*, SimpleH2O<typename GET_PROP_TYPE(TypeTag, Scalar)>, true*/>);
SET_PROP(ChannelNCTestTypeTag, PhaseIdx)
{
private:
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
public:
static constexpr int value = FluidSystem::wPhaseIdx;
};
SET_INT_PROP(ChannelNCTestTypeTag, PhaseIdx,
GET_PROP_TYPE(TypeTag, FluidSystem)::wPhaseIdx);
SET_INT_PROP(ChannelNCTestTypeTag, ReplaceCompEqIdx, 0);
SET_INT_PROP(ChannelNCTestTypeTag, ReplaceCompEqIdx, GET_PROP_VALUE(TypeTag, PhaseIdx));
// Set the grid type
SET_TYPE_PROP(ChannelNCTestTypeTag, Grid, Dune::YaspGrid<2>);
......@@ -102,8 +97,8 @@ class ChannelNCTestProblem : public NavierStokesProblem<TypeTag>
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum { dimWorld = GridView::dimensionworld };
enum {
massBalanceIdx = Indices::massBalanceIdx,
transportEqIdx = 1,
totalMassBalanceIdx = Indices::totalMassBalanceIdx,
transportEqIdx = 1-GET_PROP_VALUE(TypeTag, PhaseIdx),
momentumBalanceIdx = Indices::momentumBalanceIdx,
pressureIdx = Indices::pressureIdx,
velocityXIdx = Indices::velocityXIdx,
......@@ -112,7 +107,7 @@ class ChannelNCTestProblem : public NavierStokesProblem<TypeTag>
temperatureIdx = Indices::temperatureIdx,
energyBalanceIdx = Indices::energyBalanceIdx,
#endif
transportCompIdx = 1/*FluidSystem::wCompIdx*/
transportCompIdx = 1-GET_PROP_VALUE(TypeTag, PhaseIdx)
};
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
......@@ -142,7 +137,6 @@ public:
*/
// \{
bool shouldWriteRestartFile() const
{
return false;
......@@ -165,6 +159,7 @@ public:
{
return NumEqVector(0.0);
}
// \}
/*!
* \name Boundary conditions
......@@ -184,7 +179,7 @@ public:
if(isInlet(globalPos))
{
values.setDirichlet(momentumBalanceIdx);
values.setOutflow(massBalanceIdx);
values.setOutflow(totalMassBalanceIdx);
values.setDirichlet(transportEqIdx);
#if NONISOTHERMAL
values.setDirichlet(energyBalanceIdx);
......@@ -193,18 +188,17 @@ public:
else if(isOutlet(globalPos))
{
values.setOutflow(momentumBalanceIdx);
values.setDirichlet(massBalanceIdx);
values.setDirichlet(totalMassBalanceIdx);
values.setOutflow(transportEqIdx);
#if NONISOTHERMAL
values.setOutflow(energyBalanceIdx);
#endif
}
else
{
// set Dirichlet values for the velocity everywhere
values.setDirichlet(momentumBalanceIdx);
values.setOutflow(massBalanceIdx);
values.setOutflow(totalMassBalanceIdx);
values.setOutflow(transportEqIdx);
#if NONISOTHERMAL
values.setOutflow(energyBalanceIdx);
......
test/freeflow/navierstokesnc/densityflowproblem.hh
View file @ 0934d808
......@@ -49,15 +49,10 @@ NEW_PROP_TAG(FluidSystem);
SET_TYPE_PROP(DensityDrivenFlowTypeTag, FluidSystem,
FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)/*, SimpleH2O<typename GET_PROP_TYPE(TypeTag, Scalar)>, false*/>);
SET_PROP(DensityDrivenFlowTypeTag, PhaseIdx)
{
private:
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
public:
static constexpr int value = FluidSystem::wPhaseIdx;
};
SET_INT_PROP(DensityDrivenFlowTypeTag, PhaseIdx,
GET_PROP_TYPE(TypeTag, FluidSystem)::wPhaseIdx);
SET_INT_PROP(DensityDrivenFlowTypeTag, ReplaceCompEqIdx, 0);
SET_INT_PROP(DensityDrivenFlowTypeTag, ReplaceCompEqIdx, GET_PROP_VALUE(TypeTag, PhaseIdx));
// Set the grid type
SET_TYPE_PROP(DensityDrivenFlowTypeTag, Grid, Dune::YaspGrid<2>);
......@@ -96,13 +91,13 @@ class DensityDrivenFlowProblem : public NavierStokesProblem<TypeTag>
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum { dimWorld = GridView::dimensionworld };
enum {
massBalanceIdx = Indices::massBalanceIdx,
transportEqIdx = 1,
totalMassBalanceIdx = Indices::totalMassBalanceIdx,
transportEqIdx = 1-GET_PROP_VALUE(TypeTag, PhaseIdx),
momentumBalanceIdx = Indices::momentumBalanceIdx,
pressureIdx = Indices::pressureIdx,
velocityXIdx = Indices::velocityXIdx,
velocityYIdx = Indices::velocityYIdx,
transportCompIdx = 1/*FluidSystem::wCompIdx*/
transportCompIdx = 1-GET_PROP_VALUE(TypeTag, PhaseIdx)
};
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
......@@ -156,6 +151,7 @@ public:
{
return NumEqVector(0.0);
}
// \}
/*!
* \name Boundary conditions
......@@ -175,10 +171,10 @@ public:
// set Dirichlet values for the velocity everywhere
values.setDirichlet(momentumBalanceIdx);
values.setOutflow(transportEqIdx);
values.setOutflow(massBalanceIdx);
values.setOutflow(totalMassBalanceIdx);
if(isLowerLeftCell_(globalPos))
values.setDirichletCell(massBalanceIdx);
values.setDirichletCell(totalMassBalanceIdx);
if(globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_)
{
......
test/freeflow/navierstokesnc/msfreeflowtestproblem.hh
View file @ 0934d808
......@@ -184,7 +184,7 @@ class MaxwellStefanNCTestProblem : public NavierStokesProblem<TypeTag>
// copy some indices for convenience
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum {
massBalanceIdx = Indices::massBalanceIdx,
totalMassBalanceIdx = Indices::totalMassBalanceIdx,
compTwoIdx = FluidSystem::N2Idx,
compThreeIdx = FluidSystem::CO2Idx,
momentumBalanceIdx = Indices::momentumBalanceIdx,
......@@ -352,7 +352,7 @@ public:
values.setDirichlet(momentumBalanceIdx);
values.setOutflow(compTwoIdx);
values.setOutflow(compThreeIdx);
values.setOutflow(massBalanceIdx);
values.setOutflow(totalMassBalanceIdx);
return values;
}
......
test/freeflow/rans/pipelauferproblem.hh
View file @ 0934d808
......@@ -82,7 +82,7 @@ class PipeLauferProblem : public ZeroEqProblem<TypeTag>
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum { dimWorld = GridView::dimensionworld };
enum {
massBalanceIdx = Indices::massBalanceIdx,
totalMassBalanceIdx = Indices::totalMassBalanceIdx,
momentumBalanceIdx = Indices::momentumBalanceIdx,
pressureIdx = Indices::pressureIdx,
velocityXIdx = Indices::velocityXIdx,
......@@ -169,11 +169,11 @@ public:
// set a fixed pressure in one cell
if (isOutlet(globalPos))
{
values.setDirichlet(massBalanceIdx);
values.setDirichlet(totalMassBalanceIdx);
values.setOutflow(momentumBalanceIdx);
}
else
values.setOutflow(massBalanceIdx);
values.setOutflow(totalMassBalanceIdx);
return values;
}
......
test/freeflow/ransnc/channeltestproblem.hh
View file @ 0934d808
......@@ -91,7 +91,7 @@ class ChannelNCTestProblem : public RANSProblem<TypeTag>
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum { dimWorld = GridView::dimensionworld };
enum {
massBalanceIdx = Indices::massBalanceIdx,
totalMassBalanceIdx = Indices::totalMassBalanceIdx,
transportEqIdx = 1,
momentumBalanceIdx = Indices::momentumBalanceIdx,
pressureIdx = Indices::pressureIdx,
......@@ -173,7 +173,7 @@ public:
if(isInlet(globalPos))
{
values.setDirichlet(momentumBalanceIdx);
values.setOutflow(massBalanceIdx);
values.setOutflow(totalMassBalanceIdx);
values.setDirichlet(transportEqIdx);
#if NONISOTHERMAL
values.setDirichlet(energyBalanceIdx);
......@@ -182,7 +182,7 @@ public:
else if(isOutlet(globalPos))
{
values.setOutflow(momentumBalanceIdx);
values.setDirichlet(massBalanceIdx);
values.setDirichlet(totalMassBalanceIdx);
values.setOutflow(transportEqIdx);
#if NONISOTHERMAL
values.setOutflow(energyBalanceIdx);
......@@ -193,7 +193,7 @@ public:
{
// set Dirichlet values for the velocity everywhere
values.setDirichlet(momentumBalanceIdx);
values.setOutflow(massBalanceIdx);
values.setOutflow(totalMassBalanceIdx);
values.setOutflow(transportEqIdx);
#if NONISOTHERMAL
values.setOutflow(energyBalanceIdx);
......
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