diff --git a/dumux/freeflow/staggered/localresidual.hh b/dumux/freeflow/staggered/localresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..5045cea0455864883214698db93df741d890ba5f
--- /dev/null
+++ b/dumux/freeflow/staggered/localresidual.hh
@@ -0,0 +1,184 @@
+// -*- 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 Calculates the residual of models based on the box scheme element-wise.
+ */
+#ifndef DUMUX_STAGGERED_NAVIERSTOKES_LOCAL_RESIDUAL_HH
+#define DUMUX_STAGGERED_NAVIERSTOKES_LOCAL_RESIDUAL_HH
+
+#include <dune/istl/matrix.hh>
+
+#include <dumux/common/valgrind.hh>
+#include <dumux/implicit/staggered/localresidual.hh>
+
+#include "properties.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup CCModel
+ * \ingroup StaggeredLocalResidual
+ * \brief Element-wise calculation of the residual for models
+ * based on the fully implicit cell-centered scheme.
+ *
+ * \todo Please doc me more!
+ */
+template<class TypeTag>
+class StaggeredNavierStokesResidual : public Dumux::StaggeredLocalResidual<TypeTag>
+{
+ using ParentType = StaggeredLocalResidual<TypeTag>;
+ friend class StaggeredLocalResidual<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ 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);
+ using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using CellCenterSolutionVector = typename GET_PROP_TYPE(TypeTag, CellCenterSolutionVector);
+ using FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
+ using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+ using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+
+
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+ enum {
+ // grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
+
+public:
+ // copying the local residual class is not a good idea
+ StaggeredNavierStokesResidual(const StaggeredNavierStokesResidual &) = delete;
+
+ StaggeredNavierStokesResidual() = default;
+
+
+ static CellCenterPrimaryVariables computeFluxForCellCenter(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const GlobalFaceVars& globalFaceVars,
+ const SubControlVolumeFace &scvf,
+ const FluxVariablesCache& fluxVarsCache)
+ {
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto& insideVolVars = elemVolVars[insideScv];
+ const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
+ const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndexSelf()).velocity();
+
+ const Scalar eps = 1e-6;
+
+ for(int direction = 0; direction < dim; ++direction)
+ {
+ if(scvf.unitOuterNormal()[direction] > eps && velocity > 0) // positive x-direction
+ return insideVolVars.density(0) * velocity;
+
+ if(scvf.unitOuterNormal()[direction] < eps && velocity < 0)
+ return outsideVolVars.density(0) * velocity;
+ }
+ return CellCenterPrimaryVariables(0.0);
+ }
+
+ static CellCenterPrimaryVariables computeSourceForCellCenter(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const GlobalFaceVars& globalFaceVars,
+ const SubControlVolume &scv)
+ {
+ return CellCenterPrimaryVariables(0.0);
+ }
+
+
+ /*!
+ * \brief Evaluate the rate of change of all conservation
+ * quantites (e.g. phase mass) within a sub-control
+ * volume of a finite volume element for the immiscible models.
+ * \param scv The sub control volume
+ * \param volVars The current or previous volVars
+ * \note This function should not include the source and sink terms.
+ * \note The volVars can be different to allow computing
+ * the implicit euler time derivative here
+ */
+ static CellCenterPrimaryVariables computeStorageForCellCenter(const SubControlVolume& scv,
+ const VolumeVariables& volVars)
+ {
+ CellCenterPrimaryVariables storage;
+ storage[0] = volVars.density(0);
+ return storage;
+ }
+
+ /*!
+ * \brief Evaluate the rate of change of all conservation
+ * quantites (e.g. phase mass) within a sub-control
+ * volume of a finite volume element for the immiscible models.
+ * \param scvf The sub control volume
+ * \param volVars The current or previous volVars
+ * \note This function should not include the source and sink terms.
+ * \note The volVars can be different to allow computing
+ * the implicit euler time derivative here
+ */
+ static FacePrimaryVariables computeStorageForFace(const SubControlVolumeFace& scvf,
+ const VolumeVariables& volVars,
+ const GlobalFaceVars& globalFaceVars)
+ {
+ FacePrimaryVariables storage;
+ const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndexSelf()).velocity();
+ storage[0] = volVars.density(0) * velocity;
+ return storage;
+ }
+
+ static FacePrimaryVariables computeSourceForFace(const SubControlVolumeFace& scvf,
+ const ElementVolumeVariables& elemVolVars,
+ const GlobalFaceVars& globalFaceVars)
+ {
+ return FacePrimaryVariables(0.0);
+ }
+
+ static FacePrimaryVariables computeFluxForFace(const SubControlVolumeFace& scvf,
+ const ElementVolumeVariables& elemVolVars,
+ const GlobalFaceVars& globalFaceVars)
+ {
+ return FacePrimaryVariables(0.0);
+ }
+
+
+
+};
+
+}
+
+#endif // DUMUX_CC_LOCAL_RESIDUAL_HH
diff --git a/dumux/freeflow/staggered/propertydefaults.hh b/dumux/freeflow/staggered/propertydefaults.hh
index 498bb0e398a642a5abc8faf37272d34242477103..b751c593673db51fdd9cb48c575966a2a896644b 100644
--- a/dumux/freeflow/staggered/propertydefaults.hh
+++ b/dumux/freeflow/staggered/propertydefaults.hh
@@ -33,8 +33,9 @@
#include "volumevariables.hh"
#include "indices.hh"
#include "problem.hh"
+#include "localresidual.hh"
-#include <dumux/porousmediumflow/immiscible/localresidual.hh>
+#include <dumux/implicit/staggered/localresidual.hh>
#include <dumux/material/fluidsystems/gasphase.hh>
#include <dumux/material/fluidsystems/liquidphase.hh>
#include <dumux/material/components/nullcomponent.hh>
@@ -52,7 +53,7 @@ SET_INT_PROP(NavierStokes, NumEq, 1); //!< set the number of equations to 1
SET_INT_PROP(NavierStokes, NumPhases, 1); //!< The number of phases in the 1p model is 1
//! The local residual function
-SET_TYPE_PROP(NavierStokes, LocalResidual, ImmiscibleLocalResidual<TypeTag>);
+SET_TYPE_PROP(NavierStokes, LocalResidual, StaggeredNavierStokesResidual<TypeTag>);
//! the Model property
SET_TYPE_PROP(NavierStokes, Model, NavierStokesModel<TypeTag>);
diff --git a/dumux/implicit/staggered/localjacobian.hh b/dumux/implicit/staggered/localjacobian.hh
index 939de60938f69bce58fbdee27d152b3342c05965..4d11f2582478b08039afb3bf202536f8ff44fd13 100644
--- a/dumux/implicit/staggered/localjacobian.hh
+++ b/dumux/implicit/staggered/localjacobian.hh
@@ -145,18 +145,25 @@ public:
ElementBoundaryTypes elemBcTypes;
elemBcTypes.update(this->problem_(), element, fvGeometry);
+ auto& curGlobalFaceVars = this->model_().curGlobalFaceVars();
+ auto& prevGlobalFaceVars = this->model_().prevGlobalFaceVars();
+
// calculate the local residual
if (isGhost)
{
- this->residual_ = 0.0;
- residual[cellCenterIdx][globalI_] = 0.0;
+// this->residual_ = 0.0;
+// residual[cellCenterIdx][globalI_] = 0.0;
}
else
{
- this->localResidual().eval(element, fvGeometry, prevElemVolVars, curElemVolVars, elemBcTypes, elemFluxVarsCache);
- this->residual_ = this->localResidual().residual();
+ this->localResidual().eval(element, fvGeometry,
+ prevElemVolVars, curElemVolVars,
+ prevGlobalFaceVars, curGlobalFaceVars,
+ elemBcTypes, elemFluxVarsCache);
+// this->localResidual().eval(element, fvGeometry, prevElemVolVars, curElemVolVars, elemBcTypes, elemFluxVarsCache);
+// this->residual_ = this->localResidual().residual();
// store residual in global container as well
- residual[cellCenterIdx][globalI_] = this->localResidual().residual(0);
+// residual[cellCenterIdx][globalI_] = this->localResidual().residual(0);
}
this->model_().updatePVWeights(fvGeometry);
diff --git a/dumux/implicit/staggered/localresidual.hh b/dumux/implicit/staggered/localresidual.hh
index 84e1fd1ff8ebd4a7b979bc8560088ed114482f58..38e371d3929af428a54f313b90538434ccb8da2f 100644
--- a/dumux/implicit/staggered/localresidual.hh
+++ b/dumux/implicit/staggered/localresidual.hh
@@ -26,7 +26,6 @@
#include <dune/istl/matrix.hh>
#include <dumux/common/valgrind.hh>
-#include <dumux/implicit/localresidual.hh>
#include "properties.hh"
@@ -41,7 +40,7 @@ namespace Dumux
* \todo Please doc me more!
*/
template<class TypeTag>
-class StaggeredLocalResidual : public ImplicitLocalResidual<TypeTag>
+class StaggeredLocalResidual
{
using ParentType = ImplicitLocalResidual<TypeTag>;
friend class ImplicitLocalResidual<TypeTag>;
@@ -49,6 +48,9 @@ class StaggeredLocalResidual : public ImplicitLocalResidual<TypeTag>
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ 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);
@@ -59,40 +61,200 @@ class StaggeredLocalResidual : public ImplicitLocalResidual<TypeTag>
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using CellCenterSolutionVector = typename GET_PROP_TYPE(TypeTag, CellCenterSolutionVector);
+ using FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
+ using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
+ using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+ using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+
+
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+ enum {
+ // grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
public:
// copying the local residual class is not a good idea
StaggeredLocalResidual(const StaggeredLocalResidual &) = delete;
- StaggeredLocalResidual() : ParentType() {}
+ StaggeredLocalResidual() = default;
+
+
+ /*!
+ * \brief Initialize the local residual.
+ *
+ * This assumes that all objects of the simulation have been fully
+ * allocated but not necessarily initialized completely.
+ *
+ * \param problem The representation of the physical problem to be
+ * solved.
+ */
+ void init(Problem &problem)
+ { problemPtr_ = &problem; }
+
+
+ /*!
+ * \name User interface
+ * \note The following methods are usually expensive to evaluate
+ * They are useful for outputting residual information.
+ */
+ // \{
+
+ /*!
+ * \brief Compute the local residual, i.e. the deviation of the
+ * equations from zero.
+ *
+ * \param element The DUNE Codim<0> entity for which the residual
+ * ought to be calculated
+ */
+ void eval(const Element &element)
+ {
+ // make sure FVElementGeometry and volume variables are bound to the element
+ auto fvGeometry = localView(this->problem().model().globalFvGeometry());
+ fvGeometry.bind(element);
+
+ auto curElemVolVars = localView(problem().model().curGlobalVolVars());
+ curElemVolVars.bind(element, fvGeometry, problem().model().curSol());
+
+ auto prevElemVolVars = localView(problem().model().prevGlobalVolVars());
+ prevElemVolVars.bindElement(element, fvGeometry, problem().model().prevSol());
+
+ auto elemFluxVarsCache = localView(problem().model().globalFluxVarsCache());
+ elemFluxVarsCache.bindElement(element, fvGeometry, curElemVolVars);
+
+ ElementBoundaryTypes bcTypes;
+ bcTypes.update(problem(), element, fvGeometry);
+
+ auto& curGlobalFaceVars = problem().model().curGlobalFaceVars();
+ auto& prevGlobalFaceVars = problem().model().prevGlobalFaceVars();
+
+
+ asImp_().eval(element, fvGeometry,
+ prevElemVolVars, curElemVolVars,
+ prevGlobalFaceVars, curGlobalFaceVars,
+ bcTypes, elemFluxVarsCache);
+ }
+
+ /*!
+ * \brief Compute the local residual, i.e. the deviation of the
+ * equations from zero.
+ *
+ * \param element The DUNE Codim<0> entity for which the residual
+ * ought to be calculated
+ * \param fvGeometry The finite-volume geometry of the element
+ * \param prevVolVars The volume averaged variables for all
+ * sub-control volumes of the element at the previous
+ * time level
+ * \param curVolVars The volume averaged variables for all
+ * sub-control volumes of the element at the current
+ * time level
+ * \param bcTypes The types of the boundary conditions for all
+ * vertices of the element
+ */
+ void eval(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& prevElemVolVars,
+ const ElementVolumeVariables& curElemVolVars,
+ const GlobalFaceVars& prevFaceVars,
+ const GlobalFaceVars& curFaceVars,
+ const ElementBoundaryTypes &bcTypes,
+ const ElementFluxVariablesCache& elemFluxVarsCache)
+ {
+ // resize the vectors for all terms
+ const auto numScv = fvGeometry.numScv();
+ const auto numScvf = fvGeometry.numScvf();
+ ccResidual_.resize(numScv, 0.0);
+ ccStorageTerm_.resize(numScv, 0.0);
+
+ faceResiduals_.resize(numScvf);
+ faceStorageTerms_.resize(numScvf);
+ for(int i = 0; i < numScvf; ++i)
+ {
+ faceResiduals_[i].resize(1, 0.0);
+ faceStorageTerms_[i].resize(1, 0.0);
+ }
+
+ asImp_().evalVolumeTerms_(element, fvGeometry, prevElemVolVars, curElemVolVars, prevFaceVars, curFaceVars, bcTypes);
+ asImp_().evalFluxes_(element, fvGeometry, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
+ asImp_().evalBoundary_(element, fvGeometry, curElemVolVars, bcTypes, elemFluxVarsCache);
+ }
+
+ /*!
+ * \brief Return the problem we are solving. Only call this after init()!
+ */
+ const Problem& problem() const
+ { return *problemPtr_; }
+
+ /*!
+ * \brief Return the problem we are solving. Only call this after init()!
+ */
+ Problem& problem()
+ { return *problemPtr_; }
+
protected:
void evalFluxes_(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
+ const GlobalFaceVars& faceVars,
const ElementBoundaryTypes& bcTypes,
const ElementFluxVariablesCache& elemFluxVarsCache)
{
- // calculate the mass flux over the scv faces
- for (auto&& scvf : scvfs(fvGeometry))
- {
- this->residual_[0] += this->asImp_().computeFlux_(element, fvGeometry, elemVolVars, scvf, bcTypes, elemFluxVarsCache[scvf]);
- }
+ evalFluxesForCellCenter_(element, fvGeometry, elemVolVars, faceVars, bcTypes, elemFluxVarsCache);
+ evalFluxesForFaces_(element, fvGeometry, elemVolVars, faceVars, bcTypes, elemFluxVarsCache);
}
- PrimaryVariables computeFlux_(const Element &element,
+ void evalFluxesForCellCenter_(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace &scvf,
+ const GlobalFaceVars& faceVars,
const ElementBoundaryTypes& bcTypes,
- const FluxVariablesCache& fluxVarsCache)
+ const ElementFluxVariablesCache& elemFluxVarsCache)
{
- if (!scvf.boundary() /*TODO: || GET_PROP_VALUE(TypeTag, BoundaryReconstruction)*/)
- return this->asImp_().computeFlux(element, fvGeometry, elemVolVars, scvf, fluxVarsCache);
- else
- return PrimaryVariables(0.0);
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ ccResidual_[0] += asImp_().computeFluxForCellCenter_(element, fvGeometry, elemVolVars, faceVars, scvf, bcTypes, elemFluxVarsCache[scvf]);
+ }
+ }
+
+ void evalFluxesForFaces_(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const GlobalFaceVars& globalFaceVars,
+ const ElementBoundaryTypes& bcTypes,
+ const ElementFluxVariablesCache& elemFluxVarsCache)
+ {
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ if(!scvf.boundary())
+ faceResiduals_[scvf.localFaceIdx()][0] += asImp_().computeFluxForFace(scvf, elemVolVars, globalFaceVars);
+ }
+ }
+
+
+ CellCenterPrimaryVariables computeFluxForCellCenter_(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const GlobalFaceVars& faceVars,
+ const SubControlVolumeFace &scvf,
+ const ElementBoundaryTypes& bcTypes,
+ const FluxVariablesCache& fluxVarsCache)
+ {
+// if (!scvf.boundary() /*TODO: || GET_PROP_VALUE(TypeTag, BoundaryReconstruction)*/)
+// return /*this->asImp_().*/asImp_().computeFluxForCellCenter(element, fvGeometry, elemVolVars, faceVars, scvf, fluxVarsCache);
+// else
+// {
+// return CellCenterPrimaryVariables(0.0);
+// }
+
+ return asImp_().computeFluxForCellCenter(element, fvGeometry, elemVolVars, faceVars, scvf, fluxVarsCache);
}
void evalBoundary_(const Element &element,
@@ -105,8 +267,8 @@ protected:
{
if (scvf.boundary())
{
- auto bcTypes = this->problem().boundaryTypes(element, scvf);
- this->residual_[0] += evalBoundaryFluxes_(element, fvGeometry, elemVolVars, scvf, bcTypes, elemFluxVarsCache[scvf]);
+// auto bcTypes = this->problem().boundaryTypes(element, scvf);
+// this->residual_[0] += evalBoundaryFluxes_(element, fvGeometry, elemVolVars, scvf, bcTypes, elemFluxVarsCache[scvf]);
}
}
@@ -225,20 +387,20 @@ protected:
const BoundaryTypes &bcTypes)
{
// temporary vector to store the Dirichlet values
- PrimaryVariables dirichletValues = this->problem().dirichlet(element, scvf);
-
- // get the primary variables
- const auto& priVars = elemVolVars[scvf.insideScvIdx()].priVars();
-
- // set Dirichlet conditions in a strong sense
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
- {
- if (bcTypes.isDirichlet(eqIdx))
- {
- auto pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
- this->residual_[0][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
- }
- }
+// PrimaryVariables dirichletValues = this->problem().dirichlet(element, scvf);
+//
+// // get the primary variables
+// const auto& priVars = elemVolVars[scvf.insideScvIdx()].priVars();
+//
+// // set Dirichlet conditions in a strong sense
+// for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+// {
+// if (bcTypes.isDirichlet(eqIdx))
+// {
+// auto pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
+// this->residual_[0][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
+// }
+// }
}
/*!
@@ -299,6 +461,128 @@ protected:
// restore the pointer to the FVElementGeometry
this->fvElemGeomPtr_ = oldFVGeometryPtr;
}*/
+
+
+ /*!
+ * \brief Add the change the source term for stationary problems
+ * to the local residual of all sub-control volumes of the
+ * current element.
+ */
+ template<class P = Problem>
+ typename std::enable_if<Dumux::Capabilities::isStationary<P>::value, void>::type
+ evalVolumeTerms_(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& prevElemVolVars,
+ const ElementVolumeVariables& curElemVolVars,
+ const GlobalFaceVars& prevFaceVars,
+ const GlobalFaceVars& curFaceVars,
+ const ElementBoundaryTypes &bcTypes)
+ {
+ // evaluate the volume terms (storage + source terms)
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ auto curExtrusionFactor = curElemVolVars[scv].extrusionFactor();
+
+ // subtract the source term from the local rate
+ CellCenterPrimaryVariables source = asImp_().computeSourceForCellCenter(element, fvGeometry, curElemVolVars, curFaceVars, scv);
+ source *= scv.volume()*curExtrusionFactor;
+
+ ccResidual_[0] -= source;
+
+ // now, treat the dofs on the facets:
+ for(auto&& scvf : scvfs(fvGeometry))
+ {
+ // the source term:
+ auto faceSource = asImp_().computeSourceForFace(scvf, curElemVolVars, curFaceVars);
+ faceSource *= 0.5*scv.volume()*curExtrusionFactor;
+ faceResiduals_[scvf.localFaceIdx()][0] -= faceSource;
+ }
+ }
+ }
+
+ /*!
+ * \brief Add the change in the storage terms and the source term
+ * to the local residual of all sub-control volumes of the
+ * current element.
+ */
+ template<class P = Problem>
+ typename std::enable_if<!Dumux::Capabilities::isStationary<P>::value, void>::type
+ evalVolumeTerms_(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& prevElemVolVars,
+ const ElementVolumeVariables& curElemVolVars,
+ const GlobalFaceVars& prevFaceVars,
+ const GlobalFaceVars& curFaceVars,
+ const ElementBoundaryTypes &bcTypes)
+ {
+ // evaluate the volume terms (storage + source terms)
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto& curVolVars = curElemVolVars[scv];
+ const auto& prevVolVars = prevElemVolVars[scv];
+
+ // mass balance within the element. this is the
+ // \f$\frac{m}{\partial t}\f$ term if using implicit
+ // euler as time discretization.
+ //
+ // We might need a more explicit way for
+ // doing the time discretization...
+ auto prevCCStorage = asImp_().computeStorageForCellCenter(scv, prevVolVars);
+ auto curCCStorage = asImp_().computeStorageForCellCenter(scv, curVolVars);
+
+ prevCCStorage *= prevVolVars.extrusionFactor();
+ curCCStorage *= curVolVars.extrusionFactor();
+
+ ccStorageTerm_[0] = std::move(curCCStorage);
+ ccStorageTerm_[0] -= std::move(prevCCStorage);
+ ccStorageTerm_[0] *= scv.volume();
+ ccStorageTerm_[0] /= problem().timeManager().timeStepSize();
+
+ // add the storage term to the residual
+ ccResidual_[0] += ccStorageTerm_[0];
+
+ // subtract the source term from the local rate
+ CellCenterPrimaryVariables source = asImp_().computeSourceForCellCenter(element, fvGeometry, curElemVolVars, curFaceVars, scv);
+ source *= scv.volume()*curVolVars.extrusionFactor();
+
+ ccResidual_[0] -= source;
+
+
+ // now, treat the dofs on the facets:
+ for(auto&& scvf : scvfs(fvGeometry))
+ {
+ // the storage term:
+ auto prevFaceStorage = asImp_().computeStorageForFace(scvf, prevVolVars, prevFaceVars);
+ auto curFaceStorage = asImp_().computeStorageForFace(scvf, curVolVars, prevFaceVars);
+ faceStorageTerms_[scvf.localFaceIdx()][0] = std::move(curFaceStorage);
+ faceStorageTerms_[scvf.localFaceIdx()][0] -= std::move(prevFaceStorage);
+ faceStorageTerms_[scvf.localFaceIdx()][0] *= (scv.volume()/2.0);
+ faceStorageTerms_[scvf.localFaceIdx()][0] /= problem().timeManager().timeStepSize();
+
+ faceResiduals_[scvf.localFaceIdx()][0] += faceStorageTerms_[scvf.localFaceIdx()][0];
+
+ // the source term:
+ auto faceSource = asImp_().computeSourceForFace(scvf, curElemVolVars, curFaceVars);
+ faceSource *= 0.5*scv.volume()*curVolVars.extrusionFactor();
+ faceResiduals_[scvf.localFaceIdx()][0] -= faceSource;
+ }
+ }
+ }
+
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+
+ CellCenterSolutionVector ccResidual_;
+ CellCenterSolutionVector ccStorageTerm_;
+ std::vector<FaceSolutionVector> faceResiduals_;
+ std::vector<FaceSolutionVector> faceStorageTerms_;
+
+private:
+ Problem* problemPtr_;
+
};
}
diff --git a/test/freeflow/staggered/staggeredtestproblem.hh b/test/freeflow/staggered/staggeredtestproblem.hh
index 027fb29af025adff6d2ea9370302d2db2456a00e..56d8789ed0a5fe5a79db69302d47da03cba18152 100644
--- a/test/freeflow/staggered/staggeredtestproblem.hh
+++ b/test/freeflow/staggered/staggeredtestproblem.hh
@@ -43,7 +43,7 @@ namespace Capabilities
{
template<class TypeTag>
struct isStationary<StaggeredTestProblem<TypeTag>>
- { static const bool value = true; };
+ { static const bool value = false; };
}
namespace Properties