diff --git a/dumux/discretization/box/maxwellstefanslaw.hh b/dumux/discretization/box/maxwellstefanslaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e824a7da33d98f4fef35c9b3b4710f403457fc13
--- /dev/null
+++ b/dumux/discretization/box/maxwellstefanslaw.hh
@@ -0,0 +1,280 @@
+// -*- 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 . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief This file contains the data which is required to calculate
+ * diffusive mass fluxes due to molecular diffusion with Fick's law.
+ */
+#ifndef DUMUX_DISCRETIZATION_BOX_MAXWELL_STEFAN_LAW_HH
+#define DUMUX_DISCRETIZATION_BOX_MAXWELL_STEFAN_LAW_HH
+
+#include
+
+#include
+#include
+
+#include
+#include
+
+namespace Dumux
+{
+
+namespace Properties
+{
+// forward declaration of properties
+NEW_PROP_TAG(NumPhases);
+NEW_PROP_TAG(EffectiveDiffusivityModel);
+}
+
+/*!
+ * \ingroup CCTpfaMaxwellStefansLaw
+ * \brief Specialization of Maxwell Stefan's Law for the Box method.
+ */
+template
+class MaxwellStefansLawImplementation
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ 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;
+
+ enum { dim = GridView::dimension} ;
+ enum { dimWorld = GridView::dimensionworld} ;
+ enum
+ {
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+ numComponents = GET_PROP_VALUE(TypeTag,NumComponents)
+ };
+ using DimWorldMatrix = Dune::FieldMatrix;
+ using GlobalPosition = Dune::FieldVector;
+ using ComponentFluxVector = Dune::FieldVector;
+ using ReducedComponentVector = Dune::FieldVector;
+ using ReducedComponentMatrix = Dune::FieldMatrix;
+
+public:
+ static ComponentFluxVector flux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const int phaseIdx,
+ const ElementFluxVariablesCache& elemFluxVarsCache)
+ {
+ //this is to calculate the maxwellStefan diffusion in a multicomponent system.
+ //see: Multicomponent Mass Transfer. R. Taylor u. R. Krishna. J. Wiley & Sons, New York 1993
+ ComponentFluxVector componentFlux(0.0);
+ ReducedComponentMatrix reducedDiffusionMatrix(0.0);
+ ReducedComponentVector reducedFlux(0.0);
+ ComponentFluxVector moleFrac(0.0);
+ ReducedComponentVector normalX(0.0);
+
+ // evaluate gradX at integration point and interpolate density
+ const auto& fluxVarsCache = elemFluxVarsCache[scvf];
+ const auto& jacInvT = fluxVarsCache.jacInvT();
+ const auto& shapeJacobian = fluxVarsCache.shapeJacobian();
+ const auto& shapeValues = fluxVarsCache.shapeValues();
+
+ Scalar rho(0.0);
+ std::vector gradN(fvGeometry.numScv());
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+
+ // density interpolation
+ rho += volVars.molarDensity(phaseIdx)*shapeValues[scv.indexInElement()][0];
+ // the ansatz function gradient
+ jacInvT.mv(shapeJacobian[scv.indexInElement()][0], gradN[scv.indexInElement()]);
+
+ //interpolate the mole fraction for the diffusion matrix
+ for (int compIdx = 0; compIdx < numComponents; compIdx++)
+ {
+ moleFrac[compIdx] += volVars.moleFraction(phaseIdx, compIdx)*shapeValues[scv.indexInElement()][0];
+ }
+ }
+
+ reducedDiffusionMatrix = setupMSMatrix_(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, moleFrac);
+
+ for (int compIdx = 0; compIdx < numComponents-1; compIdx++)
+ {
+ GlobalPosition gradX(0.0);
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+
+ // the mole/mass fraction gradient
+ gradX.axpy(volVars.moleFraction(phaseIdx, compIdx), gradN[scv.indexInElement()]);
+ }
+
+ normalX[compIdx] = gradX *scvf.unitOuterNormal();
+ }
+ reducedDiffusionMatrix.solve(reducedFlux,normalX);
+ reducedFlux *= -1.0*rho*scvf.area();
+
+
+ for (int compIdx = 0; compIdx < numComponents-1; compIdx++)
+ {
+ componentFlux[compIdx] = reducedFlux[compIdx];
+ componentFlux[numComponents-1] -= reducedFlux[compIdx];
+ }
+ return componentFlux ;
+ }
+
+private:
+
+ static ReducedComponentMatrix setupMSMatrix_(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const int phaseIdx,
+ const ComponentFluxVector moleFrac)
+ {
+
+ ReducedComponentMatrix reducedDiffusionMatrix(0.0);
+
+ const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+ const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
+ const auto insideScvIdx = scvf.insideScvIdx();
+ const auto outsideScvIdx = scvf.outsideScvIdx();
+ //this is to not devide by 0 if the saturation in 0 and the effectiveDiffusivity becomes zero due to that
+ if(insideVolVars.saturation(phaseIdx) == 0 || outsideVolVars.saturation(phaseIdx) == 0)
+ return reducedDiffusionMatrix;
+
+ for (int compIIdx = 0; compIIdx < numComponents-1; compIIdx++)
+ {
+ // effective diffusion tensors
+ using EffDiffModel = typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel);
+
+ const auto xi = moleFrac[compIIdx];
+
+ //calculate diffusivity for i,numComponents
+
+ auto tinInside = getDiffusionCoefficient(phaseIdx, compIIdx, numComponents-1, problem, element, insideVolVars, fvGeometry.scv(insideScvIdx));
+ tinInside = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), tinInside);
+ auto tinOutside = getDiffusionCoefficient(phaseIdx, compIIdx, numComponents-1, problem, element, outsideVolVars, fvGeometry.scv(outsideScvIdx));
+ tinOutside = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), tinOutside);
+
+ // scale by extrusion factor
+ tinInside *= insideVolVars.extrusionFactor();
+ tinOutside *= outsideVolVars.extrusionFactor();
+
+ // the resulting averaged diffusion tensor
+ const auto tin = problem.spatialParams().harmonicMean(tinInside, tinOutside, scvf.unitOuterNormal());
+
+ //set the entrys of the diffusion matrix of the diagonal
+ reducedDiffusionMatrix[compIIdx][compIIdx] += xi/tin;
+
+ for (int compkIdx = 0; compkIdx < numComponents; compkIdx++)
+ {
+ if (compkIdx == compIIdx)
+ continue;
+
+ const auto xk = moleFrac[compkIdx];
+
+ auto tikInside = getDiffusionCoefficient(phaseIdx, compIIdx, compkIdx, problem, element, insideVolVars, fvGeometry.scv(insideScvIdx));
+ tikInside = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), tikInside);
+ auto tikOutside = getDiffusionCoefficient(phaseIdx, compIIdx, compkIdx, problem, element, outsideVolVars, fvGeometry.scv(outsideScvIdx));
+ tikOutside = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), tikOutside);
+
+ // scale by extrusion factor
+ tikInside *= insideVolVars.extrusionFactor();
+ tikOutside *= outsideVolVars.extrusionFactor();
+
+ // the resulting averaged diffusion tensor
+ const auto tik = problem.spatialParams().harmonicMean(tikInside, tikOutside, scvf.unitOuterNormal());
+
+ reducedDiffusionMatrix[compIIdx][compIIdx] += xk/tik;
+ }
+
+ // now set the rest of the entries (off-diagonal)
+ for (int compJIdx = 0; compJIdx < numComponents-1; compJIdx++)
+ {
+ //we don't want to calculate e.g. water in water diffusion
+ if (compIIdx == compJIdx)
+ continue;
+ //calculate diffusivity for compIIdx, compJIdx
+ auto tijInside = getDiffusionCoefficient(phaseIdx, compIIdx, compJIdx, problem, element, insideVolVars, fvGeometry.scv(insideScvIdx));
+ tijInside = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), tijInside);
+ auto tijOutside = getDiffusionCoefficient(phaseIdx, compIIdx, compJIdx, problem, element, outsideVolVars, fvGeometry.scv(outsideScvIdx));
+ tijOutside = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), outsideVolVars.saturation(phaseIdx), tijOutside);
+
+ // scale by extrusion factor
+ tijInside *= insideVolVars.extrusionFactor();
+ tijOutside *= outsideVolVars.extrusionFactor();
+
+ // the resulting averaged diffusion tensor
+ const auto tij = problem.spatialParams().harmonicMean(tijInside, tijOutside, scvf.unitOuterNormal());
+
+ reducedDiffusionMatrix[compIIdx][compJIdx] +=xi*(1/tin - 1/tij);
+ }
+ }
+
+ return reducedDiffusionMatrix;
+ }
+
+private:
+ template =0 >
+ static Scalar getDiffusionCoefficient(const int phaseIdx,
+ const int compIIdx,
+ const int compJIdx,
+ const Problem& problem,
+ const Element& element,
+ const VolumeVariables& volVars,
+ const SubControlVolume& scv)
+ {
+ return FluidSystem::binaryDiffusionCoefficient(compIIdx,
+ compJIdx,
+ problem,
+ element,
+ scv);
+ }
+
+ template =0 >
+ static Scalar getDiffusionCoefficient(const int phaseIdx,
+ const int compIIdx,
+ const int compJIdx,
+ const Problem& problem,
+ const Element& element,
+ const VolumeVariables& volVars,
+ const SubControlVolume& scv)
+ {
+ auto fluidState = volVars.fluidState();
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updateAll(fluidState);
+ return FluidSystem::binaryDiffusionCoefficient(fluidState,
+ paramCache,
+ phaseIdx,
+ compIIdx,
+ compJIdx);
+ }
+
+};
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/cellcentered/tpfa/maxwellstefanslaw.hh b/dumux/discretization/cellcentered/tpfa/maxwellstefanslaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..82e2f9cee59d7044187182335d3672c5d1a61b87
--- /dev/null
+++ b/dumux/discretization/cellcentered/tpfa/maxwellstefanslaw.hh
@@ -0,0 +1,298 @@
+// -*- 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 . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief This file contains the data which is required to calculate
+ * diffusive mass fluxes due to molecular diffusion with Fick's law.
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_TPFA_MAXWELL_STEFAN_LAW_HH
+#define DUMUX_DISCRETIZATION_CC_TPFA_MAXWELL_STEFAN_LAW_HH
+
+#include
+
+#include
+#include
+
+#include
+#include
+
+namespace Dumux
+{
+
+namespace Properties
+{
+// forward declaration of properties
+NEW_PROP_TAG(NumPhases);
+NEW_PROP_TAG(EffectiveDiffusivityModel);
+}
+
+/*!
+ * \ingroup CCTpfaMaxwellStefansLaw
+ * \brief Specialization of Maxwell Stefan's Law for the CCTpfa method.
+ */
+template
+class MaxwellStefansLawImplementation
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using EffDiffModel = typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+ static const int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+ static const int numComponents = GET_PROP_VALUE(TypeTag,NumComponents);
+
+ using DimWorldMatrix = Dune::FieldMatrix;
+ using GlobalPosition = Dune::FieldVector;
+ using ComponentFluxVector = Dune::FieldVector;
+ using ReducedComponentVector = Dune::FieldVector;
+ using ReducedComponentMatrix = Dune::FieldMatrix;
+
+public:
+ // state the discretization method this implementation belongs to
+ static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::CCTpfa;
+
+ //! state the type for the corresponding cache and its filler
+ //! We don't cache anything for this law
+ using Cache = FluxVariablesCaching::EmptyDiffusionCache;
+ using CacheFiller = FluxVariablesCaching::EmptyCacheFiller;
+
+ static ComponentFluxVector flux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const int phaseIdx,
+ const ElementFluxVariablesCache& elemFluxVarsCache)
+ {
+ //this is to calculate the maxwellStefan diffusion in a multicomponent system.
+ //see: Multicomponent Mass Transfer. R. Taylor u. R. Krishna. J. Wiley & Sons, New York 1993
+ ComponentFluxVector componentFlux(0.0);
+ ReducedComponentVector moleFracInside(0.0);
+ ReducedComponentVector moleFracOutside(0.0);
+ ReducedComponentVector reducedFlux(0.0);
+ ReducedComponentMatrix reducedDiffusionMatrixInside(0.0);
+ ReducedComponentMatrix reducedDiffusionMatrixOutside(0.0);
+
+ // get inside/outside volume variables
+ const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+ const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
+ const auto rhoInside = insideVolVars.molarDensity(phaseIdx);
+ const auto rhoOutside = outsideVolVars.molarDensity(phaseIdx);
+ //calculate the mole fraction vectors
+ for (int compIdx = 0; compIdx < numComponents-1; compIdx++)
+ {
+ //calculate x_inside
+ const auto xInside = insideVolVars.moleFraction(phaseIdx, compIdx);
+ //calculate outside molefraction with the respective transmissibility
+ const auto xOutside = outsideVolVars.moleFraction(phaseIdx, compIdx);
+
+ moleFracInside[compIdx] = xInside;
+ moleFracOutside[compIdx] = xOutside;
+ }
+
+ const auto insideScvIdx = scvf.insideScvIdx();
+ const auto& insideScv = fvGeometry.scv(insideScvIdx);
+ const auto outsideScvIdx = scvf.outsideScvIdx();
+ const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
+
+ //now we have to do the tpfa: J_i = J_j which leads to: tij(xi -xj) = -rho Bi^-1 omegai(x*-xi) with x* = (omegai Bi^-1 + omegaj Bj^-1)^-1 (xi omegai Bi^-1 + xj omegaj Bj^-1) with i inside and j outside
+ reducedDiffusionMatrixInside = setupMSMatrix_(problem, element, fvGeometry, insideVolVars, insideScv, phaseIdx);
+
+ reducedDiffusionMatrixOutside = setupMSMatrix_(problem, element, fvGeometry, outsideVolVars, outsideScv, phaseIdx);
+
+ //we cannot solve that if the matrix is 0 everywhere
+ if(!(insideVolVars.saturation(phaseIdx) == 0 || outsideVolVars.saturation(phaseIdx) == 0))
+ {
+
+
+ const Scalar omegai = calculateOmega_(scvf,
+ insideScv,
+ insideVolVars.extrusionFactor());
+ //if on boundary
+ if (scvf.boundary() || scvf.numOutsideScvs() > 1)
+ {
+ moleFracOutside -= moleFracInside;
+ reducedDiffusionMatrixInside.solve(reducedFlux, moleFracOutside);
+ reducedFlux *= omegai;
+ }
+ else
+ {
+ Scalar omegaj;
+ if (dim == dimWorld)
+ // assume the normal vector from outside is anti parallel so we save flipping a vector
+ omegaj = -1.0*calculateOmega_(scvf,
+ outsideScv,
+ outsideVolVars.extrusionFactor());
+ else
+ omegaj = calculateOmega_(fvGeometry.flipScvf(scvf.index()),
+ outsideScv,
+ outsideVolVars.extrusionFactor());
+
+ reducedDiffusionMatrixInside.invert();
+ reducedDiffusionMatrixOutside.invert();
+ reducedDiffusionMatrixInside *= omegai;
+ reducedDiffusionMatrixOutside *= omegaj;
+
+ //in the helpervector we store the values for x*
+ ReducedComponentVector helperVector(0.0);
+ ReducedComponentVector gradientVectori(0.0);
+ ReducedComponentVector gradientVectorj(0.0);
+
+ reducedDiffusionMatrixInside.mv(moleFracInside, gradientVectori);
+ reducedDiffusionMatrixOutside.mv(moleFracOutside, gradientVectorj);
+
+ auto gradientVectorij = (gradientVectori + gradientVectorj);
+
+ //add the two matrixes to each other
+ reducedDiffusionMatrixOutside += reducedDiffusionMatrixInside;
+
+ reducedDiffusionMatrixOutside.solve(helperVector, gradientVectorij);
+
+ //Bi^-1 omegai rho(x*-xi)
+ helperVector -=moleFracInside;
+ reducedDiffusionMatrixInside.mv(helperVector, reducedFlux);
+ }
+
+ reducedFlux *= -0.5*(rhoInside+rhoOutside)*scvf.area();
+ for (int compIdx = 0; compIdx < numComponents-1; compIdx++)
+ {
+ componentFlux[compIdx] = reducedFlux[compIdx];
+ componentFlux[numComponents-1] -=reducedFlux[compIdx];
+ }
+ }
+ return componentFlux ;
+ }
+
+private:
+ static Scalar calculateOmega_(const SubControlVolumeFace& scvf,
+ const SubControlVolume &scv,
+ const Scalar extrusionFactor)
+ {
+ auto distanceVector = scvf.ipGlobal();
+ distanceVector -= scv.center();
+ distanceVector /= distanceVector.two_norm2();
+
+ Scalar omega = (distanceVector * scvf.unitOuterNormal());
+ omega *= extrusionFactor;
+
+ return omega;
+ }
+
+ static ReducedComponentMatrix setupMSMatrix_(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const VolumeVariables& volVars,
+ const SubControlVolume& scv,
+ const int phaseIdx)
+ {
+ ReducedComponentMatrix reducedDiffusionMatrix(0.0);
+
+ //this is to not devide by 0 if the saturation in 0 and the effectiveDiffusivity becomes zero due to that
+ if(volVars.saturation(phaseIdx) == 0)
+ return reducedDiffusionMatrix;
+
+ for (int compIIdx = 0; compIIdx < numComponents-1; compIIdx++)
+ {
+ const auto xi = volVars.moleFraction(phaseIdx, compIIdx);
+
+ //calculate diffusivity for i,numComponents
+ Scalar tin = getDiffusionCoefficient(phaseIdx, compIIdx, numComponents-1, problem, element, volVars, scv);
+ tin = EffDiffModel::effectiveDiffusivity(volVars.porosity(), volVars.saturation(phaseIdx), tin);
+ //set the entrys of the diffusion matrix of the diagonal
+ reducedDiffusionMatrix[compIIdx][compIIdx] += xi/tin;
+ for (int compkIdx = 0; compkIdx < numComponents; compkIdx++)
+ {
+ if (compkIdx == compIIdx)
+ continue;
+
+ const auto xk = volVars.moleFraction(phaseIdx, compkIdx);
+ Scalar tik = getDiffusionCoefficient(phaseIdx, compIIdx, compkIdx, problem, element, volVars, scv);
+ tik = EffDiffModel::effectiveDiffusivity(volVars.porosity(), volVars.saturation(phaseIdx), tik);
+ reducedDiffusionMatrix[compIIdx][compIIdx] += xk/tik;
+ }
+
+ // now set the rest of the entries (off-diagonal)
+ for (int compJIdx = 0; compJIdx < numComponents-1; compJIdx++)
+ {
+ //we don't want to calculate e.g. water in water diffusion
+ if (compIIdx == compJIdx)
+ continue;
+ //calculate diffusivity for compIIdx, compJIdx
+ Scalar tij = getDiffusionCoefficient(phaseIdx, compIIdx, compJIdx, problem, element, volVars, scv);
+ tij = EffDiffModel::effectiveDiffusivity(volVars.porosity(), volVars.saturation(phaseIdx), tij);
+ reducedDiffusionMatrix[compIIdx][compJIdx] +=xi*(1/tin - 1/tij);
+ }
+ }
+ return reducedDiffusionMatrix;
+ }
+
+ template =0 >
+ static Scalar getDiffusionCoefficient(const int phaseIdx,
+ const int compIIdx,
+ const int compJIdx,
+ const Problem& problem,
+ const Element& element,
+ const VolumeVariables& volVars,
+ const SubControlVolume& scv)
+ {
+ return FluidSystem::binaryDiffusionCoefficient(compIIdx,
+ compJIdx,
+ problem,
+ element,
+ scv);
+ }
+
+ template =0 >
+ static Scalar getDiffusionCoefficient(const int phaseIdx,
+ const int compIIdx,
+ const int compJIdx,
+ const Problem& problem,
+ const Element& element,
+ const VolumeVariables& volVars,
+ const SubControlVolume& scv)
+ {
+ auto fluidState = volVars.fluidState();
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updateAll(fluidState);
+ return FluidSystem::binaryDiffusionCoefficient(fluidState,
+ paramCache,
+ phaseIdx,
+ compIIdx,
+ compJIdx);
+ }
+
+
+
+};
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/maxwellstefanslaw.hh b/dumux/discretization/maxwellstefanslaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..624489d7f99259ed5300d4bbe42f0fcc92f68478
--- /dev/null
+++ b/dumux/discretization/maxwellstefanslaw.hh
@@ -0,0 +1,49 @@
+// -*- 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 . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief This file contains the data which is required to calculate
+ * diffusive mass fluxes due to molecular diffusion with Fick's law.
+ */
+#ifndef DUMUX_DISCRETIZATION_MAXWELL_STEFAN_LAW_HH
+#define DUMUX_DISCRETIZATION_MAXWELL_STEFAN_LAW_HH
+
+#include
+
+namespace Dumux
+{
+// forward declaration
+template
+class MaxwellStefansLawImplementation
+{};
+
+/*!
+ * \ingroup MaxwellStefansLaw
+ * \brief Evaluates the diffusive mass flux according to Maxwell Stafan's law
+ */
+template
+using MaxwellStefansLaw = MaxwellStefansLawImplementation;
+
+} // end namespace
+
+#include
+#include
+#include
+
+#endif
diff --git a/dumux/discretization/staggered/freeflow/maxwellstefanslaw.hh b/dumux/discretization/staggered/freeflow/maxwellstefanslaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ac537b0c1b358472f034cf01ee350173330d4d6e
--- /dev/null
+++ b/dumux/discretization/staggered/freeflow/maxwellstefanslaw.hh
@@ -0,0 +1,265 @@
+// -*- 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 . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief This file contains the data which is required to calculate
+ * diffusive mass fluxes due to molecular diffusion with Fick's law.
+ */
+#ifndef DUMUX_DISCRETIZATION_STAGGERED_MAXWELL_STEFAN_LAW_HH
+#define DUMUX_DISCRETIZATION_STAGGERED_MAXWELL_STEFAN_LAW_HH
+
+#include
+
+#include
+#include
+
+#include
+#include
+
+#include
+namespace Dumux
+{
+
+/*!
+ * \ingroup StaggeredMaxwellStefansLaw
+ * \brief Specialization of Maxwell Stefan's Law for the Staggered method.
+ */
+template
+class MaxwellStefansLawImplementation
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ 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 FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+
+ static const int numComponents = GET_PROP_VALUE(TypeTag,NumComponents);
+ static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+
+ using ReducedComponentVector = Dune::FieldVector;
+ using ReducedComponentMatrix = Dune::FieldMatrix;
+
+ static_assert(GET_PROP_VALUE(TypeTag, NumPhases) == 1, "Only one phase allowed supported!");
+
+ enum {
+ pressureIdx = Indices::pressureIdx,
+ conti0EqIdx = Indices::conti0EqIdx,
+ replaceCompEqIdx = Indices::replaceCompEqIdx,
+ phaseIdx = Indices::phaseIdx
+ };
+
+public:
+ // state the discretization method this implementation belongs to
+ static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::Staggered;
+
+ //! state the type for the corresponding cache and its filler
+ //! We don't cache anything for this law
+ using Cache = FluxVariablesCaching::EmptyDiffusionCache;
+ using CacheFiller = FluxVariablesCaching::EmptyCacheFiller;
+
+ static CellCenterPrimaryVariables diffusiveFluxForCellCenter(const Problem& problem,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf)
+ {
+ //this is to calculate the maxwellStefan diffusion in a multicomponent system.
+ //see: Multicomponent Mass Transfer. R. Taylor u. R. Krishna. J. Wiley & Sons, New York 1993
+ CellCenterPrimaryVariables componentFlux(0.0);
+ ReducedComponentVector moleFracInside(0.0);
+ ReducedComponentVector moleFracOutside(0.0);
+ ReducedComponentVector reducedFlux(0.0);
+ ReducedComponentMatrix reducedDiffusionMatrixInside(0.0);
+ ReducedComponentMatrix reducedDiffusionMatrixOutside(0.0);
+
+ // get inside/outside volume variables
+ const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+ const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
+ const auto rhoInside = insideVolVars.molarDensity();
+ const auto rhoOutside = scvf.boundary() ? insideVolVars.molarDensity() : outsideVolVars.molarDensity();
+ //calculate the mole fraction vectors
+ for (int compIdx = 0; compIdx < numComponents-1; compIdx++)
+ {
+ //calculate x_inside
+ const auto xInside = insideVolVars.moleFraction(phaseIdx, compIdx);
+ //calculate outside molefraction with the respective transmissibility
+ const auto xOutside = outsideVolVars.moleFraction(phaseIdx, compIdx);
+
+ moleFracInside[compIdx] = xInside;
+ moleFracOutside[compIdx] = xOutside;
+
+ // get equation index
+ const auto eqIdx = conti0EqIdx + compIdx;
+
+ if(scvf.boundary())
+ {
+ const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
+ if(bcTypes.isOutflow(eqIdx) && eqIdx != pressureIdx)
+ return componentFlux;
+ }
+ }
+
+ //now we have to do the tpfa: J_i = J_j which leads to: tij(xi -xj) = -rho Bi^-1 omegai(x*-xi) with x* = (omegai Bi^-1 + omegaj Bj^-1)^-1 (xi omegai Bi^-1 + xj omegaj Bj^-1) with i inside and j outside
+ reducedDiffusionMatrixInside = setupMSMatrix_(problem, fvGeometry, insideVolVars, scvf);
+
+ reducedDiffusionMatrixOutside = setupMSMatrix_(problem, fvGeometry, outsideVolVars, scvf);
+
+ const auto insideScvIdx = scvf.insideScvIdx();
+ const auto& insideScv = fvGeometry.scv(insideScvIdx);
+ const auto outsideScvIdx = scvf.outsideScvIdx();
+ const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
+
+ const Scalar omegai = calculateOmega_(scvf,
+ insideScv,
+ 1);
+
+ //if on boundary
+ if (scvf.boundary())
+ {
+ moleFracOutside -= moleFracInside;
+ reducedDiffusionMatrixInside.solve(reducedFlux, moleFracOutside);
+ reducedFlux *= omegai;
+ }
+ else
+ {
+ Scalar omegaj;
+ omegaj = -1.0*calculateOmega_(scvf,
+ outsideScv,
+ 1);
+
+ reducedDiffusionMatrixInside.invert();
+ reducedDiffusionMatrixOutside.invert();
+ reducedDiffusionMatrixInside *= omegai;
+ reducedDiffusionMatrixOutside *= omegaj;
+
+ //in the helpervector we store the values for x*
+ ReducedComponentVector helperVector(0.0);
+ ReducedComponentVector gradientVectori(0.0);
+ ReducedComponentVector gradientVectorj(0.0);
+
+ reducedDiffusionMatrixInside.mv(moleFracInside, gradientVectori);
+ reducedDiffusionMatrixOutside.mv(moleFracOutside, gradientVectorj);
+
+ auto gradientVectorij = (gradientVectori + gradientVectorj);
+
+ //add the two matrixes to each other
+ reducedDiffusionMatrixOutside += reducedDiffusionMatrixInside;
+
+ reducedDiffusionMatrixOutside.solve(helperVector, gradientVectorij);
+
+ //Bi^-1 omegai rho(x*-xi)
+ helperVector -=moleFracInside;
+ reducedDiffusionMatrixInside.mv(helperVector, reducedFlux);
+ }
+
+ reducedFlux *= -0.5*(rhoInside+rhoOutside)*scvf.area();
+
+ for (int compIdx = 0; compIdx < numComponents-1; compIdx++)
+ {
+ componentFlux[compIdx] = reducedFlux[compIdx];
+ componentFlux[numComponents-1] -=reducedFlux[compIdx];
+ }
+
+ if(useMoles && replaceCompEqIdx <= numComponents)
+ componentFlux[replaceCompEqIdx] = 0.0;
+
+ return componentFlux ;
+ }
+
+private:
+ static Scalar calculateOmega_(const SubControlVolumeFace& scvf,
+ const SubControlVolume &scv,
+ const Scalar extrusionFactor)
+ {
+ auto distanceVector = scvf.ipGlobal();
+ distanceVector -= scv.center();
+ distanceVector /= distanceVector.two_norm2();
+
+ Scalar omega = (distanceVector * scvf.unitOuterNormal());
+ omega *= extrusionFactor;
+
+ return omega;
+ }
+
+ static ReducedComponentMatrix setupMSMatrix_(const Problem& problem,
+ const FVElementGeometry& fvGeometry,
+ const VolumeVariables& volVars,
+ const SubControlVolumeFace& scvf)
+ {
+ ReducedComponentMatrix reducedDiffusionMatrix(0.0);
+
+ for (int compIIdx = 0; compIIdx < numComponents-1; compIIdx++)
+ {
+ const auto xi = volVars.moleFraction(phaseIdx, compIIdx);
+
+ //calculate diffusivity for i,numComponents
+ auto fluidState = volVars.fluidState();
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updateAll(fluidState);
+ auto tin = FluidSystem::binaryDiffusionCoefficient(fluidState,
+ paramCache,
+ phaseIdx,
+ compIIdx,
+ numComponents-1);
+ //set the entrys of the diffusion matrix of the diagonal
+ reducedDiffusionMatrix[compIIdx][compIIdx] += xi/tin;
+
+ for (int compkIdx = 0; compkIdx < numComponents; compkIdx++)
+ {
+ if (compkIdx == compIIdx)
+ continue;
+
+ const auto xk = volVars.moleFraction(phaseIdx, compkIdx);
+ Scalar tik = FluidSystem::binaryDiffusionCoefficient(fluidState,
+ paramCache,
+ phaseIdx,
+ compIIdx,
+ compkIdx);
+ reducedDiffusionMatrix[compIIdx][compIIdx] += xk/tik;
+ }
+
+ // now set the rest of the entries (off-diagonal)
+ for (int compJIdx = 0; compJIdx < numComponents-1; compJIdx++)
+ {
+ //we don't want to calculate e.g. water in water diffusion
+ if (compIIdx == compJIdx)
+ continue;
+ //calculate diffusivity for compIIdx, compJIdx
+ Scalar tij = FluidSystem::binaryDiffusionCoefficient(fluidState,
+ paramCache,
+ phaseIdx,
+ compIIdx,
+ compJIdx);
+ reducedDiffusionMatrix[compIIdx][compJIdx] +=xi*(1/tin - 1/tij);
+ }
+ }
+ return reducedDiffusionMatrix;
+ }
+};
+} // end namespace
+
+#endif
diff --git a/test/freeflow/navierstokesnc/CMakeLists.txt b/test/freeflow/navierstokesnc/CMakeLists.txt
index 1c3742efe47cf20e71688e7c40b0361b6e794769..07cf0742c2936e0b3145250f53712d070dddde20 100644
--- a/test/freeflow/navierstokesnc/CMakeLists.txt
+++ b/test/freeflow/navierstokesnc/CMakeLists.txt
@@ -54,3 +54,12 @@ dune_add_test(NAME test_stokes2cni_diffusion
--files ${CMAKE_SOURCE_DIR}/test/references/stokes2cni-diffusion.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_stokes2cni_diffusion-00014.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}//test_stokes2cni test_stokes2cni_diffusion.input")
+
+dune_add_test(NAME test_msfreeflow
+ SOURCES test_msfreeflow.cc
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/maxwellstefanfreeflow-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_msfreeflow-00022.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_msfreeflow test_msfreeflow.input")
diff --git a/test/freeflow/navierstokesnc/msfreeflowtestproblem.hh b/test/freeflow/navierstokesnc/msfreeflowtestproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a3696e846aec6941b112d1ab0532529b8d644caf
--- /dev/null
+++ b/test/freeflow/navierstokesnc/msfreeflowtestproblem.hh
@@ -0,0 +1,463 @@
+// -*- 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 . *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Channel flow test for the multi-component staggered grid (Navier-)Stokes model
+ */
+#ifndef DUMUX_CHANNEL_MAXWELL_STEFAN_TEST_PROBLEM_HH
+#define DUMUX_CHANNEL_MAXWELL_STEFAN_TEST_PROBLEM_HH
+
+#include
+#include
+
+#include
+#include
+
+#include
+#include
+
+#include
+
+namespace Dumux
+{
+template
+class MaxwellStefanNCTestProblem;
+
+namespace Properties {
+
+#if !NONISOTHERMAL
+NEW_TYPE_TAG(MaxwellStefanNCTestProblem, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNC));
+#else
+NEW_TYPE_TAG(MaxwellStefanNCTestProblem, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNCNI));
+#endif
+
+NEW_PROP_TAG(FluidSystem);
+
+SET_INT_PROP(MaxwellStefanNCTestProblem, ReplaceCompEqIdx, 0);
+
+// Set the grid type
+SET_TYPE_PROP(MaxwellStefanNCTestProblem, Grid, Dune::YaspGrid<2>);
+
+// Set the problem property
+SET_TYPE_PROP(MaxwellStefanNCTestProblem, Problem, Dumux::MaxwellStefanNCTestProblem );
+
+SET_BOOL_PROP(MaxwellStefanNCTestProblem, EnableFVGridGeometryCache, true);
+
+SET_BOOL_PROP(MaxwellStefanNCTestProblem, EnableGridFluxVariablesCache, true);
+SET_BOOL_PROP(MaxwellStefanNCTestProblem, EnableGridVolumeVariablesCache, true);
+
+SET_BOOL_PROP(MaxwellStefanNCTestProblem, UseMoles, true);
+
+// #if ENABLE_NAVIERSTOKES
+SET_BOOL_PROP(MaxwellStefanNCTestProblem, EnableInertiaTerms, true);
+
+//! Here we set FicksLaw or MaxwellStefansLaw
+SET_TYPE_PROP(MaxwellStefanNCTestProblem, MolecularDiffusionType, MaxwellStefansLaw);
+
+
+//! A simple fluid system with one MaxwellStefan component
+template
+class MaxwellStefanFluidSystem: public FluidSystems::BaseFluidSystem>
+
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ typedef MaxwellStefanFluidSystem ThisType;
+ typedef FluidSystems::BaseFluidSystem Base;
+
+public:
+ //! The number of components
+ static constexpr int numPhases = 1;
+ static constexpr int numComponents = 3;
+
+ static constexpr int H2Idx = 0;//first major component
+ static constexpr int N2Idx = 1;//second major component
+ static constexpr int CO2Idx = 2;//secondary component
+
+ //! Human readable component name (index compIdx) (for vtk output)
+ static std::string componentName(int compIdx)
+ { return "MaxwellStefan_" + std::to_string(compIdx); }
+
+ //! Human readable phase name (index phaseIdx) (for velocity vtk output)
+ static std::string phaseName(int phaseIdx = 0)
+ { return "Gas"; }
+
+ //! Molar mass in kg/mol of the component with index compIdx
+ static Scalar molarMass(unsigned int compIdx)
+ { return 0.02896; }
+
+
+ using Base::binaryDiffusionCoefficient;
+ /*!
+ * \brief Given a phase's composition, temperature and pressure,
+ * return the binary diffusion coefficient \f$\mathrm{[m^2/s]}\f$ for components
+ * \f$i\f$ and \f$j\f$ in this phase.
+ *
+ * \param fluidState An arbitrary fluid state
+ * \param phaseIdx The index of the fluid phase to consider
+ * \param compIIdx The index of the first component to consider
+ * \param compJIdx The index of the second component to consider
+ */
+ template
+ static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
+ int phaseIdx,
+ unsigned int compIIdx,
+ unsigned int compJIdx)
+ {
+ if (compIIdx > compJIdx)
+ {
+ using std::swap;
+ swap(compIIdx, compJIdx);
+ }
+
+ if (compIIdx == H2Idx && compJIdx == N2Idx)
+ return 83.3e-6;
+ if (compIIdx == H2Idx && compJIdx == CO2Idx)
+ return 68.0e-6;
+ if (compIIdx == N2Idx && compJIdx == CO2Idx)
+ return 16.8e-6;
+ DUNE_THROW(Dune::InvalidStateException,
+ "Binary diffusion coefficient of components "
+ << compIIdx << " and " << compJIdx << " is undefined!\n");
+ }
+ using Base::density;
+ /*!
+ * \brief Given a phase's composition, temperature, pressure, and
+ * the partial pressures of all components, return its
+ * density \f$\mathrm{[kg/m^3]}\f$.
+ * \param phaseIdx index of the phase
+ * \param fluidState the fluid state
+ *
+ */
+ template
+ static Scalar density(const FluidState &fluidState,
+ const int phaseIdx)
+ {
+ return 1;
+ }
+
+ using Base::viscosity;
+ /*!
+ * \brief Calculate the dynamic viscosity of a fluid phase \f$\mathrm{[Pa*s]}\f$
+ * \param fluidState An arbitrary fluid state
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ template
+ static Scalar viscosity(const FluidState &fluidState,
+ int phaseIdx)
+ {
+ return 1e-6;
+ }
+};
+
+SET_TYPE_PROP(MaxwellStefanNCTestProblem, FluidSystem, MaxwellStefanFluidSystem);
+
+} //end namespace Property
+/*!
+ * \brief Test problem for the maxwell stefan model
+ \todo doc me!
+ */
+template
+class MaxwellStefanNCTestProblem : public NavierStokesProblem
+{
+ using ParentType = NavierStokesProblem;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ // copy some indices for convenience
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ enum {
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+ enum {
+ massBalanceIdx = Indices::massBalanceIdx,
+ compTwoIdx = FluidSystem::N2Idx,
+ compThreeIdx = FluidSystem::CO2Idx,
+ momentumBalanceIdx = Indices::momentumBalanceIdx,
+ momentumXBalanceIdx = Indices::momentumXBalanceIdx,
+ momentumYBalanceIdx = Indices::momentumYBalanceIdx,
+ pressureIdx = Indices::pressureIdx,
+ velocityXIdx = Indices::velocityXIdx,
+ velocityYIdx = Indices::velocityYIdx,
+ };
+
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+
+ using GlobalPosition = Dune::FieldVector;
+
+ using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+ using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+
+public:
+ MaxwellStefanNCTestProblem(std::shared_ptr fvGridGeometry)
+ : ParentType(fvGridGeometry), eps_(1e-6)
+ {
+ name_ = getParam("Problem.Name");
+ plotOutput_ = false;
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ std::string name() const
+ {
+ return name_;
+ }
+
+ bool shouldWriteRestartFile() const
+ {
+ return false;
+ }
+
+ //! Called after every time step
+ //! Output the diffusion rates from left to right
+ void postTimeStep(const SolutionVector& curSol,
+ const GridVariables& gridVariables,
+ Scalar time)
+ {
+
+ if (plotOutput_)
+ {
+ Scalar x_co2_left = 0.0;
+ Scalar x_n2_left = 0.0;
+ Scalar x_co2_right = 0.0;
+ Scalar x_n2_right = 0.0;
+ Scalar x_h2_left = 0.0;
+ Scalar x_h2_right = 0.0;
+ Scalar i = 0.0;
+ Scalar j = 0.0;
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(gridVariables.curGridVolVars());
+ elemVolVars.bind(element, fvGeometry, curSol);
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto globalPos = scv.dofPosition();
+
+ if (globalPos[0] < 0.5)
+ {
+ x_co2_left += elemVolVars[scv].moleFraction(0,2);
+
+ x_n2_left += elemVolVars[scv].moleFraction(0,1);
+ x_h2_left += elemVolVars[scv].moleFraction(0,0);
+ i +=1;
+ }
+ else
+ {
+ x_co2_right += elemVolVars[scv].moleFraction(0,2);
+ x_n2_right += elemVolVars[scv].moleFraction(0,1);
+ x_h2_right += elemVolVars[scv].moleFraction(0,0);
+ j +=1;
+ }
+
+ }
+ }
+ x_co2_left /= i;
+ x_n2_left /= i;
+ x_h2_left /= i;
+ x_co2_right /= j;
+ x_n2_right /= j;
+ x_h2_right /= j;
+
+ //do a gnuplot
+ x_.push_back(time); // in seconds
+ y_.push_back(x_n2_left);
+ y2_.push_back(x_n2_right);
+ y3_.push_back(x_co2_left);
+ y4_.push_back(x_co2_right);
+ y5_.push_back(x_h2_left);
+ y6_.push_back(x_h2_right);
+
+ gnuplot_.resetPlot();
+ gnuplot_.setXRange(0, std::max(time, 72000.0));
+ gnuplot_.setYRange(0.4, 0.6);
+ gnuplot_.setXlabel("time [s]");
+ gnuplot_.setYlabel("mole fraction mol/mol");
+ gnuplot_.addDataSetToPlot(x_, y_, "N2 left");
+ gnuplot_.addDataSetToPlot(x_, y2_, "N2 right");
+ gnuplot_.plot("mole_fraction_N2");
+
+ gnuplot2_.resetPlot();
+ gnuplot2_.setXRange(0, std::max(time, 72000.0));
+ gnuplot2_.setYRange(0.0, 0.6);
+ gnuplot2_.setXlabel("time [s]");
+ gnuplot2_.setYlabel("mole fraction mol/mol");
+ gnuplot2_.addDataSetToPlot(x_, y3_, "C02 left");
+ gnuplot2_.addDataSetToPlot(x_, y4_, "C02 right");
+ gnuplot2_.plot("mole_fraction_C02");
+
+ gnuplot3_.resetPlot();
+ gnuplot3_.setXRange(0, std::max(time, 72000.0));
+ gnuplot3_.setYRange(0.0, 0.6);
+ gnuplot3_.setXlabel("time [s]");
+ gnuplot3_.setYlabel("mole fraction mol/mol");
+ gnuplot3_.addDataSetToPlot(x_, y5_, "H2 left");
+ gnuplot3_.addDataSetToPlot(x_, y6_, "H2 right");
+ gnuplot3_.plot("mole_fraction_H2");
+ }
+
+ }
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature() const
+ { return 273.15 + 10; } // 10C
+
+ /*!
+ * \brief Return the sources within the domain.
+ *
+ * \param globalPos The global position
+ */
+ SourceValues sourceAtPos(const GlobalPosition &globalPos) const
+ {
+ return SourceValues(0.0);
+ }
+
+ // \}
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary control volume.
+ *
+ * \param globalPos The position of the center of the finite volume
+ */
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+ {
+ BoundaryTypes values;
+ values.setAllNeumann();
+ // set Dirichlet values for the velocity everywhere
+ values.setDirichlet(momentumBalanceIdx);
+ values.setDirichlet(massBalanceIdx);
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * control volume.
+ *
+ * \param globalPos The center of the finite volume which ought to be set.
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values = initialAtPos(globalPos);
+ return values;
+
+ }
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables initialValues(0.0);
+ if (globalPos[0] < 0.5)
+ {
+ initialValues[compTwoIdx] = 0.50086;
+ initialValues[compThreeIdx] = 0.49914;
+ }
+ else
+ {
+ initialValues[compTwoIdx] = 0.49879;
+ initialValues[compThreeIdx] = 0.0;
+ }
+
+ initialValues[pressureIdx] = 1.1e+5;
+ initialValues[velocityYIdx] = 0.0;
+
+ return initialValues;
+ }
+
+private:
+
+ bool isInlet(const GlobalPosition& globalPos) const
+ {
+ return globalPos[0] < eps_;
+ }
+
+ bool isOutlet(const GlobalPosition& globalPos) const
+ {
+ return globalPos[0] > this->bBoxMax()[0] - eps_;
+ }
+
+ bool isWall(const GlobalPosition& globalPos) const
+ {
+ return globalPos[0] > eps_ || globalPos[0] < this->bBoxMax()[0] - eps_;
+ }
+
+ bool plotOutput_;
+
+ const Scalar eps_{1e-6};
+ std::string name_;
+
+ Dumux::GnuplotInterface gnuplot_;
+ Dumux::GnuplotInterface gnuplot2_;
+ Dumux::GnuplotInterface gnuplot3_;
+
+ std::vector x_;
+ std::vector y_;
+ std::vector y2_;
+ std::vector y3_;
+ std::vector y4_;
+ std::vector y5_;
+ std::vector y6_;
+};
+} //end namespace
+
+#endif
diff --git a/test/freeflow/navierstokesnc/test_msfreeflow.cc b/test/freeflow/navierstokesnc/test_msfreeflow.cc
new file mode 100644
index 0000000000000000000000000000000000000000..af2672b58acf11b6e1196a7884da99047219ce2a
--- /dev/null
+++ b/test/freeflow/navierstokesnc/test_msfreeflow.cc
@@ -0,0 +1,261 @@
+// -*- 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 . *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Test for the staggered grid multi-component (Navier-)Stokes model
+ */
+ #include
+
+ #include
+ #include
+
+ #include
+ #include
+ #include
+ #include
+ #include
+
+#include "msfreeflowtestproblem.hh"
+
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+#include
+
+#include
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ * reading in parameters.
+ *
+ * \param progName The name of the program, that was tried to be started.
+ * \param errorMsg The error message that was issued by the start function.
+ * Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+ if (errorMsg.size() > 0) {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\n\nThe list of mandatory arguments for this program is:\n"
+ "\t-TimeManager.TEnd End of the simulation [s] \n"
+ "\t-TimeManager.DtInitial Initial timestep size [s] \n"
+ "\t-Grid.File Name of the file containing the grid \n"
+ "\t definition in DGF format\n"
+ "\t-SpatialParams.LensLowerLeftX x-coordinate of the lower left corner of the lens [m] \n"
+ "\t-SpatialParams.LensLowerLeftY y-coordinate of the lower left corner of the lens [m] \n"
+ "\t-SpatialParams.LensUpperRightX x-coordinate of the upper right corner of the lens [m] \n"
+ "\t-SpatialParams.LensUpperRightY y-coordinate of the upper right corner of the lens [m] \n"
+ "\t-SpatialParams.Permeability Permeability of the domain [m^2] \n"
+ "\t-SpatialParams.PermeabilityLens Permeability of the lens [m^2] \n";
+
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(MaxwellStefanNCTestProblem);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv, usage);
+
+ // try to create a grid (from the given grid file or the input file)
+ using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+ GridCreator::makeGrid();
+ GridCreator::loadBalance();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = GridCreator::grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared(fvGridGeometry);
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam("TimeLoop.TEnd");
+ const auto maxDivisions = getParam("TimeLoop.MaxTimeStepDivisions");
+ const auto maxDt = getParam("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam("TimeLoop.Restart");
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+ const auto numDofsCellCenter = leafGridView.size(0);
+ const auto numDofsFace = leafGridView.size(1);
+ SolutionVector x;
+ x[cellCenterIdx].resize(numDofsCellCenter);
+ x[faceIdx].resize(numDofsFace);
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // intialize the vtk output module
+ using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+ StaggeredVtkOutputModule vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+ VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+ vtkWriter.write(0.0);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = StaggeredFVAssembler;
+ auto assembler = std::make_shared(problem, fvGridGeometry, gridVariables, timeLoop);
+
+ // the linear solver
+ using LinearSolver = Dumux::UMFPackBackend;
+ auto linearSolver = std::make_shared();
+
+ // the non-linear solver
+ using NewtonController = StaggeredNewtonController;
+ using NewtonMethod = Dumux::NewtonMethod;
+ auto newtonController = std::make_shared(leafGridView.comm(), timeLoop);
+ NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(xOld);
+
+ // try solving the non-linear system
+ for (int i = 0; i < maxDivisions; ++i)
+ {
+ // linearize & solve
+ auto converged = nonLinearSolver.solve(x);
+
+ if (converged)
+ break;
+
+ if (!converged && i == maxDivisions-1)
+ DUNE_THROW(Dune::MathError,
+ "Newton solver didn't converge after "
+ << maxDivisions
+ << " time-step divisions. dt="
+ << timeLoop->timeStepSize()
+ << ".\nThe solutions of the current and the previous time steps "
+ << "have been saved to restart files.");
+ }
+
+ // make the new solution the old solution
+ xOld = x;
+ gridVariables->advanceTimeStep();
+ problem->postTimeStep(x, *gridVariables,timeLoop->time()+timeLoop->timeStepSize());
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // write vtk output
+ vtkWriter.write(timeLoop->time());
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by newton controller
+ timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/test/freeflow/navierstokesnc/test_msfreeflow.input b/test/freeflow/navierstokesnc/test_msfreeflow.input
new file mode 100644
index 0000000000000000000000000000000000000000..b82489c0658c1cf67bd196a13be95fca1435d29b
--- /dev/null
+++ b/test/freeflow/navierstokesnc/test_msfreeflow.input
@@ -0,0 +1,19 @@
+[TimeLoop]
+DtInitial = 1 # [s]
+TEnd = 70000 # [s]
+MaxTimeStepSize = 1e6
+
+[Grid]
+UpperRight = 1 1
+Cells = 30 30
+
+[Problem]
+Name = test_msfreeflow # name passed to the output routines
+EnableGravity = false
+
+[ Newton ]
+MaxSteps = 10
+MaxRelativeShift = 1e-8
+
+[Vtk]
+WriteFaceData = false
diff --git a/test/porousmediumflow/2pnc/implicit/2pncdiffusion.hh b/test/porousmediumflow/2pnc/implicit/2pncdiffusion.hh
new file mode 100644
index 0000000000000000000000000000000000000000..6cffe8fc595f82d2aee6093fe3e262c59789b60c
--- /dev/null
+++ b/test/porousmediumflow/2pnc/implicit/2pncdiffusion.hh
@@ -0,0 +1,264 @@
+// -*- 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 . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup TwoPNCTest
+ * \brief Problem where air is injected under a low permeable layer in a depth of 2700m.
+ */
+#ifndef DUMUX_TWOPNC_DIFFUSION_PROBLEM_HH
+#define DUMUX_TWOPNC_DIFFUSION_PROBLEM_HH
+
+#include
+#include
+#include
+#include
+#include
+
+#include "2pncdiffusionspatialparams.hh"
+#include
+
+namespace Dumux
+{
+
+template
+class TwoPNCDiffusionProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(TwoPNCDiffusionProblem, INHERITS_FROM(TwoPNC, TwoPNCDiffusionSpatialParams));
+NEW_TYPE_TAG(TwoPNCDiffusionCCProblem, INHERITS_FROM(CCTpfaModel, TwoPNCDiffusionProblem));
+
+// Set the grid type
+SET_TYPE_PROP(TwoPNCDiffusionProblem, Grid, Dune::YaspGrid<2>);
+
+// Set the problem property
+SET_TYPE_PROP(TwoPNCDiffusionProblem, Problem, TwoPNCDiffusionProblem);
+
+// // Set fluid configuration
+SET_TYPE_PROP(TwoPNCDiffusionProblem,
+ FluidSystem,
+ FluidSystems::H2ON2);
+
+// Define whether mole(true) or mass (false) fractions are used
+SET_BOOL_PROP(TwoPNCDiffusionProblem, UseMoles, true);
+
+//! Here we set FicksLaw or TwoPNCDiffusionsLaw
+SET_TYPE_PROP(TwoPNCDiffusionProblem, MolecularDiffusionType, DIFFUSIONTYPE);
+
+//! Set the default formulation to pw-Sn: This can be over written in the problem.
+SET_INT_PROP(TwoPNCDiffusionProblem, Formulation, TwoPNCFormulation::pwsn);
+}
+
+
+/*!
+ * \ingroup TwoPNCTest
+ * \brief DOC_ME
+ */
+template
+class TwoPNCDiffusionProblem : public PorousMediumFlowProblem
+{
+ using ParentType = PorousMediumFlowProblem;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ enum {
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+
+ // copy some indices for convenience
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ enum {
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+
+ wCompIdx = FluidSystem::wCompIdx,
+ nCompIdx = FluidSystem::nCompIdx,
+
+ contiH2OEqIdx = Indices::contiWEqIdx,
+ contiN2EqIdx = Indices::contiNEqIdx
+ };
+
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using GlobalPosition = Dune::FieldVector;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+
+ //! property that defines whether mole or mass fractions are used
+ static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+
+public:
+ /*!
+ * \brief The constructor
+ *
+ * \param timeManager The time manager
+ * \param gridView The grid view
+ */
+ TwoPNCDiffusionProblem(std::shared_ptr fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ {
+ // initialize the tables of the fluid system
+ FluidSystem::init();
+
+ name_ = getParam("Problem.Name");
+
+ //stateing in the console whether mole or mass fractions are used
+ if(useMoles)
+ {
+ std::cout<<"problem uses mole-fractions"<fvGridGeometry().bBoxMin()[0] + eps_)
+ priVars[Indices::switchIdx] = 1e-3;
+
+ return priVars;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann
+ * boundary segment.
+ *
+ * For this method, the \a priVars parameter stores the mass flux
+ * in normal direction of each component. Negative values mean
+ * influx.
+ *
+ * The units must be according to either using mole or mass fractions. (mole/(m^2*s) or kg/(m^2*s))
+ */
+ NeumannFluxes neumannAtPos(const GlobalPosition& globalPos) const
+ {
+ NeumannFluxes values(0.0);
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluates the initial values for a control volume
+ *
+ * \param values Stores the initial values for the conservation equations in
+ * \f$ [ \textnormal{unit of primary variables} ] \f$
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ return initial_(globalPos);
+ }
+
+ // \}
+
+private:
+ /*!
+ * \brief Evaluates the initial values for a control volume
+ *
+ * The internal method for the initial condition
+ *
+ * \param values Stores the initial values for the conservation equations in
+ * \f$ [ \textnormal{unit of primary variables} ] \f$
+ * \param globalPos The global position
+ */
+ PrimaryVariables initial_(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables priVars(0.0);
+ priVars.setState(Indices::wPhaseOnly);
+
+ //mole-fraction formulation
+ priVars[Indices::switchIdx] = 1e-5;
+ priVars[Indices::pressureIdx] = 1e5;
+ return priVars;
+ }
+
+ Scalar temperature_;
+ static constexpr Scalar eps_ = 1e-6;
+
+ std::string name_ ;
+
+
+};
+
+} //end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/2pnc/implicit/2pncdiffusionspatialparams.hh b/test/porousmediumflow/2pnc/implicit/2pncdiffusionspatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4ea817a75577408987553db2c218167c57de666a
--- /dev/null
+++ b/test/porousmediumflow/2pnc/implicit/2pncdiffusionspatialparams.hh
@@ -0,0 +1,149 @@
+// -*- 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 . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup TwoPNCTest
+ * \brief Definition of the spatial parameters for the fuel cell
+ * problem which uses the isothermal/non-insothermal 2pnc box model
+ */
+
+#ifndef DUMUX_TWOPNC_DIFFUSION_SPATIAL_PARAMS_HH
+#define DUMUX_TWOPNC_DIFFUSION_SPATIAL_PARAMS_HH
+
+#include
+#include
+#include
+#include
+#include
+
+namespace Dumux
+{
+
+//forward declaration
+template
+class TwoPNCDiffusionSpatialParams;
+
+namespace Properties
+{
+// The spatial parameters TypeTag
+NEW_TYPE_TAG(TwoPNCDiffusionSpatialParams);
+
+// Set the spatial parameters
+SET_TYPE_PROP(TwoPNCDiffusionSpatialParams, SpatialParams, TwoPNCDiffusionSpatialParams);
+
+// Set the material Law
+SET_PROP(TwoPNCDiffusionSpatialParams, MaterialLaw)
+{
+private:
+ // define the material law which is parameterized by effective
+ // saturations
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using EffectiveLaw = RegularizedBrooksCorey;
+public:
+ // define the material law parameterized by absolute saturations
+ using type = EffToAbsLaw;
+};
+
+} // end namespace Properties
+
+/*!
+ * \ingroup TwoPNCTest
+ * \brief Definition of the spatial parameters for the TwoPNCDiffusion
+ * problem which uses the isothermal 2p2c box model
+ */
+template
+class TwoPNCDiffusionSpatialParams : public FVSpatialParams
+{
+ using ParentType = FVSpatialParams;
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using CoordScalar = typename GridView::ctype;
+
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+
+ using GlobalPosition = Dune::FieldVector;
+ using DimWorldMatrix = Dune::FieldMatrix;
+
+public:
+ using PermeabilityType = DimWorldMatrix;
+
+ using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+ using MaterialLawParams = typename MaterialLaw::Params;
+
+ /*!
+ * \brief The constructor
+ *
+ * \param gridView The grid view
+ */
+ TwoPNCDiffusionSpatialParams(const Problem& problem)
+ : ParentType(problem), K_(0)
+ {
+ // intrinsic permeabilities
+ K_[0][0] = 5e-11;
+ K_[1][1] = 5e-11;
+
+ // porosities
+ porosity_ = 0.2;
+
+ // residual saturations
+ materialParams_.setSwr(0.02); //here water, see philtophoblaw
+ materialParams_.setSnr(0.0);
+
+ //parameters for the vanGenuchten law
+ materialParams_.setPe(1e2); // alpha = 1/pcb
+ materialParams_.setLambda(1.3);
+ }
+
+ /*!
+ * \brief Returns the hydraulic conductivity \f$[m^2]\f$
+ *
+ * \param globalPos The global position
+ */
+ DimWorldMatrix permeabilityAtPos(const GlobalPosition& globalPos) const
+ { return K_; }
+
+ /*!
+ * \brief Define the porosity \f$[-]\f$ of the spatial parameters
+ *
+ * \param globalPos The global position
+ */
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
+ { return 0.2; }
+
+ /*!
+ * \brief return the parameter object for the Brooks-Corey material law which depends on the position
+ *
+ * \param globalPos The global position
+ */
+ const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
+ { return materialParams_; }
+
+private:
+ DimWorldMatrix K_;
+ Scalar porosity_;
+ static constexpr Scalar eps_ = 1e-6;
+ MaterialLawParams materialParams_;
+};
+
+}//end namespace
+
+#endif
diff --git a/test/porousmediumflow/2pnc/implicit/CMakeLists.txt b/test/porousmediumflow/2pnc/implicit/CMakeLists.txt
index e533a50b90b0b19dc49528a4ca3ff47bbff3f344..6b5d5ba73229a4cf3451254e458ed665e50ed17d 100644
--- a/test/porousmediumflow/2pnc/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/2pnc/implicit/CMakeLists.txt
@@ -1,4 +1,4 @@
-dune_symlink_to_source_files(FILES test_2pnc_fv.input)
+dune_symlink_to_source_files(FILES test_2pnc_fv.input test_2pnc_diffusion.input)
# isothermal tests
dune_add_test(NAME test_2pnc_box
@@ -19,10 +19,32 @@ dune_add_test(NAME test_2pnc_tpfa
${CMAKE_CURRENT_BINARY_DIR}/fuelcell_tpfa-00015.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_2pnc_tpfa -ParameterFile test_2pnc_fv.input -Problem.Name fuelcell_tpfa")
+dune_add_test(NAME test_cc2pnc_maxwellstefan
+ SOURCES test_cc2pnc_diffusion.cc
+ COMPILE_DEFINITIONS TYPETAG=TwoPNCDiffusionProblem DIFFUSIONTYPE=MaxwellStefansLaw
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/2pncdiffusioncc-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_cc2pnc_maxwellstefan-00026.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_cc2pnc_maxwellstefan test_2pnc_diffusion.input -Problem.Name test_cc2pnc_maxwellstefan")
+
+dune_add_test(NAME test_cc2pnc_fickslaw
+ SOURCES test_cc2pnc_diffusion.cc
+ COMPILE_DEFINITIONS TYPETAG=TwoPNCDiffusionProblem DIFFUSIONTYPE=FicksLaw
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/2pncdiffusioncc-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_cc2pnc_fickslaw-00026.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_cc2pnc_fickslaw test_2pnc_diffusion.input -Problem.Name test_cc2pnc_fickslaw")
+
#install sources
install(FILES
fuelcellproblem.hh
fuelcellspatialparams.hh
test_2pnc_fv.cc
test_2pnc_fv.input
+test_cc2pnc_diffusion.cc
+maxwellstefandiffusion.hh
+maxwellstefandiffusionspatialparams.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/implicit/2pnc)
+set(CMAKE_BUILD_TYPE Release)
diff --git a/test/porousmediumflow/2pnc/implicit/test_2pnc_diffusion.input b/test/porousmediumflow/2pnc/implicit/test_2pnc_diffusion.input
new file mode 100644
index 0000000000000000000000000000000000000000..b658cedc1385c9db7beeee5d4cb8e7e06510140a
--- /dev/null
+++ b/test/porousmediumflow/2pnc/implicit/test_2pnc_diffusion.input
@@ -0,0 +1,11 @@
+[TimeLoop]
+DtInitial = 250 # [s]
+TEnd = 1e8 # [s]
+
+[Grid]
+UpperRight = 60 40
+Cells = 50 5
+
+[Problem]
+Name = 2pnc_diffusion
+EnableGravity = false
diff --git a/test/porousmediumflow/2pnc/implicit/test_cc2pnc_diffusion.cc b/test/porousmediumflow/2pnc/implicit/test_cc2pnc_diffusion.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c039c680287c658e8009484ca46509f808e7359c
--- /dev/null
+++ b/test/porousmediumflow/2pnc/implicit/test_cc2pnc_diffusion.cc
@@ -0,0 +1,243 @@
+// -*- 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 . *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Test for the 2pnc cc model used for water management in PEM fuel cells.
+ */
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+
+#include "2pncdiffusion.hh"
+
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+#include
+#include
+
+#include
+
+#include
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ * reading in parameters.
+ *
+ * \param progName The name of the program, that was tried to be started.
+ * \param errorMsg The error message that was issued by the start function.
+ * Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+ if (errorMsg.size() > 0) {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
+ "\t-ParameterFile Parameter file (Input file) \n";
+
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(TwoPNCDiffusionCCProblem);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv, usage);
+
+ // try to create a grid (from the given grid file or the input file)
+ using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+ GridCreator::makeGrid();
+ GridCreator::loadBalance();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = GridCreator::grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam("TimeLoop.TEnd");
+ const auto maxDivisions = getParam("TimeLoop.MaxTimeStepDivisions");
+ const auto maxDt = getParam("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam("TimeLoop.Restart");
+
+ // intialize the vtk output module
+ using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+ VtkOutputModule vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+ VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+ vtkWriter.write(0.0);
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = FVAssembler;
+ auto assembler = std::make_shared(problem, fvGridGeometry, gridVariables, timeLoop);
+
+ // the linear solver
+ using LinearSolver = AMGBackend;
+ auto linearSolver = std::make_shared(leafGridView, fvGridGeometry->dofMapper());
+
+ // the non-linear solver
+ using NewtonController = PriVarSwitchNewtonController;
+ using NewtonMethod = NewtonMethod;
+ auto newtonController = std::make_shared(leafGridView.comm(), timeLoop);
+ NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(xOld);
+
+ // try solving the non-linear system
+ for (int i = 0; i < maxDivisions; ++i)
+ {
+ // linearize & solve
+ auto converged = nonLinearSolver.solve(x);
+
+ if (converged)
+ break;
+
+ if (!converged && i == maxDivisions-1)
+ DUNE_THROW(Dune::MathError,
+ "Newton solver didn't converge after "
+ << maxDivisions
+ << " time-step divisions. dt="
+ << timeLoop->timeStepSize()
+ << ".\nThe solutions of the current and the previous time steps "
+ << "have been saved to restart files.");
+ }
+
+ // make the new solution the old solution
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // write vtk output
+ vtkWriter.write(timeLoop->time());
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by newton controller
+ timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/test/porousmediumflow/tracer/CMakeLists.txt b/test/porousmediumflow/tracer/CMakeLists.txt
index 7770e5721c70a9cfe7886bf13faf30360902c875..213b2b806b7107612a60ce5a10fdda093cf7dc0c 100644
--- a/test/porousmediumflow/tracer/CMakeLists.txt
+++ b/test/porousmediumflow/tracer/CMakeLists.txt
@@ -1,2 +1,3 @@
add_subdirectory("constvel")
add_subdirectory("1ptracer")
+add_subdirectory("multicomp")
diff --git a/test/porousmediumflow/tracer/multicomp/CMakeLists.txt b/test/porousmediumflow/tracer/multicomp/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a7b2e68e70ea0db68461119deff645ba97036dec
--- /dev/null
+++ b/test/porousmediumflow/tracer/multicomp/CMakeLists.txt
@@ -0,0 +1,25 @@
+add_input_file_links()
+
+# tracer tests
+dune_add_test(NAME test_maxwellstefan_tpfa
+ SOURCES test_tracer_maxwellstefan.cc
+ COMPILE_DEFINITIONS TYPETAG=MaxwellStefanTestCCProblem
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/cctracermaxwellstefan-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_ccmaxwellstefan-00026.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_maxwellstefan_tpfa test_tracer_maxwellstefan.input -Problem.Name test_ccmaxwellstefan")
+
+dune_add_test(NAME test_maxwellstefan_box
+ SOURCES test_tracer_maxwellstefan.cc
+ COMPILE_DEFINITIONS TYPETAG=MaxwellStefanTestBoxProblem
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/boxtracermaxwellstefan-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_boxmaxwellstefan-00026.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_maxwellstefan_box test_tracer_maxwellstefan.input -Problem.Name test_boxmaxwellstefan")
+#install sources
+install(FILES
+maxwellstefantestproblem.hh
+maxwellstefantestspatialparams.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/porousmediumflow/tracer)
diff --git a/test/porousmediumflow/tracer/multicomp/maxwellstefantestproblem.hh b/test/porousmediumflow/tracer/multicomp/maxwellstefantestproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7be6fc9431991acf4d8e737a532fd2c6f0f569e4
--- /dev/null
+++ b/test/porousmediumflow/tracer/multicomp/maxwellstefantestproblem.hh
@@ -0,0 +1,400 @@
+// -*- 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 . *
+ *****************************************************************************/
+/**
+ * \file
+ * \ingroup TracerTest
+ * \brief Definition of a problem, for the MaxwellStefan problem:
+ * A rotating velocity field mixes a MaxwellStefan band in a porous groundwater reservoir.
+ */
+#ifndef DUMUX_MAXWELL_STEFAN_TEST_PROBLEM_HH
+#define DUMUX_MAXWELL_STEFAN_TEST_PROBLEM_HH
+
+#include
+#include
+
+#include
+#include
+
+#include "maxwellstefantestspatialparams.hh"
+#include
+
+#include
+#include
+
+namespace Dumux
+{
+template
+class MaxwellStefanTestProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(MaxwellStefanTestProblem, INHERITS_FROM(Tracer));
+NEW_TYPE_TAG(MaxwellStefanTestCCProblem, INHERITS_FROM(CCTpfaModel, MaxwellStefanTestProblem));
+NEW_TYPE_TAG(MaxwellStefanTestBoxProblem, INHERITS_FROM(BoxModel, MaxwellStefanTestProblem));
+
+// Set the grid type
+SET_TYPE_PROP(MaxwellStefanTestProblem, Grid, Dune::YaspGrid<2>);
+
+// Set the problem property
+SET_TYPE_PROP(MaxwellStefanTestProblem, Problem, MaxwellStefanTestProblem);
+
+// Set the spatial parameters
+SET_TYPE_PROP(MaxwellStefanTestProblem, SpatialParams, MaxwellStefanTestSpatialParams);
+
+// Define whether mole(true) or mass (false) fractions are used
+SET_BOOL_PROP(MaxwellStefanTestProblem, UseMoles, true);
+
+//! Here we set FicksLaw or MaxwellStefansLaw
+SET_TYPE_PROP(MaxwellStefanTestProblem, MolecularDiffusionType, MaxwellStefansLaw);
+
+//! A simple fluid system with one MaxwellStefan component
+template
+class H2N2CO2FluidSystem: public FluidSystems::BaseFluidSystem>
+
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+
+public:
+ static constexpr bool isTracerFluidSystem()
+ { return true; }
+ //! The number of components
+ static constexpr int numComponents = 3;
+
+ static constexpr int H2Idx = 0;//first major component
+ static constexpr int N2Idx = 1;//second major component
+ static constexpr int CO2Idx = 2;//secondary component
+
+ //! Human readable component name (index compIdx) (for vtk output)
+ static std::string componentName(int compIdx)
+ {
+ switch (compIdx)
+ {
+ case H2Idx: return "H2";
+ case N2Idx: return "N2";
+ case CO2Idx:return "CO2";
+ }
+ DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
+ }
+
+ //! Human readable phase name (index phaseIdx) (for velocity vtk output)
+ static std::string phaseName(int phaseIdx = 0)
+ { return "Gas"; }
+
+ //! Molar mass in kg/mol of the component with index compIdx
+ static Scalar molarMass(unsigned int compIdx)
+ {
+ switch (compIdx)
+ {
+ case H2Idx: return 0.002;
+ case N2Idx: return 0.028;
+ case CO2Idx:return 0.044;
+ }
+ DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
+ }
+
+ //! binary diffusion coefficient
+ //! (might depend on spatial parameters like pressure / temperature)
+ static Scalar binaryDiffusionCoefficient(unsigned int compIdx,
+ const Problem& problem,
+ const Element& element,
+ const SubControlVolume& scv)
+ {
+ if (compIdx == H2Idx)
+ return 0;
+ if (compIdx == N2Idx)
+ return 83.3e-6;
+ if (compIdx == CO2Idx)
+ return 68.0e-6;
+ DUNE_THROW(Dune::InvalidStateException,
+ "Binary diffusion coefficient of component "
+ << compIdx <<" is undefined!\n");
+ }
+
+ //! binary diffusion coefficient
+ //! (might depend on spatial parameters like pressure / temperature)
+ static Scalar binaryDiffusionCoefficient(unsigned int compIIdx,
+ unsigned int compJIdx,
+ const Problem& problem,
+ const Element& element,
+ const SubControlVolume& scv)
+ {
+ if (compIIdx > compJIdx)
+ {
+ using std::swap;
+ swap(compIIdx, compJIdx);
+ }
+
+ if (compIIdx == H2Idx && compJIdx == N2Idx)
+ return 83.3e-6;
+ if (compIIdx == H2Idx && compJIdx == CO2Idx)
+ return 68.0e-6;
+ if (compIIdx == N2Idx && compJIdx == CO2Idx)
+ return 16.8e-6;
+ DUNE_THROW(Dune::InvalidStateException,
+ "Binary diffusion coefficient of components "
+ << compIIdx << " and " << compJIdx << " is undefined!\n");
+ }
+};
+
+SET_TYPE_PROP(MaxwellStefanTestProblem, FluidSystem, H2N2CO2FluidSystem);
+
+} // end namespace Properties
+
+
+/*!
+ * \ingroup TracerTest
+ * \brief Definition of a problem, for the MaxwellStefan problem *
+ * This problem uses the \ref MaxwellStefanModel model.
+ *
+ * To run the simulation execute the following line in shell:
+ * ./test_boxmaxwellstefan -ParameterFile ./test_boxmaxwellstefan.input or
+ * ./test_ccmaxwellstefan -ParameterFile ./test_ccMaxwellstefan.input
+ */
+template
+class MaxwellStefanTestProblem : public PorousMediumFlowProblem
+{
+ using ParentType = PorousMediumFlowProblem;
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+
+ //! property that defines whether mole or mass fractions are used
+ static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+
+ static const int dimWorld = GridView::dimensionworld;
+ using GlobalPosition = Dune::FieldVector;
+
+public:
+ MaxwellStefanTestProblem(std::shared_ptr fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ {
+ name_ = getParam("Problem.Name");
+
+ // stating in the console whether mole or mass fractions are used
+ if(useMoles)
+ std::cout<<"problem uses mole fractions" << '\n';
+ else
+ std::cout<<"problem uses mass fractions" << '\n';
+
+ plotOutput_ = false;
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const std::string& name() const
+ { return name_; }
+
+ //! Called after every time step
+ void postTimeStep(const SolutionVector& curSol, Scalar time)
+ {
+ if (plotOutput_)
+ {
+ Scalar x_co2_left = 0.0;
+ Scalar x_n2_left = 0.0;
+ Scalar x_co2_right = 0.0;
+ Scalar x_n2_right = 0.0;
+ Scalar x_h2_left = 0.0;
+ Scalar x_h2_right = 0.0;
+ Scalar i = 0.0;
+ Scalar j = 0.0;
+ if (!(time < 0.0))
+ {
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+
+ ElementSolutionVector elemSol(element, curSol, this->fvGridGeometry());
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto& globalPos = scv.dofPosition();
+ VolumeVariables volVars;
+ volVars.update(elemSol, *this, element, scv);
+
+ if (globalPos[0] < 0.5)
+ {
+ x_co2_left += volVars.moleFraction(0,2);
+
+ x_n2_left += volVars.moleFraction(0,1);
+ x_h2_left += volVars.moleFraction(0,0);
+ i +=1;
+ }
+ else
+ {
+ x_co2_right += volVars.moleFraction(0,2);
+ x_n2_right += volVars.moleFraction(0,1);
+ x_h2_right += volVars.moleFraction(0,0);
+ j +=1;
+ }
+
+ }
+ }
+ x_co2_left /= i;
+ x_n2_left /= i;
+ x_h2_left /= i;
+ x_co2_right /= j;
+ x_n2_right /= j;
+ x_h2_right /= j;
+
+ //do a gnuplot
+ x_.push_back(time); // in seconds
+ y_.push_back(x_n2_left);
+ y2_.push_back(x_n2_right);
+ y3_.push_back(x_co2_left);
+ y4_.push_back(x_co2_right);
+ y5_.push_back(x_h2_left);
+ y6_.push_back(x_h2_right);
+
+ gnuplot_.resetPlot();
+ gnuplot_.setXRange(0, std::min(time, 72000.0));
+ gnuplot_.setYRange(0.4, 0.6);
+ gnuplot_.setXlabel("time [s]");
+ gnuplot_.setYlabel("mole fraction mol/mol");
+ gnuplot_.addDataSetToPlot(x_, y_, "N2_left");
+ gnuplot_.addDataSetToPlot(x_, y2_, "N2_right");
+ gnuplot_.plot("mole_fraction_N2");
+
+ gnuplot2_.resetPlot();
+ gnuplot2_.setXRange(0, std::min(time, 72000.0));
+ gnuplot2_.setYRange(0.0, 0.6);
+ gnuplot2_.setXlabel("time [s]");
+ gnuplot2_.setYlabel("mole fraction mol/mol");
+ gnuplot2_.addDataSetToPlot(x_, y3_, "C02_left");
+ gnuplot2_.addDataSetToPlot(x_, y4_, "C02_right");
+ gnuplot2_.plot("mole_fraction_C02");
+
+ gnuplot3_.resetPlot();
+ gnuplot3_.setXRange(0, std::min(time, 72000.0));
+ gnuplot3_.setYRange(0.0, 0.6);
+ gnuplot3_.setXlabel("time [s]");
+ gnuplot3_.setYlabel("mole fraction mol/mol");
+ gnuplot3_.addDataSetToPlot(x_, y5_, "H2_left");
+ gnuplot3_.addDataSetToPlot(x_, y6_, "H2_right");
+ gnuplot3_.plot("mole_fraction_H2");
+ }
+ }
+ }
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ *
+ * \param globalPos The position for which the bc type should be evaluated
+ */
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+ {
+ BoundaryTypes values;
+ values.setAllNeumann();
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann
+ * boundary segment.
+ *
+ * \param globalPos The position for which the bc type should be evaluated
+ * The units must be according to either using mole or mass fractions. (mole/(m^2*s) or kg/(m^2*s))
+ */
+ PrimaryVariables neumannAtPos(const GlobalPosition& globalPos) const
+ { return PrimaryVariables(0.0); }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The position for which the initial condition should be evaluated
+ *
+ * For this method, the \a values parameter stores primary
+ * variables.
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables initialValues(0.0);
+ if (globalPos[0] < 0.5)
+ {
+ initialValues[FluidSystem::H2Idx] = 0.0;
+ initialValues[FluidSystem::N2Idx] = 0.50086;
+ initialValues[FluidSystem::CO2Idx] = 0.49914;
+ }
+ else
+ {
+ initialValues[FluidSystem::H2Idx] = 0.50121;
+ initialValues[FluidSystem::N2Idx] = 0.49879;
+ initialValues[FluidSystem::CO2Idx] = 0.0;
+ }
+ return initialValues;
+ }
+
+ // \}
+
+private:
+ static constexpr Scalar eps_ = 1e-6;
+ std::string name_;
+
+ Dumux::GnuplotInterface gnuplot_;
+ Dumux::GnuplotInterface gnuplot2_;
+ Dumux::GnuplotInterface gnuplot3_;
+
+ std::vector x_;
+ std::vector y_;
+ std::vector y2_;
+ std::vector y3_;
+ std::vector y4_;
+ std::vector y5_;
+ std::vector y6_;
+
+ bool plotOutput_;
+};
+
+} //end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/tracer/multicomp/maxwellstefantestspatialparams.hh b/test/porousmediumflow/tracer/multicomp/maxwellstefantestspatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0eadc6d948ffc41a967ab0308a834d763174ef30
--- /dev/null
+++ b/test/porousmediumflow/tracer/multicomp/maxwellstefantestspatialparams.hh
@@ -0,0 +1,111 @@
+// -*- 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 . *
+ *****************************************************************************/
+/*!
+ * \file
+ *\ingroup TracerTest
+ * \brief Definition of the spatial parameters for the MaxwellStefan problem
+ */
+#ifndef DUMUX_MAXWELL_STEFAN_TEST_SPATIAL_PARAMS_HH
+#define DUMUX_MAXWELL_STEFAN_TEST_SPATIAL_PARAMS_HH
+
+#include
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TracerTest
+ * \brief Definition of the spatial parameters for the MaxwellStefan problem
+ */
+template
+class MaxwellStefanTestSpatialParams : public FVSpatialParamsOneP
+{
+ using ParentType = FVSpatialParamsOneP;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+
+ static const int dimWorld = GridView::dimensionworld;
+ using GlobalPosition = typename Dune::FieldVector;
+
+public:
+
+ MaxwellStefanTestSpatialParams(const Problem& problem)
+ : ParentType(problem) {}
+
+ /*!
+ * \brief Define the porosity \f$\mathrm{[-]}\f$.
+ *
+ * \param globalPos The global position
+ */
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
+ { return 0.4; }
+
+ /*!
+ * \brief Define the dispersivity.
+ *
+ * \param element The finite element
+ * \param scv The sub-control volume
+ * \param elemSol The solution for all dofs of the element
+ */
+ Scalar dispersivity(const Element &element,
+ const SubControlVolume& scv,
+ const ElementSolutionVector& elemSol) const
+ { return 0; }
+
+ //! Fluid properties that are spatial params in the tracer model
+ //! They can possible vary with space but are usually constants
+
+ //! fluid density
+ Scalar fluidDensity(const Element &element,
+ const SubControlVolume& scv) const
+ { return 1; }
+
+ //! fluid molar mass
+ Scalar fluidMolarMass(const Element &element,
+ const SubControlVolume& scv) const
+ { return 0.02896; /*air*/}
+
+ //! velocity field
+ GlobalPosition velocity(const SubControlVolumeFace& scvf) const
+ {
+ GlobalPosition vel(0.0);
+
+ return vel;
+ }
+
+ //! velocity field
+ Scalar volumeFlux(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ return 0;
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/tracer/multicomp/test_tracer_maxwellstefan.cc b/test/porousmediumflow/tracer/multicomp/test_tracer_maxwellstefan.cc
new file mode 100644
index 0000000000000000000000000000000000000000..3d90cd96979c4e6b73fb33028ffe6e677c2e379d
--- /dev/null
+++ b/test/porousmediumflow/tracer/multicomp/test_tracer_maxwellstefan.cc
@@ -0,0 +1,253 @@
+// -*- 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 . *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief test for the two-phase porousmedium flow model
+ */
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+
+#include "maxwellstefantestproblem.hh"
+
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+#include
+#include
+
+#include
+
+#include
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ * reading in parameters.
+ *
+ * \param progName The name of the program, that was tried to be started.
+ * \param errorMsg The error message that was issued by the start function.
+ * Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+ if (errorMsg.size() > 0) {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\n\nThe list of mandatory arguments for this program is:\n"
+ "\t-TimeManager.TEnd End of the simulation [s] \n"
+ "\t-TimeManager.DtInitial Initial timestep size [s] \n"
+ "\t-Grid.LowerLeft Lower left corner coordinates\n"
+ "\t-Grid.UpperRight Upper right corner coordinates\n"
+ "\t-Grid.Cells Number of cells in respective coordinate directions\n"
+ "\t definition in DGF format\n"
+ "\t-SpatialParams.LensLowerLeft coordinates of the lower left corner of the lens [m] \n"
+ "\t-SpatialParams.LensUpperRight coordinates of the upper right corner of the lens [m] \n"
+ "\t-SpatialParams.Permeability Permeability of the domain [m^2] \n"
+ "\t-SpatialParams.PermeabilityLens Permeability of the lens [m^2] \n";
+
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(TYPETAG);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv, usage);
+
+ // try to create a grid (from the given grid file or the input file)
+ using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+ GridCreator::makeGrid();
+ GridCreator::loadBalance();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = GridCreator::grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam("TimeLoop.TEnd");
+ const auto maxDivisions = getParam("TimeLoop.MaxTimeStepDivisions");
+ const auto maxDt = getParam("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam("TimeLoop.Restart");
+
+ // intialize the vtk output module
+ using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+ VtkOutputModule vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+ VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+ vtkWriter.write(0.0);
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared