diff --git a/dumux/CMakeLists.txt b/dumux/CMakeLists.txt
index 1464208a7f77935e6c923a6d23b8faf6b4f0d378..a92131ff57c90861e0de0de7d22028c350ac3dcf 100644
--- a/dumux/CMakeLists.txt
+++ b/dumux/CMakeLists.txt
@@ -3,7 +3,6 @@ add_subdirectory("assembly")
add_subdirectory("common")
add_subdirectory("discretization")
add_subdirectory("freeflow")
-add_subdirectory("implicit")
add_subdirectory("io")
add_subdirectory("linear")
add_subdirectory("material")
diff --git a/dumux/assembly/staggeredlocalassembler.hh b/dumux/assembly/staggeredlocalassembler.hh
index 98278964dd3f2e52e9d5109905fc01a53e13bd66..8912594d50d028c81ca56d4670e00d04aa508603 100644
--- a/dumux/assembly/staggeredlocalassembler.hh
+++ b/dumux/assembly/staggeredlocalassembler.hh
@@ -30,9 +30,16 @@
#include <dumux/common/properties.hh>
#include <dumux/assembly/diffmethod.hh>
-#include <dumux/implicit/staggered/primaryvariables.hh>
#include <dumux/discretization/staggered/facesolution.hh>
+#include <dune/common/version.hh>
+
+#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
+ #include <dune/common/hybridutilities.hh>
+#else
+ #include <dumux/common/intrange.hh>
+#endif
+
namespace Dumux {
/*!
@@ -75,7 +82,6 @@ class StaggeredLocalAssembler<TypeTag,
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
- using PriVarIndices = typename Dumux::PriVarIndices<TypeTag>;
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
@@ -93,8 +99,11 @@ public:
* to the global matrix. The element residual is written into the right hand side.
*/
template<class Assembler>
- static void assembleJacobianAndResidual(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
+ static void assembleJacobianAndResidual(Assembler& assembler,
+ JacobianMatrix& jac,
+ SolutionVector& res,
+ const Element& element,
+ const SolutionVector& curSol)
{
using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
typename DofTypeIndices::CellCenterIdx cellCenterIdx;
@@ -324,7 +333,7 @@ static void dCCdCC_(Assembler& assembler,
auto& curVolVars = getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scvJ);
VolumeVariables origVolVars(curVolVars);
- for(auto pvIdx : PriVarIndices(cellCenterIdx))
+ for(auto pvIdx : priVarIndices_(cellCenterIdx))
{
CellCenterPrimaryVariables priVars(curSol[cellCenterIdx][globalJ]);
@@ -386,7 +395,7 @@ static void dCCdFace_(Assembler& assembler,
auto& faceVars = getFaceVarAccess(gridVariables.curGridFaceVars(), curElemFaceVars, scvfJ);
const auto origFaceVars(faceVars);
- for(auto pvIdx : PriVarIndices(faceIdx))
+ for(auto pvIdx : priVarIndices_(faceIdx))
{
FacePrimaryVariables facePriVars(curSol[faceIdx][globalJ]);
const Scalar eps = numericEpsilon(facePriVars[pvIdx], cellCenterIdx, faceIdx);
@@ -451,7 +460,7 @@ static void dFacedCC_(Assembler& assembler,
auto& curVolVars = getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scvJ);
VolumeVariables origVolVars(curVolVars);
- for(auto pvIdx : PriVarIndices(cellCenterIdx))
+ for(auto pvIdx : priVarIndices_(cellCenterIdx))
{
CellCenterPrimaryVariables priVars(curSol[cellCenterIdx][globalJ]);
@@ -515,7 +524,7 @@ static void dFacedFace_(Assembler& assembler,
auto& faceVars = getFaceVarAccess(gridVariables.curGridFaceVars(), curElemFaceVars, scvf);
const auto origFaceVars(faceVars);
- for(auto pvIdx : PriVarIndices(faceIdx))
+ for(auto pvIdx : priVarIndices_(faceIdx))
{
auto faceSolution = FaceSolution(scvf, curSol[faceIdx], assembler.fvGridGeometry());
@@ -593,9 +602,31 @@ static void updateGlobalJacobian_(SubMatrix& matrix,
}
}
+//! Helper function that returns an iterable range of primary variable indices.
+//! Specialization for cell center dofs.
+static auto priVarIndices_(typename GET_PROP(TypeTag, DofTypeIndices)::CellCenterIdx)
+{
+ constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
+#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
+ return Dune::range(0, numEqCellCenter);
+#else
+ return IntRange(0, numEqCellCenter);
+#endif
+}
-
+//! Helper function that returns an iterable range of primary variable indices.
+//! Specialization for face dofs.
+static auto priVarIndices_(typename GET_PROP(TypeTag, DofTypeIndices)::FaceIdx)
+{
+ constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
+ constexpr auto numEqFace = GET_PROP_VALUE(TypeTag, NumEqFace);
+#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
+ return Dune::range(numEqCellCenter, numEqCellCenter + numEqFace);
+#else
+ return IntRange(numEqCellCenter, numEqCellCenter + numEqFace);
+#endif
+}
private:
template<class T = TypeTag>
diff --git a/dumux/assembly/staggeredlocalresidual.hh b/dumux/assembly/staggeredlocalresidual.hh
index 7be67fc539dab4369c17d41d506dbd1d209ed58b..5c0179a669ffa64cefde07b5a3984395777a34df 100644
--- a/dumux/assembly/staggeredlocalresidual.hh
+++ b/dumux/assembly/staggeredlocalresidual.hh
@@ -24,7 +24,6 @@
#define DUMUX_STAGGERED_LOCAL_RESIDUAL_HH
#include <dumux/common/valgrind.hh>
-#include <dumux/common/capabilities.hh>
#include <dumux/common/timeloop.hh>
namespace Dumux
@@ -53,7 +52,6 @@ class StaggeredLocalResidual
using Implementation = typename GET_PROP_TYPE(TypeTag, LocalResidual);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Element = typename GridView::template Codim<0>::Entity;
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
@@ -104,7 +102,7 @@ public:
/*!
* \brief Compute the local residual, i.e. the deviation of the
- * equations from zero.
+ * equations from zero for a transient problem.
*
* \param element The DUNE Codim<0> entity for which the residual
* ought to be calculated
@@ -128,12 +126,16 @@ public:
const ElementBoundaryTypes &bcTypes,
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
- // reset all terms
- CellCenterResidual residual;
- residual = 0.0;
- // ccStorageTerm_ = 0.0;
+ assert( ( (timeLoop_ && !isStationary()) || (!timeLoop_ && isStationary() ) ) && "no time loop set for storage term evaluation");
+ assert( ( (prevSol_ && !isStationary()) || (!prevSol_ && isStationary() ) ) && "no solution set for storage term evaluation");
+
+ CellCenterResidual residual(0.0);
+
+ if(isStationary())
+ asImp_().evalVolumeTermForCellCenter_(residual, problem, element, fvGeometry, curElemVolVars, curElemFaceVars, bcTypes);
+ else
+ asImp_().evalVolumeTermForCellCenter_(residual, problem, element, fvGeometry, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, bcTypes);
- asImp_().evalVolumeTermForCellCenter_(residual, problem, element, fvGeometry, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, bcTypes);
asImp_().evalFluxesForCellCenter_(residual, problem, element, fvGeometry, curElemVolVars, curElemFaceVars, bcTypes, elemFluxVarsCache);
asImp_().evalBoundaryForCellCenter_(residual, problem, element, fvGeometry, curElemVolVars, curElemFaceVars, bcTypes, elemFluxVarsCache);
@@ -142,7 +144,7 @@ public:
/*!
* \brief Compute the local residual, i.e. the deviation of the
- * equations from zero.
+ * equations from zero for a transient problem.
*
* \param element The DUNE Codim<0> entity for which the residual
* ought to be calculated
@@ -165,12 +167,18 @@ public:
const ElementFaceVariables& prevElemFaceVars,
const ElementFaceVariables& curElemFaceVars,
const ElementBoundaryTypes &bcTypes,
- const ElementFluxVariablesCache& elemFluxVarsCache,
- const bool resizeResidual = false) const
+ const ElementFluxVariablesCache& elemFluxVarsCache) const
{
- FaceResidual residual;
+ assert( ( (timeLoop_ && !isStationary()) || (!timeLoop_ && isStationary() ) ) && "no time loop set for storage term evaluation");
+ assert( ( (prevSol_ && !isStationary()) || (!prevSol_ && isStationary() ) ) && "no solution set for storage term evaluation");
+
+ FaceResidual residual(0.0);
+
+ if(isStationary())
+ asImp_().evalVolumeTermForFace_(residual, problem, element, fvGeometry, scvf, curElemVolVars, curElemFaceVars, bcTypes);
+ else
+ asImp_().evalVolumeTermForFace_(residual, problem, element, fvGeometry, scvf, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, bcTypes);
- asImp_().evalVolumeTermForFace_(residual, problem, element, fvGeometry, scvf, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, bcTypes);
asImp_().evalFluxesForFace_(residual, problem, element, fvGeometry, scvf, curElemVolVars, curElemFaceVars, bcTypes, elemFluxVarsCache);
asImp_().evalBoundaryForFace_(residual, problem, element, fvGeometry, scvf, curElemVolVars, curElemFaceVars, bcTypes, elemFluxVarsCache);
@@ -202,9 +210,9 @@ public:
protected:
- /*!
- * \brief Evaluate the flux terms for cell center dofs
- */
+ /*!
+ * \brief Evaluate the flux terms for cell center dofs
+ */
void evalFluxesForCellCenter_(CellCenterResidual& residual,
const Problem& problem,
const Element& element,
@@ -221,9 +229,9 @@ protected:
}
}
- /*!
- * \brief Evaluate the flux terms for face dofs
- */
+ /*!
+ * \brief Evaluate the flux terms for face dofs
+ */
void evalFluxesForFace_(FaceResidual& residual,
const Problem& problem,
const Element& element,
@@ -238,9 +246,9 @@ protected:
residual += asImp_().computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
}
- /*!
- * \brief Evaluate boundary conditions
- */
+ /*!
+ * \brief Evaluate boundary conditions
+ */
void evalBoundary_(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
@@ -254,20 +262,16 @@ protected:
"a evalBoundary_() method.");
}
- /*!
- * \brief Evaluate the volume term for a cell center dof for a stationary problem
- */
- template<class P = Problem>
- typename std::enable_if<Dumux::Capabilities::isStationary<P>::value, void>::type
- evalVolumeTermForCellCenter_(CellCenterResidual& residual,
- const Problem& problem,
- const Element &element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& prevElemVolVars,
- const ElementVolumeVariables& curElemVolVars,
- const ElementFaceVariables& prevFaceVars,
- const ElementFaceVariables& curFaceVars,
- const ElementBoundaryTypes &bcTypes) const
+ /*!
+ * \brief Evaluate the volume term for a cell center dof for a stationary problem
+ */
+ void evalVolumeTermForCellCenter_(CellCenterResidual& residual,
+ const Problem& problem,
+ const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& curElemVolVars,
+ const ElementFaceVariables& curFaceVars,
+ const ElementBoundaryTypes &bcTypes) const
{
for(auto&& scv : scvs(fvGeometry))
{
@@ -280,21 +284,17 @@ protected:
}
}
- /*!
- * \brief Evaluate the volume term for a face dof for a stationary problem
- */
- template<class P = Problem>
- typename std::enable_if<Dumux::Capabilities::isStationary<P>::value, void>::type
- evalVolumeTermForFace_(FaceResidual& residual,
- const Problem& problem,
- const Element &element,
- const FVElementGeometry& fvGeometry,
- const SubControlVolumeFace& scvf,
- const ElementVolumeVariables& prevElemVolVars,
- const ElementVolumeVariables& curElemVolVars,
- const ElementFaceVariables& prevFaceVars,
- const ElementFaceVariables& curFaceVars,
- const ElementBoundaryTypes &bcTypes) const
+ /*!
+ * \brief Evaluate the volume term for a face dof for a stationary problem
+ */
+ void evalVolumeTermForFace_(FaceResidual& residual,
+ const Problem& problem,
+ const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf,
+ const ElementVolumeVariables& curElemVolVars,
+ const ElementFaceVariables& curFaceVars,
+ const ElementBoundaryTypes &bcTypes) const
{
// the source term:
auto faceSource = asImp_().computeSourceForFace(problem, scvf, curElemVolVars, curFaceVars);
@@ -304,21 +304,22 @@ protected:
residual -= faceSource;
}
- /*!
- * \brief Evaluate the volume term for a cell center dof for a transient problem
- */
- template<class P = Problem>
- typename std::enable_if<!Dumux::Capabilities::isStationary<P>::value, void>::type
- evalVolumeTermForCellCenter_(CellCenterResidual& residual,
- const Problem& problem,
- const Element &element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& prevElemVolVars,
- const ElementVolumeVariables& curElemVolVars,
- const ElementFaceVariables& prevFaceVars,
- const ElementFaceVariables& curFaceVars,
- const ElementBoundaryTypes &bcTypes) const
+ /*!
+ * \brief Evaluate the volume term for a cell center dof for a transient problem
+ */
+ void evalVolumeTermForCellCenter_(CellCenterResidual& residual,
+ const Problem& problem,
+ const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& prevElemVolVars,
+ const ElementVolumeVariables& curElemVolVars,
+ const ElementFaceVariables& prevFaceVars,
+ const ElementFaceVariables& curFaceVars,
+ const ElementBoundaryTypes &bcTypes) const
{
+ assert(timeLoop_ && "no time loop set for storage term evaluation");
+ assert(prevSol_ && "no solution set for storage term evaluation");
+
for(auto&& scv : scvs(fvGeometry))
{
const auto& curVolVars = curElemVolVars[scv];
@@ -354,22 +355,23 @@ protected:
}
}
- /*!
- * \brief Evaluate the volume term for a face dof for a transient problem
- */
- template<class P = Problem>
- typename std::enable_if<!Dumux::Capabilities::isStationary<P>::value, void>::type
- evalVolumeTermForFace_(FaceResidual& residual,
- const Problem& problem,
- const Element &element,
- const FVElementGeometry& fvGeometry,
- const SubControlVolumeFace& scvf,
- const ElementVolumeVariables& prevElemVolVars,
- const ElementVolumeVariables& curElemVolVars,
- const ElementFaceVariables& prevFaceVars,
- const ElementFaceVariables& curFaceVars,
- const ElementBoundaryTypes &bcTypes) const
+ /*!
+ * \brief Evaluate the volume term for a face dof for a transient problem
+ */
+ void evalVolumeTermForFace_(FaceResidual& residual,
+ const Problem& problem,
+ const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf,
+ const ElementVolumeVariables& prevElemVolVars,
+ const ElementVolumeVariables& curElemVolVars,
+ const ElementFaceVariables& prevFaceVars,
+ const ElementFaceVariables& curFaceVars,
+ const ElementBoundaryTypes &bcTypes) const
{
+ assert(timeLoop_ && "no time loop set for storage term evaluation");
+ assert(prevSol_ && "no solution set for storage term evaluation");
+
const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
const auto& curVolVars = curElemVolVars[scv];
const auto& prevVolVars = prevElemVolVars[scv];
diff --git a/dumux/common/capabilities.hh b/dumux/common/capabilities.hh
deleted file mode 100644
index dbcf98988517918bccdc8bda12ed261d5c93686d..0000000000000000000000000000000000000000
--- a/dumux/common/capabilities.hh
+++ /dev/null
@@ -1,49 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \ingroup Common
- * \file
- * \brief Defines capabilities recognized for problems in Dumux
- * \note Specialize the subsequent capabilities for your problem to optimize
- * your program's effeciency.
- */
-#ifndef DUMUX_CAPABILITIES_HH
-#define DUMUX_CAPABILITIES_HH
-
-#include <dune/common/deprecated.hh>
-
-namespace Dumux
-{
-
-namespace Capabilities
-{
-
-//! If a problem is stationary (not time-dependent)
-template <class Problem>
-struct DUNE_DEPRECATED_MSG("isStationary is deprecated and will be removed!") isStationary
-{
- //! by default all problems are instationary
- static const bool value = false;
-};
-
-} // namespace Capabilities
-
-} // namespace Dumux
-
-#endif
diff --git a/dumux/common/staggeredfvproblem.hh b/dumux/common/staggeredfvproblem.hh
index fb13dc50a06077fcfa821feda1c55829d28c2a50..a311c1e5dcf5d443d570fecc88ee3a4a1f1cfaab 100644
--- a/dumux/common/staggeredfvproblem.hh
+++ b/dumux/common/staggeredfvproblem.hh
@@ -46,6 +46,11 @@ class StaggeredFVProblem : public FVProblem<TypeTag>
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
enum {
dim = GridView::dimension,
@@ -112,20 +117,38 @@ public:
{
// let the problem do the dirty work of nailing down
// the initial solution.
- auto initPriVars = asImp_().initial(scv)[cellCenterIdx];
- auto dofIdxGlobal = scv.dofIndex();
- sol[cellCenterIdx][dofIdxGlobal] += initPriVars;
+ auto initPriVars = asImp_().initial(scv);
+ asImp_().applyInititalCellCenterSolution(sol, scv, initPriVars);
}
// loop over faces
for(auto&& scvf : scvfs(fvGeometry))
{
- auto initPriVars = asImp_().initial(scvf)[faceIdx][scvf.directionIndex()];
- sol[faceIdx][scvf.dofIndex()] = initPriVars;
+ auto initPriVars = asImp_().initial(scvf);
+ asImp_().applyInititalFaceSolution(sol, scvf, initPriVars);
}
}
}
+
+ //! Applys the initial cell center solution
+ void applyInititalCellCenterSolution(SolutionVector& sol,
+ const SubControlVolume& scv,
+ const PrimaryVariables& initSol) const
+ {
+ for(auto&& i : priVarIndices_(cellCenterIdx))
+ sol[cellCenterIdx][scv.dofIndex()][i] = initSol[i];
+ }
+
+ //! Applys the initial face solution
+ void applyInititalFaceSolution(SolutionVector& sol,
+ const SubControlVolumeFace& scvf,
+ const PrimaryVariables& initSol) const
+ {
+ for(auto&& i : priVarIndices_(faceIdx))
+ sol[faceIdx][scvf.dofIndex()][i] = initSol[i];
+ }
+
protected:
//! Returns the implementation of the problem (i.e. static polymorphism)
Implementation &asImp_()
@@ -135,6 +158,32 @@ protected:
const Implementation &asImp_() const
{ return *static_cast<const Implementation *>(this); }
+ //! Helper function that returns an iterable range of primary variable indices.
+ //! Specialization for cell center dofs.
+ static auto priVarIndices_(typename GET_PROP(TypeTag, DofTypeIndices)::CellCenterIdx)
+ {
+ constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
+
+#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
+ return Dune::range(0, numEqCellCenter);
+#else
+ return IntRange(0, numEqCellCenter);
+#endif
+ }
+
+ //! Helper function that returns an iterable range of primary variable indices.
+ //! Specialization for face dofs.
+ static auto priVarIndices_(typename GET_PROP(TypeTag, DofTypeIndices)::FaceIdx)
+ {
+ constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
+ constexpr auto numEq = GET_PROP_VALUE(TypeTag, NumEq);
+#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
+ return Dune::range(numEqCellCenter, numEq);
+#else
+ return IntRange(numEqCellCenter, numEq);
+#endif
+ }
+
};
} // end namespace Dumux
diff --git a/dumux/discretization/staggered/freeflow/facevariables.hh b/dumux/discretization/staggered/freeflow/facevariables.hh
index dc09d3d43fff90e03ce5a56a23fc85893c0da2c2..26b07b34dfd48859865a8a48131ff6b22fb29544 100644
--- a/dumux/discretization/staggered/freeflow/facevariables.hh
+++ b/dumux/discretization/staggered/freeflow/facevariables.hh
@@ -43,6 +43,7 @@ class StaggeredFaceVariables
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
static constexpr int dimWorld = GridView::dimensionworld;
static constexpr int numPairs = (dimWorld == 2) ? 2 : 4;
@@ -120,13 +121,13 @@ public:
{
const auto& normalFace = fvGeometry.scvf(scvf.insideScvIdx(), subFaceData.localNormalFaceIdx);
const auto normalDirIdx = normalFace.directionIndex();
- velocityNormalOutside_[i] = problem.dirichlet(element, makeGhostFace(subFaceData.virtualOuterNormalFaceDofPos))[faceIdx][normalDirIdx];
+ velocityNormalOutside_[i] = problem.dirichlet(element, makeGhostFace(subFaceData.virtualOuterNormalFaceDofPos))[Indices::velocity(normalDirIdx)];
}
// treat the velocity parallel to the face
velocityParallel_[i] = hasParallelNeighbor(subFaceData) ?
velocityParallel_[i] = faceSol[subFaceData.outerParallelFaceDofIdx] :
- problem.dirichlet(element, makeGhostFace(subFaceData.virtualOuterParallelFaceDofPos))[faceIdx][scvf.directionIndex()];
+ problem.dirichlet(element, makeGhostFace(subFaceData.virtualOuterParallelFaceDofPos))[Indices::velocity(scvf.directionIndex())];
}
}
diff --git a/dumux/discretization/staggered/freeflow/fickslaw.hh b/dumux/discretization/staggered/freeflow/fickslaw.hh
index 58fe0a71ce8a3f07729bd24bfb796be587943447..85eb80e9631619ca5d48bb2da6e7588656c580cf 100644
--- a/dumux/discretization/staggered/freeflow/fickslaw.hh
+++ b/dumux/discretization/staggered/freeflow/fickslaw.hh
@@ -61,21 +61,16 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::Staggered >
static constexpr int dimWorld = GridView::dimensionworld;
static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+ static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
static_assert(GET_PROP_VALUE(TypeTag, NumPhases) == 1, "Only one phase allowed supported!");
enum {
-
pressureIdx = Indices::pressureIdx,
- velocityIdx = Indices::velocityIdx,
-
- massBalanceIdx = Indices::massBalanceIdx,
- momentumBalanceIdx = Indices::momentumBalanceIdx,
conti0EqIdx = Indices::conti0EqIdx,
mainCompIdx = Indices::mainCompIdx,
replaceCompEqIdx = Indices::replaceCompEqIdx,
- phaseIdx = Indices::phaseIdx
};
public:
diff --git a/dumux/discretization/staggered/freeflow/properties.hh b/dumux/discretization/staggered/freeflow/properties.hh
index e734421301bb0cf2fd50c5936e6c4e2a869c4984..8dbac411bc996e4d47c60042bb89d57935838ad0 100644
--- a/dumux/discretization/staggered/freeflow/properties.hh
+++ b/dumux/discretization/staggered/freeflow/properties.hh
@@ -44,10 +44,6 @@ namespace Properties
//! Type tag for the staggered scheme specialized for free flow.
NEW_TYPE_TAG(StaggeredFreeFlowModel, INHERITS_FROM(StaggeredModel));
-// TODO: Ugly hack. How can this be improved? This is needed, because otherwise the physical model overwrites the properties set here.
-// This requires to include the physical model before the discretization, otherwise the type tag FreeFlow is undefined.
-UNSET_PROP(FreeFlow, NumEqVector);
-
/*!
* \brief Set the number of equations on the faces to 1. We only consider scalar values because the velocity vector
* is normal to the face.
@@ -103,30 +99,6 @@ public:
//! The variables living on the faces
SET_TYPE_PROP(StaggeredFreeFlowModel, FaceVariables, StaggeredFaceVariables<TypeTag>);
-//! A container class used to specify values for boundary/initial conditions
-SET_PROP(StaggeredFreeFlowModel, PrimaryVariables)
-{
-private:
- using CellCenterBoundaryValues = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using FaceBoundaryValues = Dune::FieldVector<typename GET_PROP_TYPE(TypeTag, Scalar),
- GridView::dimension>;
-public:
- using type = StaggeredPrimaryVariables<TypeTag, CellCenterBoundaryValues, FaceBoundaryValues>;
-};
-
-//! A container class used to specify values for sources and Neumann BCs
-SET_PROP(StaggeredFreeFlowModel, NumEqVector)
-{
-private:
- using CellCenterBoundaryValues = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using FaceBoundaryValues = Dune::FieldVector<typename GET_PROP_TYPE(TypeTag, Scalar),
- GridView::dimension>;
-public:
- using type = StaggeredPrimaryVariables<TypeTag, CellCenterBoundaryValues, FaceBoundaryValues>;
-};
-
//! Boundary types at a single degree of freedom
SET_PROP(StaggeredFreeFlowModel, BoundaryTypes)
{
diff --git a/dumux/discretization/staggered/gridvolumevariables.hh b/dumux/discretization/staggered/gridvolumevariables.hh
index b63beb1b1fe5955415ceec5cf85bff028f6570e9..ba7bce363b3d06c95dbc1ff0919aaaccedf7b689 100644
--- a/dumux/discretization/staggered/gridvolumevariables.hh
+++ b/dumux/discretization/staggered/gridvolumevariables.hh
@@ -97,7 +97,7 @@ public:
for(int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
{
if(bcTypes.isDirichlet(eqIdx) || bcTypes.isDirichletCell(eqIdx))
- boundaryPriVars[eqIdx] = problem().dirichlet(element, scvf)[cellCenterIdx][eqIdx];
+ boundaryPriVars[eqIdx] = problem().dirichlet(element, scvf)[eqIdx];
else if(bcTypes.isNeumann(eqIdx) || bcTypes.isOutflow(eqIdx) || bcTypes.isSymmetry())
boundaryPriVars[eqIdx] = sol[cellCenterIdx][scvf.insideScvIdx()][eqIdx];
//TODO: this assumes a zero-gradient for e.g. the pressure on the boundary
diff --git a/dumux/discretization/staggered/properties.hh b/dumux/discretization/staggered/properties.hh
index ceed74b9791fe772e181e0cdce2651d3008305d9..fe45e134f6040d49b7bd58524e3a27977a5d405a 100644
--- a/dumux/discretization/staggered/properties.hh
+++ b/dumux/discretization/staggered/properties.hh
@@ -33,7 +33,6 @@
#include <dumux/discretization/cellcentered/elementboundarytypes.hh>
#include <dumux/assembly/staggeredlocalresidual.hh>
-#include <dumux/implicit/staggered/primaryvariables.hh>
#include <dumux/discretization/cellcentered/subcontrolvolume.hh>
#include <dumux/discretization/staggered/gridvariables.hh>
@@ -250,25 +249,7 @@ public:
using type = BoundaryTypes<numEqCellCenter + numEqFace>;
};
-//! A container class used to specify values for boundary/initial conditions
-SET_PROP(StaggeredModel, PrimaryVariables)
-{
-private:
- using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
- using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
-public:
- using type = StaggeredPrimaryVariables<TypeTag, CellCenterPrimaryVariables, FacePrimaryVariables>;
-};
-//! A container class used to specify values for sources and Neumann BCs
-SET_PROP(StaggeredModel, NumEqVector)
-{
-private:
- using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
- using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
-public:
- using type = StaggeredPrimaryVariables<TypeTag, CellCenterPrimaryVariables, FacePrimaryVariables>;
-};
SET_TYPE_PROP(StaggeredModel, GridVariables, StaggeredGridVariables<TypeTag>);
diff --git a/dumux/freeflow/navierstokes/indices.hh b/dumux/freeflow/navierstokes/indices.hh
index 165ec691c42e2437934796f5d0267e008172319c..5a72bc66d8c08ff08ac0f00034e8ab7e4f226354 100644
--- a/dumux/freeflow/navierstokes/indices.hh
+++ b/dumux/freeflow/navierstokes/indices.hh
@@ -55,11 +55,15 @@ struct NavierStokesCommonIndices
static constexpr int momentumXBalanceIdx = momentumBalanceIdx; //!< Index of the momentum balance equation
static constexpr int momentumYBalanceIdx = momentumBalanceIdx + 1; //!< Index of the momentum balance equation
static constexpr int momentumZBalanceIdx = momentumBalanceIdx + 2; //!< Index of the momentum balance equation
- static constexpr int velocityIdx = momentumBalanceIdx; //!< Index of the velocity in a solution vector
- static constexpr int velocityXIdx = velocityIdx; //!< Index of the velocity in a solution vector
- static constexpr int velocityYIdx = velocityIdx + 1; //!< Index of the velocity in a solution vector
- static constexpr int velocityZIdx = velocityIdx + 2; //!< Index of the velocity in a solution vector
+ static constexpr int velocityXIdx = momentumBalanceIdx; //!< Index of the velocity in a solution vector
+ static constexpr int velocityYIdx = momentumBalanceIdx + 1; //!< Index of the velocity in a solution vector
+ static constexpr int velocityZIdx = momentumBalanceIdx + 2; //!< Index of the velocity in a solution vector
+
+ static constexpr int velocity(int dirIdx) //!< Index of the velocity in a solution vector given a certain dimension
+ {
+ return dirIdx + momentumBalanceIdx;
+ }
};
// \}
diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh
index 5046e29772b33e9dde05a0a845e0d314216b365d..b755600aaedf3bdd9ed0602400be27e91d13a5af 100644
--- a/dumux/freeflow/navierstokes/problem.hh
+++ b/dumux/freeflow/navierstokes/problem.hh
@@ -63,6 +63,10 @@ class NavierStokesProblem : public NavierStokesParentProblem<TypeTag>
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum {
dim = Grid::dimension,
@@ -113,6 +117,18 @@ public:
const GlobalPosition &gravity() const
{ return gravity_; }
+ //! Applys the initial face solution. Specialization for staggered grid discretization.
+ template <class T = TypeTag>
+ typename std::enable_if<GET_PROP_VALUE(T, DiscretizationMethod) == DiscretizationMethods::Staggered, void>::type
+ applyInititalFaceSolution(SolutionVector& sol,
+ const SubControlVolumeFace& scvf,
+ const PrimaryVariables& initSol) const
+ {
+ typename GET_PROP(TypeTag, DofTypeIndices)::FaceIdx faceIdx;
+ const auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
+ sol[faceIdx][scvf.dofIndex()][numEqCellCenter] = initSol[Indices::velocity(scvf.directionIndex())];
+ }
+
// \}
private:
diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
index 06d2c25ba64fafcad867042102aa42e71f757ddf..778405f0ef4bd51e08a9b487ea6f998b5db30869 100644
--- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
+++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
@@ -76,11 +76,7 @@ class NavierStokesFluxVariablesImpl<TypeTag, DiscretizationMethods::Staggered>
dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
- pressureIdx = Indices::pressureIdx,
- velocityIdx = Indices::velocityIdx,
-
massBalanceIdx = Indices::massBalanceIdx,
- momentumBalanceIdx = Indices::momentumBalanceIdx
};
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
diff --git a/dumux/freeflow/navierstokes/staggered/localresidual.hh b/dumux/freeflow/navierstokes/staggered/localresidual.hh
index 75024913c3895c96b62ccd9aa9d164c3151e0a32..9da48e9c6f6d5e09cd200e0450bbb7e1b8d69517 100644
--- a/dumux/freeflow/navierstokes/staggered/localresidual.hh
+++ b/dumux/freeflow/navierstokes/staggered/localresidual.hh
@@ -63,7 +63,6 @@ class NavierStokesResidualImpl<TypeTag, DiscretizationMethods::Staggered>
using Implementation = typename GET_PROP_TYPE(TypeTag, LocalResidual);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Element = typename GridView::template Codim<0>::Entity;
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
@@ -87,13 +86,14 @@ class NavierStokesResidualImpl<TypeTag, DiscretizationMethods::Staggered>
using CellCenterResidual = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using FaceResidual = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+ static constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
+
enum {
// grid and world dimension
dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
pressureIdx = Indices::pressureIdx,
- velocityIdx = Indices::velocityIdx,
massBalanceIdx = Indices::massBalanceIdx,
momentumBalanceIdx = Indices::momentumBalanceIdx
@@ -153,7 +153,16 @@ public:
const ElementFaceVariables& elemFaceVars,
const SubControlVolume &scv) const
{
- return problem.sourceAtPos(scv.center())[cellCenterIdx];
+ CellCenterPrimaryVariables result(0.0);
+
+ // get the values from the problem
+ const auto sourceValues = problem.sourceAtPos(scv.center());
+
+ // copy the respective cell center related values to the result
+ for(int i = 0; i < numEqCellCenter; ++i)
+ result[i] = sourceValues[i];
+
+ return result;
}
@@ -215,7 +224,7 @@ public:
const auto& insideVolVars = elemVolVars[insideScvIdx];
source += problem.gravity()[scvf.directionIndex()] * insideVolVars.density();
- source += problem.sourceAtPos(scvf.center())[faceIdx][scvf.directionIndex()];
+ source += problem.sourceAtPos(scvf.center())[Indices::velocity(scvf.directionIndex())];
return source;
}
@@ -290,7 +299,7 @@ protected:
if(bcTypes.isNeumann(eqIdx))
{
const auto extrusionFactor = 1.0; //TODO: get correct extrusion factor
- boundaryFlux[eqIdx] = problem.neumann(element, fvGeometry, elemVolVars, scvf)[cellCenterIdx][eqIdx]
+ boundaryFlux[eqIdx] = problem.neumann(element, fvGeometry, elemVolVars, scvf)[eqIdx]
* extrusionFactor * scvf.area();
}
}
@@ -320,7 +329,7 @@ protected:
if(bcTypes.isDirichletCell(massBalanceIdx))
{
const auto& insideVolVars = elemVolVars[insideScv];
- residual[pressureIdx] = insideVolVars.pressure() - problem.dirichletAtPos(insideScv.center())[cellCenterIdx][pressureIdx];
+ residual[pressureIdx] = insideVolVars.pressure() - problem.dirichletAtPos(insideScv.center())[pressureIdx];
}
}
@@ -345,9 +354,8 @@ protected:
// set a fixed value for the velocity for Dirichlet boundary conditions
if(bcTypes.isDirichlet(momentumBalanceIdx))
{
- // const Scalar velocity = faceVars.faceVars(scvf.dofIndex()).velocity();
const Scalar velocity = elementFaceVars[scvf].velocitySelf();
- const Scalar dirichletValue = problem.dirichlet(element, scvf)[faceIdx][scvf.directionIndex()];
+ const Scalar dirichletValue = problem.dirichlet(element, scvf)[Indices::velocity(scvf.directionIndex())];
residual = velocity - dirichletValue;
}
@@ -394,14 +402,14 @@ private:
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideVolVars = elemVolVars[insideScvIdx];
- const Scalar deltaP = normalizePressure ? problem.initialAtPos(scvf.center())[cellCenterIdx][pressureIdx] : 0.0;
+ const Scalar deltaP = normalizePressure ? problem.initialAtPos(scvf.center())[pressureIdx] : 0.0;
Scalar result = (insideVolVars.pressure() - deltaP) * scvf.area() * -1.0 * sign(scvf.outerNormalScalar());
// treat outflow BCs
if(scvf.boundary())
{
- const Scalar pressure = problem.dirichlet(element, scvf)[cellCenterIdx][pressureIdx] - deltaP;
+ const Scalar pressure = problem.dirichlet(element, scvf)[pressureIdx] - deltaP;
result += pressure * scvf.area() * sign(scvf.outerNormalScalar());
}
return result;
diff --git a/dumux/freeflow/navierstokesnc/staggered/fluxvariables.hh b/dumux/freeflow/navierstokesnc/staggered/fluxvariables.hh
index bd3eccdb5f45cfe1a7d49f0b0f0a65c512825992..a2834fd481a8883292fef4e8b5a9a995a98cee2f 100644
--- a/dumux/freeflow/navierstokesnc/staggered/fluxvariables.hh
+++ b/dumux/freeflow/navierstokesnc/staggered/fluxvariables.hh
@@ -48,7 +48,6 @@ class NavierStokesNCFluxVariablesImpl<TypeTag, DiscretizationMethods::Staggered>
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Element = typename GridView::template Codim<0>::Entity;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
@@ -68,15 +67,7 @@ class NavierStokesNCFluxVariablesImpl<TypeTag, DiscretizationMethods::Staggered>
using ParentType = NavierStokesFluxVariables<TypeTag>;
- enum {
-
- pressureIdx = Indices::pressureIdx,
- velocityIdx = Indices::velocityIdx,
-
- massBalanceIdx = Indices::massBalanceIdx,
- momentumBalanceIdx = Indices::momentumBalanceIdx,
- conti0EqIdx = Indices::conti0EqIdx,
- };
+ enum { conti0EqIdx = Indices::conti0EqIdx };
public:
diff --git a/dumux/freeflow/navierstokesnc/staggered/localresidual.hh b/dumux/freeflow/navierstokesnc/staggered/localresidual.hh
index d1d585d5f2416296412dc68133fd1b71d99ef7cd..ed122681990f8e882a808dc26fcf31a5c0cb1c6a 100644
--- a/dumux/freeflow/navierstokesnc/staggered/localresidual.hh
+++ b/dumux/freeflow/navierstokesnc/staggered/localresidual.hh
@@ -150,7 +150,7 @@ protected:
{
const auto& insideVolVars = elemVolVars[insideScv];
const Scalar massOrMoleFraction = useMoles ? insideVolVars.moleFraction(phaseIdx, compIdx) : insideVolVars.massFraction(phaseIdx, compIdx);
- residual[eqIdx] = massOrMoleFraction - problem.dirichletAtPos(insideScv.center())[cellCenterIdx][eqIdx];
+ residual[eqIdx] = massOrMoleFraction - problem.dirichletAtPos(insideScv.center())[eqIdx];
}
}
diff --git a/dumux/implicit/CMakeLists.txt b/dumux/implicit/CMakeLists.txt
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/dumux/implicit/staggered/primaryvariables.hh b/dumux/implicit/staggered/primaryvariables.hh
deleted file mode 100644
index 64114a838736355d548bf11fdd2ae95a9b4668a4..0000000000000000000000000000000000000000
--- a/dumux/implicit/staggered/primaryvariables.hh
+++ /dev/null
@@ -1,135 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \brief Primary Variables for the Staggered Grid models
- */
-#ifndef DUMUX_STAGGERED_PRIMARYVARIABLES_HH
-#define DUMUX_STAGGERED_PRIMARYVARIABLES_HH
-
-#include <dune/istl/multitypeblockvector.hh>
-#include <dumux/common/intrange.hh>
-
-namespace Dumux
-{
-
-/*!
- * \ingroup NavierStokesModel
- * \brief This class inherits from DUNE's MultiTypeBlockVector and provides a specific [] operator for convenience
- */
-template<class TypeTag,
- class CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables),
- class FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables)>
-class StaggeredPrimaryVariables : public Dune::MultiTypeBlockVector<CellCenterPrimaryVariables, FacePrimaryVariables>
-{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
- using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
- typename DofTypeIndices::CellCenterIdx cellCenterIdx;
- typename DofTypeIndices::FaceIdx faceIdx;
-
- using ParentType = Dune::MultiTypeBlockVector<CellCenterPrimaryVariables, FacePrimaryVariables>;
-
- static constexpr auto faceOffset = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
-
-
-public:
- StaggeredPrimaryVariables() = default;
-
- // introduce a blocklevel variable so that this class can be used in a Dune::BlockVector
- static constexpr int blocklevel = 1;
-
-
- /*!
- * \brief Constructor to initialize all entries with the same value
- *
- * \param value The value
- */
- StaggeredPrimaryVariables(const Scalar value) noexcept
- {
- (*this)[cellCenterIdx] = value;
- (*this)[faceIdx] = value;
- }
-
- /*!
- * \brief Constructor to initialize the cellcenter and face primary values with given values
- *
- * \param ccPriVars The cellcenter primary variables used for initialization
- * \param facePriVars The face primary variables used for initialization
- */
- StaggeredPrimaryVariables(CellCenterPrimaryVariables&& ccPriVars, FacePrimaryVariables&& facePriVars) noexcept
- {
- (*this)[cellCenterIdx] = std::forward<decltype(ccPriVars)>(ccPriVars);
- (*this)[faceIdx] = std::forward<decltype(facePriVars)>(facePriVars);
- }
-
- /*!
- * \brief Operator overload which allows to automatically access the "right" priVars vector via pvIdx.
- * const version
- * \note: the ParentType (DUNE multitypeblockvector) [] operator has to be visible (using ...)
- *
- * \param pvIdx The global index of the primary variable
- */
- using ParentType::operator [];
- const Scalar& operator [](const unsigned int pvIdx) const
- {
- if(pvIdx < faceOffset)
- return ParentType::operator[](cellCenterIdx)[pvIdx];
- else
- return ParentType::operator[](faceIdx)[pvIdx - faceOffset];
- }
-
- /*!
- * \brief Operator overload which allows to automatically access the "right" priVars vector via pvIdx
- * non-const version
- *
- * \param pvIdx The global index of the primary variable
- */
- Scalar& operator [](const unsigned int pvIdx)
- {
- if(pvIdx < faceOffset)
- return ParentType::operator[](cellCenterIdx)[pvIdx];
- else
- return ParentType::operator[](faceIdx)[pvIdx - faceOffset];
- }
-};
-
-/*!
-* \brief Class which provides two ranges of indices (cc and face)
-* cc: for(auto i : PriVarIndices(cellCenterIdx)) { ... }
-* face: for(auto i : PriVarIndices(faceIdx)) { ... }
-*/
-template<class TypeTag>
-class PriVarIndices : public IntRange
-{
- using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
- using cellCenterIdx = typename DofTypeIndices::CellCenterIdx;
- using faceIdx = typename DofTypeIndices::FaceIdx;
-
- static constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
- static constexpr auto numEqFace = GET_PROP_VALUE(TypeTag, NumEqFace);
-
-public:
- PriVarIndices(cellCenterIdx) : IntRange(0, numEqCellCenter) {}
- PriVarIndices(faceIdx) : IntRange(numEqCellCenter, numEqCellCenter + numEqFace) {}
-};
-
-} // namespace Dumux
-
-#endif
diff --git a/test/freeflow/staggered/CMakeLists.txt b/test/freeflow/staggered/CMakeLists.txt
index 37f17e820269197c5fddb4427f9e2b8786ab65aa..2ebe3a3c9fcb484a134586675b15852e5456d06b 100644
--- a/test/freeflow/staggered/CMakeLists.txt
+++ b/test/freeflow/staggered/CMakeLists.txt
@@ -68,7 +68,7 @@ dune_add_test(NAME test_stokes_donea
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/stokes-donea-reference.vtu
- ${CMAKE_CURRENT_BINARY_DIR}/test_donea-00002.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_donea-00001.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_stokes_donea")
dune_add_test(NAME test_navierstokes_kovasznay
@@ -77,7 +77,7 @@ dune_add_test(NAME test_navierstokes_kovasznay
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_kovasznay-reference.vtu
- ${CMAKE_CURRENT_BINARY_DIR}/test_kovasznay-00002.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_kovasznay-00001.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_navierstokes_kovasznay")
dune_add_test(NAME test_navierstokes_angeli
diff --git a/test/freeflow/staggered/angelitestproblem.hh b/test/freeflow/staggered/angelitestproblem.hh
index 816b62c7569cdc567c591f1e19bbf0fe514de6e4..6f62b4c04ff5302bd0bf55cbd08646cc75450eaa 100644
--- a/test/freeflow/staggered/angelitestproblem.hh
+++ b/test/freeflow/staggered/angelitestproblem.hh
@@ -37,13 +37,6 @@ namespace Dumux
template <class TypeTag>
class AngeliTestProblem;
-namespace Capabilities
-{
- template<class TypeTag>
- struct isStationary<AngeliTestProblem<TypeTag>>
- { static const bool value = false; };
-}
-
namespace Properties
{
NEW_TYPE_TAG(AngeliTestProblem, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokes));
@@ -286,10 +279,10 @@ public:
// treat cell-center dofs
const auto dofIdxCellCenter = scv.dofIndex();
const auto& posCellCenter = scv.dofPosition();
- const auto analyticalSolutionCellCenter = dirichletAtPos(posCellCenter)[cellCenterIdx];
- const auto numericalSolutionCellCenter = curSol[cellCenterIdx][dofIdxCellCenter];
- sumError[cellCenterIdx] += squaredDiff_(analyticalSolutionCellCenter, numericalSolutionCellCenter) * scv.volume();
- sumReference[cellCenterIdx] += analyticalSolutionCellCenter * analyticalSolutionCellCenter * scv.volume();
+ const auto analyticalSolutionCellCenter = dirichletAtPos(posCellCenter)[pressureIdx];
+ const auto numericalSolutionCellCenter = curSol[cellCenterIdx][dofIdxCellCenter][pressureIdx];
+ sumError[pressureIdx] += squaredDiff_(analyticalSolutionCellCenter, numericalSolutionCellCenter) * scv.volume();
+ sumReference[pressureIdx] += analyticalSolutionCellCenter * analyticalSolutionCellCenter * scv.volume();
totalVolume += scv.volume();
// treat face dofs
@@ -297,7 +290,7 @@ public:
{
const int dofIdxFace = scvf.dofIndex();
const int dirIdx = scvf.directionIndex();
- const auto analyticalSolutionFace = dirichletAtPos(scvf.center())[faceIdx][dirIdx];
+ const auto analyticalSolutionFace = dirichletAtPos(scvf.center())[Indices::velocity(dirIdx)];
const auto numericalSolutionFace = curSol[faceIdx][dofIdxFace][momentumBalanceIdx];
directionIndex[dofIdxFace] = dirIdx;
errorVelocity[dofIdxFace] = squaredDiff_(analyticalSolutionFace, numericalSolutionFace);
@@ -309,8 +302,8 @@ public:
}
// get the absolute and relative discrete L2-error for cell-center dofs
- l2NormAbs[cellCenterIdx] = std::sqrt(sumError[cellCenterIdx] / totalVolume);
- l2NormRel[cellCenterIdx] = std::sqrt(sumError[cellCenterIdx] / sumReference[cellCenterIdx]);
+ l2NormAbs[pressureIdx] = std::sqrt(sumError[pressureIdx] / totalVolume);
+ l2NormRel[pressureIdx] = std::sqrt(sumError[pressureIdx] / sumReference[pressureIdx]);
// get the absolute and relative discrete L2-error for face dofs
for(int i = 0; i < numFaceDofs; ++i)
@@ -319,14 +312,14 @@ public:
const auto error = errorVelocity[i];
const auto ref = velocityReference[i];
const auto volume = staggeredVolume[i];
- sumError[faceIdx][dirIdx] += error * volume;
- sumReference[faceIdx][dirIdx] += ref * volume;
+ sumError[Indices::velocity(dirIdx)] += error * volume;
+ sumReference[Indices::velocity(dirIdx)] += ref * volume;
}
for(int dirIdx = 0; dirIdx < dimWorld; ++dirIdx)
{
- l2NormAbs[faceIdx][dirIdx] = std::sqrt(sumError[faceIdx][dirIdx] / totalVolume);
- l2NormRel[faceIdx][dirIdx] = std::sqrt(sumError[faceIdx][dirIdx] / sumReference[faceIdx][dirIdx]);
+ l2NormAbs[Indices::velocity(dirIdx)] = std::sqrt(sumError[Indices::velocity(dirIdx)] / totalVolume);
+ l2NormRel[Indices::velocity(dirIdx)] = std::sqrt(sumError[Indices::velocity(dirIdx)] / sumReference[Indices::velocity(dirIdx)]);
}
return std::make_pair(l2NormAbs, l2NormRel);
}
@@ -392,11 +385,13 @@ public:
const auto faceDofPosition = scvf.center();
const auto dirIdx = scvf.directionIndex();
const auto analyticalSolutionAtFace = analyticalSolution(faceDofPosition, time());
- analyticalVelocityOnFace_[faceDofIdx][dirIdx] = analyticalSolutionAtFace[faceIdx][dirIdx];
+ analyticalVelocityOnFace_[faceDofIdx][dirIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)];
}
analyticalPressure_[ccDofIdx] = analyticalSolutionAtCc[pressureIdx];
- analyticalVelocity_[ccDofIdx] = analyticalSolutionAtCc[faceIdx];
+
+ for(int dirIdx = 0; dirIdx < dim; ++dirIdx)
+ analyticalVelocity_[ccDofIdx][dirIdx] = analyticalSolutionAtCc[Indices::velocity(dirIdx)];
}
}
}
diff --git a/test/freeflow/staggered/channeltestproblem.hh b/test/freeflow/staggered/channeltestproblem.hh
index 101b184bad620b33156b63bde538362244e9a9f1..51c79e4a9cd4514aa2d8200f47cc910d714d031d 100644
--- a/test/freeflow/staggered/channeltestproblem.hh
+++ b/test/freeflow/staggered/channeltestproblem.hh
@@ -38,13 +38,6 @@ namespace Dumux
template <class TypeTag>
class ChannelTestProblem;
-namespace Capabilities
-{
- template<class TypeTag>
- struct isStationary<ChannelTestProblem<TypeTag>>
- { static const bool value = false; };
-}
-
namespace Properties
{
#if !NONISOTHERMAL
diff --git a/test/freeflow/staggered/closedsystemtestproblem.hh b/test/freeflow/staggered/closedsystemtestproblem.hh
index 015564f651af55eebbab6386ac2f45385e034013..36990a0d86bee76b5f2c24757336c1dac5b28e41 100644
--- a/test/freeflow/staggered/closedsystemtestproblem.hh
+++ b/test/freeflow/staggered/closedsystemtestproblem.hh
@@ -37,13 +37,6 @@ namespace Dumux
template <class TypeTag>
class ClosedSystemTestProblem;
-namespace Capabilities
-{
- template<class TypeTag>
- struct isStationary<ClosedSystemTestProblem<TypeTag>>
- { static const bool value = false; };
-}
-
namespace Properties
{
NEW_TYPE_TAG(ClosedSystemTestProblem, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokes));
diff --git a/test/freeflow/staggered/doneatestproblem.hh b/test/freeflow/staggered/doneatestproblem.hh
index d01b83681a50f19b6f8085a99449035a5d689ce3..27a9db66ee6767594e9853922d3a8737d5e4a8ea 100644
--- a/test/freeflow/staggered/doneatestproblem.hh
+++ b/test/freeflow/staggered/doneatestproblem.hh
@@ -40,13 +40,6 @@ namespace Dumux
template <class TypeTag>
class DoneaTestProblem;
-namespace Capabilities
-{
- template<class TypeTag>
- struct isStationary<DoneaTestProblem<TypeTag>>
- { static const bool value = true; };
-}
-
namespace Properties
{
NEW_TYPE_TAG(DoneaTestProblem, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokes));
@@ -303,10 +296,10 @@ public:
// treat cell-center dofs
const auto dofIdxCellCenter = scv.dofIndex();
const auto& posCellCenter = scv.dofPosition();
- const auto analyticalSolutionCellCenter = analyticalSolution(posCellCenter)[cellCenterIdx];
- const auto numericalSolutionCellCenter = curSol[cellCenterIdx][dofIdxCellCenter];
- sumError[cellCenterIdx] += squaredDiff_(analyticalSolutionCellCenter, numericalSolutionCellCenter) * scv.volume();
- sumReference[cellCenterIdx] += analyticalSolutionCellCenter * analyticalSolutionCellCenter * scv.volume();
+ const auto analyticalSolutionCellCenter = analyticalSolution(posCellCenter)[pressureIdx];
+ const auto numericalSolutionCellCenter = curSol[cellCenterIdx][dofIdxCellCenter][pressureIdx];
+ sumError[pressureIdx] += squaredDiff_(analyticalSolutionCellCenter, numericalSolutionCellCenter) * scv.volume();
+ sumReference[pressureIdx] += analyticalSolutionCellCenter * analyticalSolutionCellCenter * scv.volume();
totalVolume += scv.volume();
// treat face dofs
@@ -314,7 +307,7 @@ public:
{
const int dofIdxFace = scvf.dofIndex();
const int dirIdx = scvf.directionIndex();
- const auto analyticalSolutionFace = analyticalSolution(scvf.center())[faceIdx][dirIdx];
+ const auto analyticalSolutionFace = analyticalSolution(scvf.center())[Indices::velocity(dirIdx)];
const auto numericalSolutionFace = curSol[faceIdx][dofIdxFace][momentumBalanceIdx];
directionIndex[dofIdxFace] = dirIdx;
errorVelocity[dofIdxFace] = squaredDiff_(analyticalSolutionFace, numericalSolutionFace);
@@ -326,8 +319,8 @@ public:
}
// get the absolute and relative discrete L2-error for cell-center dofs
- l2NormAbs[cellCenterIdx] = std::sqrt(sumError[cellCenterIdx] / totalVolume);
- l2NormRel[cellCenterIdx] = std::sqrt(sumError[cellCenterIdx] / sumReference[cellCenterIdx]);
+ l2NormAbs[pressureIdx] = std::sqrt(sumError[pressureIdx] / totalVolume);
+ l2NormRel[pressureIdx] = std::sqrt(sumError[pressureIdx] / sumReference[pressureIdx]);
// get the absolute and relative discrete L2-error for face dofs
for(int i = 0; i < numFaceDofs; ++i)
@@ -336,14 +329,14 @@ public:
const auto error = errorVelocity[i];
const auto ref = velocityReference[i];
const auto volume = staggeredVolume[i];
- sumError[faceIdx][dirIdx] += error * volume;
- sumReference[faceIdx][dirIdx] += ref * volume;
+ sumError[Indices::velocity(dirIdx)] += error * volume;
+ sumReference[Indices::velocity(dirIdx)] += ref * volume;
}
for(int dirIdx = 0; dirIdx < dimWorld; ++dirIdx)
{
- l2NormAbs[faceIdx][dirIdx] = std::sqrt(sumError[faceIdx][dirIdx] / totalVolume);
- l2NormRel[faceIdx][dirIdx] = std::sqrt(sumError[faceIdx][dirIdx] / sumReference[faceIdx][dirIdx]);
+ l2NormAbs[Indices::velocity(dirIdx)] = std::sqrt(sumError[Indices::velocity(dirIdx)] / totalVolume);
+ l2NormRel[Indices::velocity(dirIdx)] = std::sqrt(sumError[Indices::velocity(dirIdx)] / sumReference[Indices::velocity(dirIdx)]);
}
return std::make_pair(l2NormAbs, l2NormRel);
}
@@ -390,9 +383,9 @@ private:
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
- auto ccDofIdx = scv.dofIndex();
- auto ccDofPosition = scv.dofPosition();
- auto analyticalSolutionAtCc = analyticalSolution(ccDofPosition);
+ const auto ccDofIdx = scv.dofIndex();
+ const auto ccDofPosition = scv.dofPosition();
+ const auto analyticalSolutionAtCc = analyticalSolution(ccDofPosition);
// velocities on faces
for (auto&& scvf : scvfs(fvGeometry))
@@ -401,11 +394,13 @@ private:
const auto faceDofPosition = scvf.center();
const auto dirIdx = scvf.directionIndex();
const auto analyticalSolutionAtFace = analyticalSolution(faceDofPosition);
- analyticalVelocityOnFace_[faceDofIdx][dirIdx] = analyticalSolutionAtFace[faceIdx][dirIdx];
+ analyticalVelocityOnFace_[faceDofIdx][dirIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)];
}
analyticalPressure_[ccDofIdx] = analyticalSolutionAtCc[pressureIdx];
- analyticalVelocity_[ccDofIdx] = analyticalSolutionAtCc[faceIdx];
+
+ for(int dirIdx = 0; dirIdx < dim; ++dirIdx)
+ analyticalVelocity_[ccDofIdx][dirIdx] = analyticalSolutionAtCc[Indices::velocity(dirIdx)];
}
}
}
diff --git a/test/freeflow/staggered/kovasznaytestproblem.hh b/test/freeflow/staggered/kovasznaytestproblem.hh
index 49415b0fa1fb7b8bc1d225123502e73ca0127984..3b04ff4342296f239fd155e874d7e903fc512d7e 100644
--- a/test/freeflow/staggered/kovasznaytestproblem.hh
+++ b/test/freeflow/staggered/kovasznaytestproblem.hh
@@ -32,22 +32,11 @@
#include <dumux/discretization/staggered/freeflow/properties.hh>
#include <dumux/freeflow/navierstokes/model.hh>
-// solve Navier-Stokes equations
-#define ENABLE_NAVIERSTOKES 1
-
-
namespace Dumux
{
template <class TypeTag>
class KovasznayTestProblem;
-namespace Capabilities
-{
- template<class TypeTag>
- struct isStationary<KovasznayTestProblem<TypeTag>>
- { static const bool value = true; };
-}
-
namespace Properties
{
NEW_TYPE_TAG(KovasznayTestProblem, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokes));
@@ -70,11 +59,7 @@ SET_BOOL_PROP(KovasznayTestProblem, EnableFVGridGeometryCache, true);
SET_BOOL_PROP(KovasznayTestProblem, EnableGridFluxVariablesCache, true);
SET_BOOL_PROP(KovasznayTestProblem, EnableGridVolumeVariablesCache, true);
-#if ENABLE_NAVIERSTOKES
SET_BOOL_PROP(KovasznayTestProblem, EnableInertiaTerms, true);
-#else
-SET_BOOL_PROP(KovasznayTestProblem, EnableInertiaTerms, false);
-#endif
}
/*!
@@ -299,10 +284,10 @@ public:
// treat cell-center dofs
const auto dofIdxCellCenter = scv.dofIndex();
const auto& posCellCenter = scv.dofPosition();
- const auto analyticalSolutionCellCenter = dirichletAtPos(posCellCenter)[cellCenterIdx];
- const auto numericalSolutionCellCenter = curSol[cellCenterIdx][dofIdxCellCenter];
- sumError[cellCenterIdx] += squaredDiff_(analyticalSolutionCellCenter, numericalSolutionCellCenter) * scv.volume();
- sumReference[cellCenterIdx] += analyticalSolutionCellCenter * analyticalSolutionCellCenter * scv.volume();
+ const auto analyticalSolutionCellCenter = dirichletAtPos(posCellCenter)[pressureIdx];
+ const auto numericalSolutionCellCenter = curSol[cellCenterIdx][dofIdxCellCenter][pressureIdx];
+ sumError[pressureIdx] += squaredDiff_(analyticalSolutionCellCenter, numericalSolutionCellCenter) * scv.volume();
+ sumReference[pressureIdx] += analyticalSolutionCellCenter * analyticalSolutionCellCenter * scv.volume();
totalVolume += scv.volume();
// treat face dofs
@@ -310,7 +295,7 @@ public:
{
const int dofIdxFace = scvf.dofIndex();
const int dirIdx = scvf.directionIndex();
- const auto analyticalSolutionFace = dirichletAtPos(scvf.center())[faceIdx][dirIdx];
+ const auto analyticalSolutionFace = dirichletAtPos(scvf.center())[Indices::velocity(dirIdx)];
const auto numericalSolutionFace = curSol[faceIdx][dofIdxFace][momentumBalanceIdx];
directionIndex[dofIdxFace] = dirIdx;
errorVelocity[dofIdxFace] = squaredDiff_(analyticalSolutionFace, numericalSolutionFace);
@@ -322,8 +307,8 @@ public:
}
// get the absolute and relative discrete L2-error for cell-center dofs
- l2NormAbs[cellCenterIdx] = std::sqrt(sumError[cellCenterIdx] / totalVolume);
- l2NormRel[cellCenterIdx] = std::sqrt(sumError[cellCenterIdx] / sumReference[cellCenterIdx]);
+ l2NormAbs[pressureIdx] = std::sqrt(sumError[pressureIdx] / totalVolume);
+ l2NormRel[pressureIdx] = std::sqrt(sumError[pressureIdx] / sumReference[pressureIdx]);
// get the absolute and relative discrete L2-error for face dofs
for(int i = 0; i < numFaceDofs; ++i)
@@ -332,14 +317,14 @@ public:
const auto error = errorVelocity[i];
const auto ref = velocityReference[i];
const auto volume = staggeredVolume[i];
- sumError[faceIdx][dirIdx] += error * volume;
- sumReference[faceIdx][dirIdx] += ref * volume;
+ sumError[Indices::velocity(dirIdx)] += error * volume;
+ sumReference[Indices::velocity(dirIdx)] += ref * volume;
}
for(int dirIdx = 0; dirIdx < dimWorld; ++dirIdx)
{
- l2NormAbs[faceIdx][dirIdx] = std::sqrt(sumError[faceIdx][dirIdx] / totalVolume);
- l2NormRel[faceIdx][dirIdx] = std::sqrt(sumError[faceIdx][dirIdx] / sumReference[faceIdx][dirIdx]);
+ l2NormAbs[Indices::velocity(dirIdx)] = std::sqrt(sumError[Indices::velocity(dirIdx)] / totalVolume);
+ l2NormRel[Indices::velocity(dirIdx)] = std::sqrt(sumError[Indices::velocity(dirIdx)] / sumReference[Indices::velocity(dirIdx)]);
}
return std::make_pair(l2NormAbs, l2NormRel);
}
@@ -396,11 +381,13 @@ private:
const auto faceDofPosition = scvf.center();
const auto dirIdx = scvf.directionIndex();
const auto analyticalSolutionAtFace = analyticalSolution(faceDofPosition);
- analyticalVelocityOnFace_[faceDofIdx][dirIdx] = analyticalSolutionAtFace[faceIdx][dirIdx];
+ analyticalVelocityOnFace_[faceDofIdx][dirIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)];
}
analyticalPressure_[ccDofIdx] = analyticalSolutionAtCc[pressureIdx];
- analyticalVelocity_[ccDofIdx] = analyticalSolutionAtCc[faceIdx];
+
+ for(int dirIdx = 0; dirIdx < dim; ++dirIdx)
+ analyticalVelocity_[ccDofIdx][dirIdx] = analyticalSolutionAtCc[Indices::velocity(dirIdx)];
}
}
}
diff --git a/test/freeflow/staggered/test_donea.cc b/test/freeflow/staggered/test_donea.cc
index 173b09a2c3c3a2a9c61164fdf3fb77650b614042..4ef7f924d94c42fdc43eff9dd8292ba3e7e4f913 100644
--- a/test/freeflow/staggered/test_donea.cc
+++ b/test/freeflow/staggered/test_donea.cc
@@ -120,7 +120,7 @@ int main(int argc, char** argv) try
auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
fvGridGeometry->update();
- // the problem (initial and boundary conditions)
+ // the problem (boundary conditions)
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
auto problem = std::make_shared<Problem>(fvGridGeometry);
@@ -134,25 +134,11 @@ int main(int argc, char** argv) try
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<GridVariables>(problem, fvGridGeometry);
- gridVariables->init(x, xOld);
-
- // get some time loop parameters
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
- const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
- const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
- auto dt = getParam<Scalar>("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<Scalar>("TimeLoop.Restart");
+ gridVariables->init(x);
// intialize the vtk output module
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
@@ -163,13 +149,9 @@ int main(int argc, char** argv) try
vtkWriter.addFaceField(problem->getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
vtkWriter.write(0.0);
- // instantiate time loop
- auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
- timeLoop->setMaxTimeStepSize(maxDt);
-
// the assembler with time loop for instationary problem
using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
- auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables);
// the linear solver
using LinearSolver = Dumux::UMFPackBackend<TypeTag>;
@@ -178,55 +160,22 @@ int main(int argc, char** argv) try
// the non-linear solver
using NewtonController = StaggeredNewtonController<TypeTag>;
using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
- auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
+ auto newtonController = std::make_shared<NewtonController>(leafGridView.comm());
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();
-
- problem->postTimeStep(x);
-
- // write vtk output
- vtkWriter.write(timeLoop->time());
-
- // report statistics of this time step
- timeLoop->reportTimeStep();
+ // linearize & solve
+ Dune::Timer timer;
+ nonLinearSolver.solve(x);
- // set new dt as suggested by newton controller
- timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
+ // write vtk output
+ vtkWriter.write(1.0);
- } while (!timeLoop->finished());
+ timer.stop();
- timeLoop->finalize(leafGridView.comm());
+ const auto& comm = Dune::MPIHelper::getCollectiveCommunication();
+ std::cout << "Simulation took " << timer.elapsed() << " seconds on "
+ << comm.size() << " processes.\n"
+ << "The cumulative CPU time was " << timer.elapsed()*comm.size() << " seconds.\n";
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
diff --git a/test/freeflow/staggered/test_kovasznay.cc b/test/freeflow/staggered/test_kovasznay.cc
index b1c855a136d9605751c732f73eb91d507b87fd47..50afa4e0a8d1ffddf64ce0be16dca47e9ab86aa9 100644
--- a/test/freeflow/staggered/test_kovasznay.cc
+++ b/test/freeflow/staggered/test_kovasznay.cc
@@ -19,7 +19,7 @@
/*!
* \file
*
- * \brief Test for the staggered grid Stokes model (Donea et al., 2003)
+ * \brief Stationary test for the staggered grid Navier-Stokes model (Kovasznay 1947)
*/
#include <config.h>
@@ -119,7 +119,7 @@ int main(int argc, char** argv) try
auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
fvGridGeometry->update();
- // the problem (initial and boundary conditions)
+ // the problem (boundary conditions)
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
auto problem = std::make_shared<Problem>(fvGridGeometry);
@@ -133,25 +133,11 @@ int main(int argc, char** argv) try
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<GridVariables>(problem, fvGridGeometry);
- gridVariables->init(x, xOld);
-
- // get some time loop parameters
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
- const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
- const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
- auto dt = getParam<Scalar>("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<Scalar>("TimeLoop.Restart");
+ gridVariables->init(x);
// intialize the vtk output module
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
@@ -162,13 +148,9 @@ int main(int argc, char** argv) try
vtkWriter.addFaceField(problem->getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
vtkWriter.write(0.0);
- // instantiate time loop
- auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
- timeLoop->setMaxTimeStepSize(maxDt);
-
// the assembler with time loop for instationary problem
using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
- auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables);
// the linear solver
using LinearSolver = Dumux::UMFPackBackend<TypeTag>;
@@ -177,55 +159,23 @@ int main(int argc, char** argv) try
// the non-linear solver
using NewtonController = StaggeredNewtonController<TypeTag>;
using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
- auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
+ auto newtonController = std::make_shared<NewtonController>(leafGridView.comm());
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();
-
- problem->postTimeStep(x);
-
- // write vtk output
- vtkWriter.write(timeLoop->time());
+ // linearize & solve
+ Dune::Timer timer;
+ nonLinearSolver.solve(x);
- // report statistics of this time step
- timeLoop->reportTimeStep();
+ // write vtk output
+ vtkWriter.write(1.0);
- // set new dt as suggested by newton controller
- timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
+ timer.stop();
- } while (!timeLoop->finished());
+ const auto& comm = Dune::MPIHelper::getCollectiveCommunication();
+ std::cout << "Simulation took " << timer.elapsed() << " seconds on "
+ << comm.size() << " processes.\n"
+ << "The cumulative CPU time was " << timer.elapsed()*comm.size() << " seconds.\n";
- timeLoop->finalize(leafGridView.comm());
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
diff --git a/test/freeflow/staggered/test_navierstokes_kovasznay.input b/test/freeflow/staggered/test_navierstokes_kovasznay.input
index 1bcc53d179651fdccc276bb87a722200ea1a70ea..03029c6f4623315826bb61f7d797442e77a78394 100644
--- a/test/freeflow/staggered/test_navierstokes_kovasznay.input
+++ b/test/freeflow/staggered/test_navierstokes_kovasznay.input
@@ -1,7 +1,3 @@
-[TimeLoop]
-DtInitial = 1 # [s]
-TEnd = 2 # [s]
-
[Grid]
LowerLeft = -0.5 -0.5
UpperRight = 2 1.5
diff --git a/test/freeflow/staggered/test_stokes_donea.input b/test/freeflow/staggered/test_stokes_donea.input
index d58cec3d52810dafda5446a0cef97ff6e0f337d8..b394c975db1e408ed2bc16acf4e312c076230d54 100644
--- a/test/freeflow/staggered/test_stokes_donea.input
+++ b/test/freeflow/staggered/test_stokes_donea.input
@@ -1,7 +1,3 @@
-[TimeLoop]
-DtInitial = 1 # [s]
-TEnd = 2 # [s]
-
[Grid]
UpperRight = 1 1
Cells = 40 40
diff --git a/test/freeflow/staggerednc/channeltestproblem.hh b/test/freeflow/staggerednc/channeltestproblem.hh
index f7e3d6e0a08955a04d42361819921cc9fe2bd411..051e0c789f71483ad6f093cb840826073e83928f 100644
--- a/test/freeflow/staggerednc/channeltestproblem.hh
+++ b/test/freeflow/staggerednc/channeltestproblem.hh
@@ -39,13 +39,6 @@ namespace Dumux
template <class TypeTag>
class ChannelNCTestProblem;
-namespace Capabilities
-{
- template<class TypeTag>
- struct isStationary<ChannelNCTestProblem<TypeTag>>
- { static const bool value = false; };
-}
-
namespace Properties
{
diff --git a/test/freeflow/staggerednc/densityflowproblem.hh b/test/freeflow/staggerednc/densityflowproblem.hh
index 6d1491b50db05246cfbb68db66dd7c8c2f3395ea..99a028b53fdf686693be40326daf396c8a87ac2a 100644
--- a/test/freeflow/staggerednc/densityflowproblem.hh
+++ b/test/freeflow/staggerednc/densityflowproblem.hh
@@ -39,13 +39,6 @@ namespace Dumux
template <class TypeTag>
class DensityDrivenFlowProblem;
-namespace Capabilities
-{
- template<class TypeTag>
- struct isStationary<DensityDrivenFlowProblem<TypeTag>>
- { static const bool value = false; };
-}
-
namespace Properties
{
NEW_TYPE_TAG(DensityDrivenFlowProblem, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNC));