From b4cf811da197cd18fce53d9c0e7700eeec449549 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Fri, 18 Nov 2016 16:56:48 +0100
Subject: [PATCH] [staggeredModel] Add more basic infrastructure
*Stencils, assembler, fluxVars, etc..
*Almost everything needs to be reworked!
---
dumux/discretization/staggered/darcyslaw.hh | 132 ++---
dumux/discretization/staggered/stencils.hh | 228 +++++++++
.../staggered/subcontrolvolumeface.hh | 5 +
dumux/freeflow/CMakeLists.txt | 1 +
dumux/freeflow/staggered/CMakeLists.txt | 10 +
dumux/freeflow/staggered/fluxvariables.hh | 466 ++++++++++++++++++
.../freeflow/staggered/fluxvariablescache.hh | 193 ++++++++
dumux/freeflow/staggered/indices.hh | 44 ++
dumux/freeflow/staggered/model.hh | 164 ++++++
dumux/freeflow/staggered/properties.hh | 76 +++
dumux/freeflow/staggered/propertydefaults.hh | 158 ++++++
dumux/freeflow/staggered/volumevariables.hh | 183 +++++++
dumux/implicit/staggered/assembler.hh | 106 ++++
dumux/implicit/staggered/localjacobian.hh | 404 +++++++++++++++
dumux/implicit/staggered/localresidual.hh | 306 ++++++++++++
dumux/implicit/staggered/propertydefaults.hh | 32 +-
test/freeflow/CMakeLists.txt | 1 +
17 files changed, 2443 insertions(+), 66 deletions(-)
create mode 100644 dumux/discretization/staggered/stencils.hh
create mode 100644 dumux/freeflow/staggered/CMakeLists.txt
create mode 100644 dumux/freeflow/staggered/fluxvariables.hh
create mode 100644 dumux/freeflow/staggered/fluxvariablescache.hh
create mode 100644 dumux/freeflow/staggered/indices.hh
create mode 100644 dumux/freeflow/staggered/model.hh
create mode 100644 dumux/freeflow/staggered/properties.hh
create mode 100644 dumux/freeflow/staggered/propertydefaults.hh
create mode 100644 dumux/freeflow/staggered/volumevariables.hh
create mode 100644 dumux/implicit/staggered/assembler.hh
create mode 100644 dumux/implicit/staggered/localjacobian.hh
create mode 100644 dumux/implicit/staggered/localresidual.hh
diff --git a/dumux/discretization/staggered/darcyslaw.hh b/dumux/discretization/staggered/darcyslaw.hh
index 06c5d58d1c..a9d6fdf3aa 100644
--- a/dumux/discretization/staggered/darcyslaw.hh
+++ b/dumux/discretization/staggered/darcyslaw.hh
@@ -80,67 +80,81 @@ public:
int phaseIdx,
const FluxVarsCache& fluxVarsCache)
{
- const auto& tij = fluxVarsCache.tij();
-
- // Get the inside and outside volume variables
- const auto& insideScv = fvGeometry.scv(scvFace.insideScvIdx());
- const auto& insideVolVars = elemVolVars[insideScv];
- const auto& outsideVolVars = elemVolVars[scvFace.outsideScvIdx()];
-
- auto hInside = insideVolVars.pressure(phaseIdx);
- auto hOutside = outsideVolVars.pressure(phaseIdx);
-
- if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
- {
- // do averaging for the density
- const auto rhoInside = insideVolVars.density(phaseIdx);
- const auto rhoOutide = outsideVolVars.density(phaseIdx);
- const auto rho = (rhoInside + rhoOutide)*0.5;
-
- // ask for the gravitational acceleration in the inside neighbor
- const auto xInside = insideScv.center();
- const auto gInside = problem.gravityAtPos(xInside);
-
- hInside -= rho*(gInside*xInside);
-
- // and the outside neighbor
- if (scvFace.boundary())
- {
- const auto xOutside = scvFace.center();
- const auto gOutside = problem.gravityAtPos(xOutside);
- hOutside -= rho*(gOutside*xOutside);
- }
- else
- {
- const auto outsideScvIdx = scvFace.outsideScvIdx();
- // as we assemble fluxes from the neighbor to our element the outside index
- // refers to the scv of our element, so we use the scv method
- const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
- const auto xOutside = outsideScv.center();
- const auto gOutside = problem.gravityAtPos(xOutside);
- hOutside -= rho*(gOutside*xOutside);
- }
- }
-
- return tij*(hInside - hOutside);
+// const auto& tij = fluxVarsCache.tij();
+//
+// // Get the inside and outside volume variables
+// const auto& insideScv = fvGeometry.scv(scvFace.insideScvIdx());
+// const auto& insideVolVars = elemVolVars[insideScv];
+// const auto& outsideVolVars = elemVolVars[scvFace.outsideScvIdx()];
+//
+// auto hInside = insideVolVars.pressure(phaseIdx);
+// auto hOutside = outsideVolVars.pressure(phaseIdx);
+//
+// if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
+// {
+// // do averaging for the density
+// const auto rhoInside = insideVolVars.density(phaseIdx);
+// const auto rhoOutide = outsideVolVars.density(phaseIdx);
+// const auto rho = (rhoInside + rhoOutide)*0.5;
+//
+// // ask for the gravitational acceleration in the inside neighbor
+// const auto xInside = insideScv.center();
+// const auto gInside = problem.gravityAtPos(xInside);
+//
+// hInside -= rho*(gInside*xInside);
+//
+// // and the outside neighbor
+// if (scvFace.boundary())
+// {
+// const auto xOutside = scvFace.center();
+// const auto gOutside = problem.gravityAtPos(xOutside);
+// hOutside -= rho*(gOutside*xOutside);
+// }
+// else
+// {
+// const auto outsideScvIdx = scvFace.outsideScvIdx();
+// // as we assemble fluxes from the neighbor to our element the outside index
+// // refers to the scv of our element, so we use the scv method
+// const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
+// const auto xOutside = outsideScv.center();
+// const auto gOutside = problem.gravityAtPos(xOutside);
+// hOutside -= rho*(gOutside*xOutside);
+// }
+// }
+//
+// const GlobalPosition vector = {1,0};
+// const Scalar angle = std::acos(scvFace.unitOuterNormal() * vector / scvFace.unitOuterNormal().two_norm() / vector.two_norm())*180.0/M_PI;
+// const Scalar threshold = 1e-6;
+//
+// if(std::abs(angle) < threshold)
+ return 1.0;
+// // else if(std::abs(angle-180) < threshold)
+// // return -1.0;
+// else
+// return 0;
+//
+// // if(std::abs(angle) < threshold || std::abs(angle-180) < threshold)
+// // return tij*(hInside - hOutside);
+// // std::cout << "normal " << scvFace.unitOuterNormal() << std::endl;
+//
}
- static Stencil stencil(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const SubControlVolumeFace& scvFace)
- {
- Stencil stencil;
- if (!scvFace.boundary())
- {
- stencil.push_back(scvFace.insideScvIdx());
- stencil.push_back(scvFace.outsideScvIdx());
- }
- else
- stencil.push_back(scvFace.insideScvIdx());
-
- return stencil;
- }
+// static Stencil pressureStencil(const Problem& problem,
+// const Element& element,
+// const FVElementGeometry& fvGeometry,
+// const SubControlVolumeFace& scvFace)
+// {
+// Stencil pressureStencil;
+// if (!scvFace.boundary())
+// {
+// pressureStencil.push_back(scvFace.insideScvIdx());
+// pressureStencil.push_back(scvFace.outsideScvIdx());
+// }
+// else
+// pressureStencil.push_back(scvFace.insideScvIdx());
+//
+// return pressureStencil;
+// }
// The flux variables cache has to be bound to an element prior to flux calculations
// During the binding, the transmissibilities will be computed and stored using the method below.
diff --git a/dumux/discretization/staggered/stencils.hh b/dumux/discretization/staggered/stencils.hh
new file mode 100644
index 0000000000..a51d5bf075
--- /dev/null
+++ b/dumux/discretization/staggered/stencils.hh
@@ -0,0 +1,228 @@
+// -*- 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 Implements the notion of stencils for cell-centered models
+ */
+#ifndef DUMUX_DISCRETIZATION_STAGGERED_STENCILS_HH
+#define DUMUX_DISCRETIZATION_STAGGERED_STENCILS_HH
+
+#include <set>
+#include <dumux/implicit/cellcentered/properties.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \brief Element-related stencils
+ */
+template<class TypeTag>
+class StaggeredCellCenterStencils
+{
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using Element = typename GridView::template Codim<0>::Entity;
+ // TODO a separate stencil class all stencils can derive from?
+ using Stencil = std::vector<IndexType>;
+public:
+ //! Update the stencil. We expect a bound fvGeometry
+ void update(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry)
+ {
+ elementStencil_.clear();
+
+ // loop over sub control faces
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ FluxVariables fluxVars;
+ const auto& centerStencil = fluxVars.computeCellCenterStencil(problem, element, fvGeometry, scvf);
+ elementStencil_.insert(elementStencil_.end(), centerStencil.begin(), centerStencil.end());
+ }
+ // make values in elementstencil unique
+ std::sort(elementStencil_.begin(), elementStencil_.end());
+ elementStencil_.erase(std::unique(elementStencil_.begin(), elementStencil_.end()), elementStencil_.end());
+
+ auto globalI = problem.elementMapper().index(element);
+ neighborStencil_ = elementStencil_;
+
+ // remove the element itself and possible ghost neighbors from the neighbor stencil
+ neighborStencil_.erase(std::remove_if(neighborStencil_.begin(), neighborStencil_.end(),
+ [globalI](int i){ return (i == globalI); }),
+ neighborStencil_.end());
+ }
+
+ //! The full element stencil (all element this element is interacting with)
+ const Stencil& elementStencil() const
+ {
+ return elementStencil_;
+ }
+
+ //! The full element stencil without this element
+ const Stencil& neighborStencil() const
+ {
+ return neighborStencil_;
+ }
+
+private:
+ Stencil elementStencil_;
+ Stencil neighborStencil_;
+};
+
+template<class TypeTag>
+class StaggeredFaceStencils
+{
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using Element = typename GridView::template Codim<0>::Entity;
+ // TODO a separate stencil class all stencils can derive from?
+ using Stencil = std::vector<IndexType>;
+
+public:
+ void update(const Problem& problem,
+ const SubControlVolumeFace& scvf)
+ {
+ faceStencil_.clear();
+ FluxVariables fluxVars;
+
+ const auto& faceStencil = fluxVars.computeFaceStencil(problem, scvf);
+ faceStencil_.insert(faceStencil_.end(), faceStencil.begin(), faceStencil.end());
+ std::sort(faceStencil_.begin(), faceStencil_.end());
+ faceStencil_.erase(std::unique(faceStencil_.begin(), faceStencil_.end()), faceStencil_.end());
+ }
+
+ //! The full face stencil (all dofs this face is interacting with)
+ const Stencil& faceStencil() const
+ {
+ return faceStencil_;
+ }
+
+private:
+ Stencil faceStencil_;
+};
+
+
+/*!
+ * \ingroup CCModel
+ * \brief The global stencil container class
+ */
+template<class TypeTag>
+class StaggeredStencilsVector
+{
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using StencilType = StaggeredCellCenterStencils<TypeTag>;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using Stencil = std::vector<IndexType>;
+
+
+public:
+ void update(Problem& problem)
+ {
+ problemPtr_ = &problem;
+
+ const IndexType numElements = problem.gridView().size(0);
+ const IndexType numFacets = problem.gridView().size(1);
+ const IndexType numBoundaryFacets = problem.gridView().grid().numBoundarySegments();
+ elementStencils_.resize(numElements);
+ faceStencils_.resize(2*numFacets - numBoundaryFacets);
+
+ std::vector<Stencil> fullFaceStencils;
+ fullFaceStencils.resize(numFacets);
+
+ for (const auto& element : elements(problem.gridView()))
+ {
+ auto eIdx = problem.elementMapper().index(element);
+ // restrict the FvGeometry locally and bind to the element
+ auto fvGeometry = localView(problem.model().globalFvGeometry());
+ fvGeometry.bindElement(element); // for TPFA bind element is enough
+ elementStencils_[eIdx].update(problem, element, fvGeometry);
+
+ // loop over sub control faces
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ faceStencils_[scvf.index()].update(problem, scvf);
+
+ const IndexType idx = scvf.dofIndexSelf() - numElements;
+ const auto& faceStencil = faceStencils_[scvf.index()].faceStencil();
+ fullFaceStencils[idx].insert(fullFaceStencils[idx].end(), faceStencil.begin(), faceStencil.end());
+ std::sort(fullFaceStencils[idx].begin(), fullFaceStencils[idx].end());
+ fullFaceStencils[idx].erase(std::unique(fullFaceStencils[idx].begin(), fullFaceStencils[idx].end()), fullFaceStencils[idx].end());
+ }
+ }
+ // TODO: is this a good idea?
+ fullFaceDofStencils_ = std::make_unique<decltype(fullFaceStencils)>(fullFaceStencils);
+ }
+
+
+ //! overload for elements
+ auto& get(const Element& entity) const
+ {
+ return elementStencils_[problemPtr_->elementMapper().index(entity)];
+ }
+
+ //! overload for faces
+ auto& get(const SubControlVolumeFace& scvFace) const
+ {
+ const IndexType numElements = problemPtr_->gridView().size(0);
+ return faceStencils_[scvFace.dofIndexSelf() - numElements];
+ }
+
+ size_t completeFaceDofStencilSize(const int idx) const
+ {
+// const IndexType numElements = problemPtr_->gridView().size(0);
+ return fullFaceDofStencils_.get()[0][idx/*-numElements*/].size();
+ // TODO: why does this not work?
+ }
+
+ void getFullFaceDofStencils(std::unique_ptr<std::vector<Stencil>>& ptr)
+ {
+ std::cout << "before move: ";
+ if(fullFaceDofStencils_)
+ std::cout << "object still there" << std::endl;
+ ptr = std::move(fullFaceDofStencils_);
+
+ std::cout << "after move: ";
+ if(fullFaceDofStencils_)
+ std::cout << "object still there" << std::endl;
+ }
+
+
+private:
+ std::vector<StaggeredCellCenterStencils<TypeTag>> elementStencils_;
+ std::vector<StaggeredFaceStencils<TypeTag>> faceStencils_;
+ std::unique_ptr<std::vector<Stencil>> fullFaceDofStencils_;
+
+
+ const Problem* problemPtr_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/staggered/subcontrolvolumeface.hh b/dumux/discretization/staggered/subcontrolvolumeface.hh
index 9ab065a281..2dd1b61f0b 100644
--- a/dumux/discretization/staggered/subcontrolvolumeface.hh
+++ b/dumux/discretization/staggered/subcontrolvolumeface.hh
@@ -214,6 +214,11 @@ public:
return pairData_[idx];
}
+ auto& pairData() const
+ {
+ return pairData_;
+ }
+
private:
Dune::GeometryType geomType_;
std::vector<GlobalPosition> corners_;
diff --git a/dumux/freeflow/CMakeLists.txt b/dumux/freeflow/CMakeLists.txt
index 2793132c66..aff2016afb 100644
--- a/dumux/freeflow/CMakeLists.txt
+++ b/dumux/freeflow/CMakeLists.txt
@@ -1,3 +1,4 @@
+add_subdirectory("staggered")
add_subdirectory("stokes")
add_subdirectory("stokesnc")
add_subdirectory("stokesncni")
diff --git a/dumux/freeflow/staggered/CMakeLists.txt b/dumux/freeflow/staggered/CMakeLists.txt
new file mode 100644
index 0000000000..8baf812d98
--- /dev/null
+++ b/dumux/freeflow/staggered/CMakeLists.txt
@@ -0,0 +1,10 @@
+
+#install headers
+install(FILES
+indices.hh
+localresidual.hh
+model.hh
+properties.hh
+propertydefaults.hh
+volumevariables.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/1p/implicit)
diff --git a/dumux/freeflow/staggered/fluxvariables.hh b/dumux/freeflow/staggered/fluxvariables.hh
new file mode 100644
index 0000000000..5107f731df
--- /dev/null
+++ b/dumux/freeflow/staggered/fluxvariables.hh
@@ -0,0 +1,466 @@
+// -*- 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 Base class for the flux variables
+ */
+#ifndef DUMUX_FREELOW_IMPLICIT_FLUXVARIABLES_HH
+#define DUMUX_FREELOW_IMPLICIT_FLUXVARIABLES_HH
+
+#include <dumux/implicit/properties.hh>
+#include <dumux/discretization/fluxvariablesbase.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+// forward declaration
+NEW_PROP_TAG(NumPhases);
+NEW_PROP_TAG(EnableAdvection);
+NEW_PROP_TAG(EnableMolecularDiffusion);
+NEW_PROP_TAG(EnableEnergyBalance);
+}
+
+// forward declaration
+template<class TypeTag, bool enableAdvection, bool enableMolecularDiffusion, bool enableEnergyBalance>
+class FreeFlowFluxVariablesImpl;
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief The flux variables class
+ * specializations are provided for combinations of physical processes
+ * \note Not all specializations are currently implemented
+ */
+template<class TypeTag>
+using FreeFlowFluxVariables = FreeFlowFluxVariablesImpl<TypeTag, GET_PROP_VALUE(TypeTag, EnableAdvection),
+ GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion),
+ GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>;
+
+
+// specialization for pure advective flow (e.g. 1p/2p/3p immiscible darcy flow)
+template<class TypeTag>
+class FreeFlowFluxVariablesImpl<TypeTag, true, false, false>
+: public FluxVariablesBase<TypeTag, FreeFlowFluxVariablesImpl<TypeTag, true, false, false>>
+{
+ using ParentType = FluxVariablesBase<TypeTag, FreeFlowFluxVariablesImpl<TypeTag, true, false, false>>;
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Stencil = std::vector<IndexType>;
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+
+public:
+
+ void initAndComputeFluxes(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace &scvFace,
+ const FluxVariablesCache& fluxVarsCache)
+ {
+ ParentType::init(problem, element, fvGeometry, elemVolVars, scvFace, fluxVarsCache);
+ }
+
+ template<typename FunctionType>
+ Scalar advectiveFlux(const int phaseIdx, const FunctionType& upwindFunction)
+ {
+ Scalar flux = AdvectionType::flux(this->problem(),
+ this->element(),
+ this->fvGeometry(),
+ this->elemVolVars(),
+ this->scvFace(),
+ phaseIdx,
+ this->fluxVarsCache());
+
+ const auto& insideScv = this->fvGeometry().scv(this->scvFace().insideScvIdx());
+ const auto& insideVolVars = this->elemVolVars()[insideScv];
+ const auto& outsideVolVars = this->elemVolVars()[this->scvFace().outsideScvIdx()];
+
+ if (std::signbit(flux))
+ return flux*upwindFunction(outsideVolVars, insideVolVars);
+ else
+ return flux*upwindFunction(insideVolVars, outsideVolVars);
+ }
+
+ Stencil computeCellCenterStencil(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvFace)
+ {
+ Stencil cellCenterStencil;
+ if (!scvFace.boundary())
+ {
+ // the cell center dof indices
+ cellCenterStencil.push_back(scvFace.insideScvIdx());
+ cellCenterStencil.push_back(scvFace.outsideScvIdx());
+
+ // the face dof indices
+ cellCenterStencil.push_back(scvFace.dofIndexSelf());
+ }
+ else
+ cellCenterStencil.push_back(scvFace.insideScvIdx());
+
+ return cellCenterStencil;
+ }
+
+ Stencil computeFaceStencil(const Problem& problem,
+ const SubControlVolumeFace& scvFace)
+ {
+ Stencil faceStencil;
+ if (!scvFace.boundary())
+ {
+ // the cell center dof indices normal to the face
+ faceStencil.push_back(scvFace.insideScvIdx());
+ faceStencil.push_back(scvFace.outsideScvIdx());
+
+ // the cell center dof indices parallel to the face
+ for(const auto& data : scvFace.pairData())
+ {
+ faceStencil.push_back(data.outerParallelElementDofIdx);
+ faceStencil.push_back(data.outerParallelFaceDofIdx);
+ faceStencil.push_back(data.normalPair.first);
+ faceStencil.push_back(data.normalPair.second);
+ }
+
+ // the face dof indices
+ faceStencil.push_back(scvFace.dofIndexSelf());
+ faceStencil.push_back(scvFace.dofIndexOpposite());
+ }
+ else
+ faceStencil.push_back(scvFace.dofIndexSelf());
+
+ return faceStencil;
+ }
+};
+
+
+// specialization for isothermal advection molecularDiffusion equations
+template<class TypeTag>
+class FreeFlowFluxVariablesImpl<TypeTag, true, true, false>
+: public FluxVariablesBase<TypeTag, FreeFlowFluxVariablesImpl<TypeTag, true, true, false>>
+{
+ using ParentType = FluxVariablesBase<TypeTag, FreeFlowFluxVariablesImpl<TypeTag, true, true, false>>;
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Stencil = std::vector<IndexType>;
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+
+ using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
+ using MolecularDiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
+
+ enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+
+public:
+
+ void initAndComputeFluxes(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace &scvFace,
+ const FluxVariablesCache& fluxVarsCache)
+ {
+ advFluxCached_.reset();
+ ParentType::init(problem, element, fvGeometry, elemVolVars, scvFace, fluxVarsCache);
+ }
+
+ template<typename FunctionType>
+ Scalar advectiveFlux(const int phaseIdx, const FunctionType& upwindFunction)
+ {
+ if (!advFluxCached_[phaseIdx])
+ {
+
+ advPreFlux_[phaseIdx] = AdvectionType::flux(this->problem(),
+ this->element(),
+ this->fvGeometry(),
+ this->elemVolVars(),
+ this->scvFace(),
+ phaseIdx,
+ this->fluxVarsCache());
+ advFluxCached_.set(phaseIdx, true);
+ }
+
+
+
+ const auto& insideScv = this->fvGeometry().scv(this->scvFace().insideScvIdx());
+ const auto& insideVolVars = this->elemVolVars()[insideScv];
+ const auto& outsideVolVars = this->elemVolVars()[this->scvFace().outsideScvIdx()];
+
+ if (std::signbit(advPreFlux_[phaseIdx]))
+ return advPreFlux_[phaseIdx]*upwindFunction(outsideVolVars, insideVolVars);
+ else
+ return advPreFlux_[phaseIdx]*upwindFunction(insideVolVars, outsideVolVars);
+ }
+
+ Scalar molecularDiffusionFlux(const int phaseIdx, const int compIdx)
+ {
+ Scalar flux = MolecularDiffusionType::flux(this->problem(),
+ this->element(),
+ this->fvGeometry(),
+ this->elemVolVars(),
+ this->scvFace(),
+ phaseIdx, compIdx,
+ this->fluxVarsCache());
+ return flux;
+ }
+
+ Stencil computeStencil(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvFace)
+ {
+ // unifiy advective and diffusive stencil
+ Stencil stencil = AdvectionType::stencil(problem, element, fvGeometry, scvFace);
+ Stencil diffusionStencil = MolecularDiffusionType::stencil(problem, element, fvGeometry, scvFace);
+
+ stencil.insert(stencil.end(), diffusionStencil.begin(), diffusionStencil.end());
+ std::sort(stencil.begin(), stencil.end());
+ stencil.erase(std::unique(stencil.begin(), stencil.end()), stencil.end());
+
+ return stencil;
+ }
+private:
+ //! simple caching if advection flux is used twice with different upwind function
+ std::bitset<numPhases> advFluxCached_;
+ std::array<Scalar, numPhases> advPreFlux_;
+};
+
+// specialization for non-isothermal advective flow (e.g. non-isothermal one-phase darcy equation)
+template<class TypeTag>
+class FreeFlowFluxVariablesImpl<TypeTag, true, false, true>
+: public FluxVariablesBase<TypeTag, FreeFlowFluxVariablesImpl<TypeTag, true, false, true>>
+{
+ using ParentType = FluxVariablesBase<TypeTag, FreeFlowFluxVariablesImpl<TypeTag, true, false, true>>;
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Stencil = std::vector<IndexType>;
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+
+ using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
+ using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
+
+ enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+
+public:
+
+ void initAndComputeFluxes(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace &scvFace,
+ const FluxVariablesCache& fluxVarsCache)
+ {
+ advFluxCached_.reset();
+ ParentType::init(problem, element, fvGeometry, elemVolVars, scvFace, fluxVarsCache);
+ }
+
+ template<typename FunctionType>
+ Scalar advectiveFlux(const int phaseIdx, const FunctionType& upwindFunction)
+ {
+ if (!advFluxCached_[phaseIdx])
+ {
+
+ advPreFlux_[phaseIdx] = AdvectionType::flux(this->problem(),
+ this->element(),
+ this->fvGeometry(),
+ this->elemVolVars(),
+ this->scvFace(),
+ phaseIdx,
+ this->fluxVarsCache());
+ advFluxCached_.set(phaseIdx, true);
+ }
+
+
+
+ const auto& insideScv = this->fvGeometry().scv(this->scvFace().insideScvIdx());
+ const auto& insideVolVars = this->elemVolVars()[insideScv];
+ const auto& outsideVolVars = this->elemVolVars()[this->scvFace().outsideScvIdx()];
+
+ if (std::signbit(advPreFlux_[phaseIdx]))
+ return advPreFlux_[phaseIdx]*upwindFunction(outsideVolVars, insideVolVars);
+ else
+ return advPreFlux_[phaseIdx]*upwindFunction(insideVolVars, outsideVolVars);
+ }
+
+ Scalar heatConductionFlux()
+ {
+ Scalar flux = HeatConductionType::flux(this->problem(),
+ this->element(),
+ this->fvGeometry(),
+ this->elemVolVars(),
+ this->scvFace(),
+ this->fluxVarsCache());
+ return flux;
+ }
+
+ Stencil computeStencil(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvFace)
+ {
+ // unifiy advective and diffusive stencil
+ Stencil stencil = AdvectionType::stencil(problem, element, fvGeometry, scvFace);
+ Stencil energyStencil = HeatConductionType::stencil(problem, element, fvGeometry, scvFace);
+
+ stencil.insert(stencil.end(), energyStencil.begin(), energyStencil.end());
+ std::sort(stencil.begin(), stencil.end());
+ stencil.erase(std::unique(stencil.begin(), stencil.end()), stencil.end());
+
+ return stencil;
+ }
+
+private:
+ //! simple caching if advection flux is used twice with different upwind function
+ std::bitset<numPhases> advFluxCached_;
+ std::array<Scalar, numPhases> advPreFlux_;
+};
+
+// specialization for non-isothermal advection and difussion equations (e.g. non-isothermal three-phase three-component flow)
+template<class TypeTag>
+class FreeFlowFluxVariablesImpl<TypeTag, true, true, true>
+: public FluxVariablesBase<TypeTag, FreeFlowFluxVariablesImpl<TypeTag, true, true, true>>
+{
+ using ParentType = FluxVariablesBase<TypeTag, FreeFlowFluxVariablesImpl<TypeTag, true, true, true>>;
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Stencil = std::vector<IndexType>;
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+
+ using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
+ using MolecularDiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
+ using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
+
+ enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+
+public:
+
+ void initAndComputeFluxes(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace &scvFace,
+ const FluxVariablesCache& fluxVarsCache)
+ {
+ advFluxCached_.reset();
+ ParentType::init(problem, element, fvGeometry, elemVolVars, scvFace, fluxVarsCache);
+ }
+
+ template<typename FunctionType>
+ Scalar advectiveFlux(const int phaseIdx, const FunctionType& upwindFunction)
+ {
+ if (!advFluxCached_[phaseIdx])
+ {
+
+ advPreFlux_[phaseIdx] = AdvectionType::flux(this->problem(),
+ this->element(),
+ this->fvGeometry(),
+ this->elemVolVars(),
+ this->scvFace(),
+ phaseIdx,
+ this->fluxVarsCache());
+ advFluxCached_.set(phaseIdx, true);
+ }
+
+
+
+ const auto& insideScv = this->fvGeometry().scv(this->scvFace().insideScvIdx());
+ const auto& insideVolVars = this->elemVolVars()[insideScv];
+ const auto& outsideVolVars = this->elemVolVars()[this->scvFace().outsideScvIdx()];
+
+ if (std::signbit(advPreFlux_[phaseIdx]))
+ return advPreFlux_[phaseIdx]*upwindFunction(outsideVolVars, insideVolVars);
+ else
+ return advPreFlux_[phaseIdx]*upwindFunction(insideVolVars, outsideVolVars);
+ }
+
+ Scalar molecularDiffusionFlux(const int phaseIdx, const int compIdx)
+ {
+ Scalar flux = MolecularDiffusionType::flux(this->problem(),
+ this->element(),
+ this->fvGeometry(),
+ this->elemVolVars(),
+ this->scvFace(),
+ phaseIdx, compIdx,
+ this->fluxVarsCache());
+ return flux;
+ }
+
+ Scalar heatConductionFlux()
+ {
+ Scalar flux = HeatConductionType::flux(this->problem(),
+ this->element(),
+ this->fvGeometry(),
+ this->elemVolVars(),
+ this->scvFace(),
+ this->fluxVarsCache());
+ return flux;
+ }
+
+ Stencil computeStencil(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvFace)
+ {
+ // unifiy advective and diffusive stencil
+ Stencil stencil = AdvectionType::stencil(problem, element, fvGeometry, scvFace);
+ Stencil diffusionStencil = MolecularDiffusionType::stencil(problem, element, fvGeometry, scvFace);
+ Stencil energyStencil = HeatConductionType::stencil(problem, element, fvGeometry, scvFace);
+
+ stencil.insert(stencil.end(), diffusionStencil.begin(), diffusionStencil.end());
+ stencil.insert(stencil.end(), energyStencil.begin(), energyStencil.end());
+ std::sort(stencil.begin(), stencil.end());
+ stencil.erase(std::unique(stencil.begin(), stencil.end()), stencil.end());
+
+ return stencil;
+ }
+
+private:
+ //! simple caching if advection flux is used twice with different upwind function
+ std::bitset<numPhases> advFluxCached_;
+ std::array<Scalar, numPhases> advPreFlux_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/staggered/fluxvariablescache.hh b/dumux/freeflow/staggered/fluxvariablescache.hh
new file mode 100644
index 0000000000..74aeabf631
--- /dev/null
+++ b/dumux/freeflow/staggered/fluxvariablescache.hh
@@ -0,0 +1,193 @@
+// -*- 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 Base class for the flux variables
+ */
+#ifndef DUMUX_FREEFLOW_IMPLICIT_FLUXVARIABLESCACHE_HH
+#define DUMUX_FREEFLOW_IMPLICIT_FLUXVARIABLESCACHE_HH
+
+#include <dumux/implicit/properties.hh>
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+#include <dumux/discretization/methods.hh>
+
+namespace Dumux
+{
+// forward declaration
+template<class TypeTag, DiscretizationMethods Method>
+class FreeFlowFluxVariablesCacheImplementation
+{};
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief The flux variables cache classes for porous media.
+ * Store flux stencils and data required for flux calculation
+ */
+template<class TypeTag>
+using FreeFlowFluxVariablesCache = FreeFlowFluxVariablesCacheImplementation<TypeTag, GET_PROP_VALUE(TypeTag, DiscretizationMethod)>;
+
+// specialization for the Box Method
+template<class TypeTag>
+class FreeFlowFluxVariablesCacheImplementation<TypeTag, DiscretizationMethods::Box>
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using Stencil = std::vector<IndexType>;
+ using TransmissibilityVector = std::vector<IndexType>;
+
+ using CoordScalar = typename GridView::ctype;
+ static const int dim = GridView::dimension;
+
+ using FeCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
+ using FeLocalBasis = typename FeCache::FiniteElementType::Traits::LocalBasisType;
+ using ShapeJacobian = typename FeLocalBasis::Traits::JacobianType;
+ using ShapeValue = typename Dune::FieldVector<Scalar, 1>;
+ using JacobianInverseTransposed = typename Element::Geometry::JacobianInverseTransposed;
+
+public:
+
+ void update(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace &scvf)
+ {
+ const auto geometry = element.geometry();
+ const auto& localBasis = fvGeometry.feLocalBasis();
+
+ // evaluate shape functions and gradients at the integration point
+ const auto ipLocal = geometry.local(scvf.center());
+ jacInvT_ = geometry.jacobianInverseTransposed(ipLocal);
+ localBasis.evaluateJacobian(ipLocal, shapeJacobian_);
+ //localBasis.evaluateFunction(ipLocal, shapeValue_); // do we need the shapeValues for the flux?
+
+ // The stencil info is obsolete for the box method.
+ // It is here for compatibility with cc methods
+ stencil_ = Stencil(0);
+ }
+
+ const std::vector<ShapeJacobian>& shapeJacobian() const
+ { return shapeJacobian_; }
+
+ /* const std::vector<ShapeValue>& shapeValue() const
+ { return shapeValue_; }*/
+
+ const JacobianInverseTransposed& jacInvT() const
+ { return jacInvT_; }
+
+ const Stencil& stencil() const
+ {
+ return stencil_;
+ }
+
+private:
+ std::vector<ShapeJacobian> shapeJacobian_;
+ //std::vector<ShapeValue> shapeValue_;
+ JacobianInverseTransposed jacInvT_;
+
+ Stencil stencil_;
+};
+
+// specialization for the cell centered tpfa method
+template<class TypeTag>
+class FreeFlowFluxVariablesCacheImplementation<TypeTag, DiscretizationMethods::CCTpfa>
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using Stencil = std::vector<IndexType>;
+
+public:
+ void update(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace &scvf)
+ {
+ FluxVariables fluxVars;
+ stencil_ = fluxVars.computeStencil(problem, element, fvGeometry, scvf);
+ tij_ = AdvectionType::calculateTransmissibilities(problem, element, fvGeometry, elemVolVars, scvf);
+ }
+
+ const Stencil& stencil() const
+ { return stencil_; }
+
+ const Scalar& tij() const
+ { return tij_; }
+
+private:
+ Stencil stencil_;
+ Scalar tij_;
+};
+
+// specialization for the cell centered tpfa method
+template<class TypeTag>
+class FreeFlowFluxVariablesCacheImplementation<TypeTag, DiscretizationMethods::Staggered>
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using Stencil = std::vector<IndexType>;
+
+public:
+ void update(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace &scvf)
+ {
+ FluxVariables fluxVars;
+ stencil_ = fluxVars.computeCellCenterStencil(problem, element, fvGeometry, scvf);
+ tij_ = AdvectionType::calculateTransmissibilities(problem, element, fvGeometry, elemVolVars, scvf);
+ }
+
+ const Stencil& stencil() const
+ { return stencil_; }
+
+ const Scalar& tij() const
+ { return tij_; }
+
+private:
+ Stencil stencil_;
+ Scalar tij_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/staggered/indices.hh b/dumux/freeflow/staggered/indices.hh
new file mode 100644
index 0000000000..59bff37618
--- /dev/null
+++ b/dumux/freeflow/staggered/indices.hh
@@ -0,0 +1,44 @@
+// -*- 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 Defines the indices for the one-phase fully implicit model.
+ */
+#ifndef DUMUX_1P_INDICES_HH
+#define DUMUX_1P_INDICES_HH
+
+namespace Dumux
+{
+// \{
+
+/*!
+ * \ingroup OnePModel
+ * \ingroup ImplicitIndices
+ * \brief Indices for the one-phase model.
+ */
+struct OnePIndices
+{
+ static const int conti0EqIdx = 0; //index for the mass balance
+ static const int pressureIdx = 0; //index of the primary variable
+};
+
+// \}
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/staggered/model.hh b/dumux/freeflow/staggered/model.hh
new file mode 100644
index 0000000000..bb18ce217f
--- /dev/null
+++ b/dumux/freeflow/staggered/model.hh
@@ -0,0 +1,164 @@
+// -*- 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 Base class for all models which use the one-phase,
+ * fully implicit model.
+ * Adaption of the fully implicit scheme to the one-phase flow model.
+ */
+
+#ifndef DUMUX_1P_MODEL_HH
+#define DUMUX_1P_MODEL_HH
+
+#include <dumux/porousmediumflow/implicit/velocityoutput.hh>
+#include "properties.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup OnePModel
+ * \brief A single-phase, isothermal flow model using the fully implicit scheme.
+ *
+ * Single-phase, isothermal flow model, which uses a standard Darcy approach as the
+ * equation for the conservation of momentum:
+ * \f[
+ v = - \frac{\textbf K}{\mu}
+ \left(\textbf{grad}\, p - \varrho {\textbf g} \right)
+ * \f]
+ *
+ * and solves the mass continuity equation:
+ * \f[
+ \phi \frac{\partial \varrho}{\partial t} + \text{div} \left\lbrace
+ - \varrho \frac{\textbf K}{\mu} \left( \textbf{grad}\, p -\varrho {\textbf g} \right) \right\rbrace = q,
+ * \f]
+ * All equations are discretized using a vertex-centered finite volume (box)
+ * or cell-centered finite volume scheme as spatial
+ * and the implicit Euler method as time discretization.
+ * The model supports compressible as well as incompressible fluids.
+ */
+template<class TypeTag >
+class OnePModel : public GET_PROP_TYPE(TypeTag, BaseModel)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+ typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ enum { dim = GridView::dimension };
+ enum { dimWorld = GridView::dimensionworld };
+
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
+ using StencilsVector = typename GET_PROP_TYPE(TypeTag, StencilsVector);
+ using Element = typename GridView::template Codim<0>::Entity;
+
+public:
+ /*!
+ * \brief \copybrief Dumux::ImplicitModel::addOutputVtkFields
+ *
+ * Specialization for the OnePModel, adding the pressure and
+ * the process rank to the VTK writer.
+ */
+ template<class MultiWriter>
+ void addOutputVtkFields(const SolutionVector &sol,
+ MultiWriter &writer)
+ {
+ // typedef Dune::BlockVector<Dune::FieldVector<double, dimWorld> > VectorField;
+
+ // create the required scalar fields
+ unsigned numDofs = this->numDofs();
+ auto *p = writer.allocateManagedBuffer(numDofs);
+ auto *K = writer.allocateManagedBuffer(numDofs);
+ // VectorField *velocity = writer.template allocateManagedBuffer<double, dimWorld>(numDofs);
+ // ImplicitVelocityOutput<TypeTag> velocityOutput(this->problem_());
+
+ // if (velocityOutput.enableOutput())
+ // {
+ // // initialize velocity field
+ // for (unsigned int i = 0; i < numDofs; ++i)
+ // {
+ // (*velocity)[i] = double(0);
+ // }
+ // }
+
+ unsigned numElements = this->gridView_().size(0);
+ auto *rank = writer.allocateManagedBuffer(numElements);
+
+ for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
+ {
+ auto eIdx = this->elementMapper().index(element);
+ (*rank)[eIdx] = this->gridView_().comm().rank();
+
+ // get the local fv geometry
+ auto fvGeometry = localView(this->globalFvGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(this->curGlobalVolVars());
+ elemVolVars.bindElement(element, fvGeometry, this->curSol());
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+ const auto& spatialParams = this->problem_().spatialParams();
+ auto dofIdxGlobal = scv.dofIndex();
+
+ (*p)[dofIdxGlobal] = volVars.pressure();
+ (*K)[dofIdxGlobal] = spatialParams.intrinsicPermeability(scv, volVars);
+ }
+
+ // velocity output
+ //velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, element, /*phaseIdx=*/0);
+ }
+
+ writer.attachDofData(*p, "p", isBox);
+ writer.attachDofData(*K, "K", isBox);
+ // if (velocityOutput.enableOutput())
+ // {
+ // writer.attachDofData(*velocity, "velocity", isBox, dim);
+ // }
+ writer.attachCellData(*rank, "process rank");
+ }
+
+ /*!
+ * \brief Returns the number of global degrees of freedoms (DOFs)
+ */
+ size_t numDofs() const
+ {
+ return this->gridView_().size(0) + this->gridView_().size(1);;
+ }
+
+ size_t completeFaceDofStencilSize(const int idx)
+ {
+ return this->stencilsVector_.completeFaceDofStencilSize(idx);
+ }
+
+ template<class Stencil>
+ void getFullFaceDofStencils(std::unique_ptr<std::vector<Stencil>>& ptr)
+ {
+ this->stencilsVector_.getFullFaceDofStencils(ptr);
+ // TODO: find better way to do this
+ }
+
+};
+}
+
+#include "propertydefaults.hh"
+
+#endif
diff --git a/dumux/freeflow/staggered/properties.hh b/dumux/freeflow/staggered/properties.hh
new file mode 100644
index 0000000000..0d648ba013
--- /dev/null
+++ b/dumux/freeflow/staggered/properties.hh
@@ -0,0 +1,76 @@
+// -*- 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 Properties
+ * \ingroup ImplicitProperties
+ * \ingroup OnePModel
+ * \file
+ *
+ * \brief Defines the properties required for the one-phase fully implicit model.
+ */
+#ifndef DUMUX_1P_PROPERTIES_HH
+#define DUMUX_1P_PROPERTIES_HH
+
+#include <dumux/implicit/box/properties.hh>
+#include <dumux/implicit/cellcentered/properties.hh>
+#include <dumux/porousmediumflow/nonisothermal/implicit/properties.hh>
+
+namespace Dumux
+{
+// \{
+///////////////////////////////////////////////////////////////////////////
+// properties for the isothermal single phase model
+///////////////////////////////////////////////////////////////////////////
+namespace Properties {
+
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tags for the implicit single-phase problems
+NEW_TYPE_TAG(OneP);
+NEW_TYPE_TAG(BoxOneP, INHERITS_FROM(BoxModel, OneP));
+NEW_TYPE_TAG(CCOneP, INHERITS_FROM(CCModel, OneP));
+
+//! The type tags for the corresponding non-isothermal problems
+NEW_TYPE_TAG(OnePNI, INHERITS_FROM(OneP, NonIsothermal));
+NEW_TYPE_TAG(BoxOnePNI, INHERITS_FROM(BoxModel, OnePNI));
+NEW_TYPE_TAG(CCOnePNI, INHERITS_FROM(CCModel, OnePNI));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system
+NEW_PROP_TAG(Indices); //!< Enumerations for the model
+NEW_PROP_TAG(SpatialParams); //!< The type of the spatial parameters object
+NEW_PROP_TAG(FluidSystem); //!< The type of the fluid system to use
+NEW_PROP_TAG(Fluid); //!< The fluid used for the default fluid system
+NEW_PROP_TAG(FluidState); //!< The type of the fluid state to use
+NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem
+NEW_PROP_TAG(ImplicitMassUpwindWeight); //!< Returns weight of the upwind cell when calculating fluxes
+NEW_PROP_TAG(ImplicitMobilityUpwindWeight); //!< Weight for the upwind mobility in the velocity calculation
+NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient
+NEW_PROP_TAG(VtkAddVelocity); //!< Returns whether velocity vectors are written into the vtk output
+// \}
+}
+
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/staggered/propertydefaults.hh b/dumux/freeflow/staggered/propertydefaults.hh
new file mode 100644
index 0000000000..babb860120
--- /dev/null
+++ b/dumux/freeflow/staggered/propertydefaults.hh
@@ -0,0 +1,158 @@
+// -*- 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 Properties
+ * \ingroup ImplicitProperties
+ * \ingroup OnePModel
+ * \file
+ *
+ * \brief Defines the properties required for the one-phase fully implicit model.
+ */
+#ifndef DUMUX_1P_PROPERTY_DEFAULTS_HH
+#define DUMUX_1P_PROPERTY_DEFAULTS_HH
+
+#include "properties.hh"
+
+#include "model.hh"
+#include "volumevariables.hh"
+#include "indices.hh"
+
+#include <dumux/porousmediumflow/immiscible/localresidual.hh>
+#include <dumux/porousmediumflow/nonisothermal/implicit/propertydefaults.hh>
+#include <dumux/material/fluidsystems/gasphase.hh>
+#include <dumux/material/fluidsystems/liquidphase.hh>
+#include <dumux/material/components/nullcomponent.hh>
+#include <dumux/material/fluidsystems/1p.hh>
+#include <dumux/material/spatialparams/implicit1p.hh>
+#include <dumux/material/fluidmatrixinteractions/1p/thermalconductivityaverage.hh>
+
+namespace Dumux
+{
+// \{
+
+///////////////////////////////////////////////////////////////////////////
+// default property values for the isothermal single phase model
+///////////////////////////////////////////////////////////////////////////
+namespace Properties {
+SET_INT_PROP(OneP, NumEq, 1); //!< set the number of equations to 1
+SET_INT_PROP(OneP, NumPhases, 1); //!< The number of phases in the 1p model is 1
+
+//! The local residual function
+SET_TYPE_PROP(OneP, LocalResidual, ImmiscibleLocalResidual<TypeTag>);
+
+//! the Model property
+SET_TYPE_PROP(OneP, Model, OnePModel<TypeTag>);
+
+//! the VolumeVariables property
+SET_TYPE_PROP(OneP, VolumeVariables, OnePVolumeVariables<TypeTag>);
+
+//! Enable advection
+SET_BOOL_PROP(OneP, EnableAdvection, true);
+
+//! The one-phase model has no molecular diffusion
+SET_BOOL_PROP(OneP, EnableMolecularDiffusion, false);
+
+//! Isothermal model by default
+SET_BOOL_PROP(OneP, EnableEnergyBalance, false);
+
+//! The indices required by the isothermal single-phase model
+SET_TYPE_PROP(OneP, Indices, OnePIndices);
+
+//! The spatial parameters to be employed.
+//! Use ImplicitSpatialParamsOneP by default.
+SET_TYPE_PROP(OneP, SpatialParams, ImplicitSpatialParamsOneP<TypeTag>);
+
+//! The weight of the upwind control volume when calculating
+//! fluxes. Use central differences by default.
+SET_SCALAR_PROP(OneP, ImplicitMassUpwindWeight, 0.5);
+
+//! weight for the upwind mobility in the velocity calculation
+//! fluxes. Use central differences by default.
+SET_SCALAR_PROP(OneP, ImplicitMobilityUpwindWeight, 0.5);
+
+//! The fluid system to use by default
+SET_TYPE_PROP(OneP, FluidSystem, Dumux::FluidSystems::OneP<typename GET_PROP_TYPE(TypeTag, Scalar), typename GET_PROP_TYPE(TypeTag, Fluid)>);
+
+SET_PROP(OneP, Fluid)
+{ private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+public:
+ typedef Dumux::LiquidPhase<Scalar, Dumux::NullComponent<Scalar> > type;
+};
+
+/*!
+ * \brief The fluid state which is used by the volume variables to
+ * store the thermodynamic state. This should be chosen
+ * appropriately for the model ((non-)isothermal, equilibrium, ...).
+ * This can be done in the problem.
+ */
+SET_PROP(OneP, FluidState){
+ private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ public:
+ typedef Dumux::ImmiscibleFluidState<Scalar, FluidSystem> type;
+};
+
+// disable velocity output by default
+SET_BOOL_PROP(OneP, VtkAddVelocity, false);
+
+// enable gravity by default
+SET_BOOL_PROP(OneP, ProblemEnableGravity, true);
+
+//! default value for the forchheimer coefficient
+// Source: Ward, J.C. 1964 Turbulent flow in porous media. ASCE J. Hydraul. Div 90.
+// Actually the Forchheimer coefficient is also a function of the dimensions of the
+// porous medium. Taking it as a constant is only a first approximation
+// (Nield, Bejan, Convection in porous media, 2006, p. 10)
+SET_SCALAR_PROP(OneP, SpatialParamsForchCoeff, 0.55);
+
+//! average is used as default model to compute the effective thermal heat conductivity
+SET_PROP(OnePNI, ThermalConductivityModel)
+{ private :
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ public:
+ typedef ThermalConductivityAverage<Scalar> type;
+};
+
+//////////////////////////////////////////////////////////////////
+// Property values for isothermal model required for the general non-isothermal model
+//////////////////////////////////////////////////////////////////
+
+// set isothermal Model
+SET_TYPE_PROP(OnePNI, IsothermalModel, OnePModel<TypeTag>);
+
+//set isothermal VolumeVariables
+SET_TYPE_PROP(OnePNI, IsothermalVolumeVariables, OnePVolumeVariables<TypeTag>);
+
+//set isothermal LocalResidual
+SET_TYPE_PROP(OnePNI, IsothermalLocalResidual, ImmiscibleLocalResidual<TypeTag>);
+
+//set isothermal Indices
+SET_TYPE_PROP(OnePNI, IsothermalIndices, OnePIndices);
+
+//set isothermal NumEq
+SET_INT_PROP(OnePNI, IsothermalNumEq, 1);
+
+// \}
+} // end namespace Properties
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/freeflow/staggered/volumevariables.hh b/dumux/freeflow/staggered/volumevariables.hh
new file mode 100644
index 0000000000..81c65cb882
--- /dev/null
+++ b/dumux/freeflow/staggered/volumevariables.hh
@@ -0,0 +1,183 @@
+// -*- 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 Quantities required by the one-phase fully implicit model defined on a vertex.
+ */
+#ifndef DUMUX_1P_VOLUME_VARIABLES_HH
+#define DUMUX_1P_VOLUME_VARIABLES_HH
+
+#include "properties.hh"
+#include <dumux/discretization/volumevariables.hh>
+
+#include <dumux/material/fluidstates/immiscible.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup OnePModel
+ * \ingroup ImplicitVolumeVariables
+ * \brief Contains the quantities which are constant within a
+ * finite volume in the one-phase model.
+ */
+template <class TypeTag>
+class OnePVolumeVariables : public ImplicitVolumeVariables<TypeTag>
+{
+ using ParentType = ImplicitVolumeVariables<TypeTag>;
+ using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+
+public:
+
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+
+ /*!
+ * \copydoc ImplicitVolumeVariables::update
+ */
+ void update(const PrimaryVariables &priVars,
+ const Problem &problem,
+ const Element &element,
+ const SubControlVolume& scv)
+ {
+ ParentType::update(priVars, problem, element, scv);
+
+ completeFluidState(priVars, problem, element, scv, fluidState_);
+ // porosity
+ porosity_ = problem.spatialParams().porosity(scv);
+ };
+
+ /*!
+ * \copydoc ImplicitModel::completeFluidState
+ */
+ static void completeFluidState(const PrimaryVariables& priVars,
+ const Problem& problem,
+ const Element& element,
+ const SubControlVolume& scv,
+ FluidState& fluidState)
+ {
+ Scalar t = ParentType::temperature(priVars, problem, element, scv);
+ fluidState.setTemperature(t);
+ fluidState.setSaturation(/*phaseIdx=*/0, 1.);
+
+ fluidState.setPressure(/*phaseIdx=*/0, priVars[Indices::pressureIdx]);
+
+ // saturation in a single phase is always 1 and thus redundant
+ // to set. But since we use the fluid state shared by the
+ // immiscible multi-phase models, so we have to set it here...
+ fluidState.setSaturation(/*phaseIdx=*/0, 1.0);
+
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState, /*phaseIdx=*/0);
+
+ Scalar value = FluidSystem::density(fluidState, paramCache, /*phaseIdx=*/0);
+ fluidState.setDensity(/*phaseIdx=*/0, value);
+
+ value = FluidSystem::viscosity(fluidState, paramCache, /*phaseIdx=*/0);
+ fluidState.setViscosity(/*phaseIdx=*/0, value);
+
+ // compute and set the enthalpy
+ value = ParentType::enthalpy(fluidState, paramCache, /*phaseIdx=*/0);
+ fluidState.setEnthalpy(/*phaseIdx=*/0, value);
+ }
+
+ /*!
+ * \brief Return temperature \f$\mathrm{[K]}\f$ inside the sub-control volume.
+ *
+ * Note that we assume thermodynamic equilibrium, i.e. the
+ * temperatures of the rock matrix and of all fluid phases are
+ * identical.
+ */
+ Scalar temperature() const
+ { return fluidState_.temperature(); }
+
+ /*!
+ * \brief Return the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within
+ * the control volume.
+ */
+ Scalar pressure(int phaseIdx = 0) const
+ { return fluidState_.pressure(phaseIdx); }
+
+ /*!
+ * \brief Return the saturation
+ */
+ Scalar saturation(int phaseIdx = 0) const
+ { return 1.0; }
+
+ /*!
+ * \brief Return the mass density \f$\mathrm{[kg/m^3]}\f$ of a given phase within the
+ * control volume.
+ */
+ Scalar density(int phaseIdx = 0) const
+ { return fluidState_.density(phaseIdx); }
+
+ /*!
+ * \brief Return the dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
+ * control volume.
+ */
+ Scalar viscosity(int phaseIdx = 0) const
+ { return fluidState_.viscosity(phaseIdx); }
+
+ /*!
+ * \brief Returns the mobility \f$\mathrm{[1/(Pa s)]}\f$.
+ *
+ * This function enables the use of ImplicitDarcyFluxVariables
+ * with the 1p fully implicit model, ALTHOUGH the term mobility is
+ * usually not employed in the one phase context.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar mobility(int phaseIdx = 0) const
+ { return 1.0/fluidState_.viscosity(phaseIdx); }
+
+ /*!
+ * \brief Return the average porosity \f$\mathrm{[-]}\f$ within the control volume.
+ */
+ Scalar porosity() const
+ { return porosity_; }
+
+ /*!
+ * \brief Return the fluid state of the control volume.
+ */
+ const FluidState &fluidState() const
+ { return fluidState_; }
+
+protected:
+ FluidState fluidState_;
+ Scalar porosity_;
+
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+};
+
+}
+
+#endif
diff --git a/dumux/implicit/staggered/assembler.hh b/dumux/implicit/staggered/assembler.hh
new file mode 100644
index 0000000000..4d36f5390e
--- /dev/null
+++ b/dumux/implicit/staggered/assembler.hh
@@ -0,0 +1,106 @@
+// -*- 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 An assembler for the global Jacobian matrix for fully implicit models.
+ */
+#ifndef DUMUX_STAGGERED_ASSEMBLER_HH
+#define DUMUX_STAGGERED_ASSEMBLER_HH
+
+#include <dumux/implicit/properties.hh>
+#include <dumux/implicit/assembler.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief An assembler for the global Jacobian matrix for fully implicit models.
+ */
+template<class TypeTag>
+class StaggeredAssembler : public ImplicitAssembler<TypeTag>
+{
+ typedef ImplicitAssembler<TypeTag> ParentType;
+ friend class ImplicitAssembler<TypeTag>;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::IndexSet::IndexType IndexType;
+
+ void setRowSizes_()
+ {
+ for (const auto& element : elements(this->gridView_()))
+ {
+ // the global index of the element at hand
+ const auto globalI = this->elementMapper_().index(element);
+ const auto& stencil = this->model_().stencils(element).elementStencil();
+
+ this->matrix().setrowsize(globalI, stencil.size());
+ }
+
+ for(int facetIdx = 0; facetIdx < this->gridView_().size(1); ++facetIdx)
+ {
+ const int faceDofxIdx = facetIdx + this->gridView_().size(0);
+ const int size = this->model_().completeFaceDofStencilSize(facetIdx);
+ this->matrix().setrowsize(faceDofxIdx, size);
+ }
+ this->matrix().endrowsizes();
+ }
+
+ void addIndices_()
+ {
+ for (const auto& element : elements(this->gridView_()))
+ {
+ // the global index of the element at hand
+ const auto globalI = this->elementMapper_().index(element);
+ const auto& stencil = this->model_().stencils(element).elementStencil();
+
+
+ for (auto&& globalJ : stencil)
+ this->matrix().addindex(globalI, globalJ);
+ }
+
+ // TODO: this is a mess. find a better way for the pointer handling!
+ using Stencil = std::vector<IndexType>;
+ std::unique_ptr<std::vector<Stencil>> ptr;
+
+ this->model_().getFullFaceDofStencils(ptr);
+
+ if(ptr)
+ std::cout << "success!!!" << std::endl;
+
+ auto& fullFaceDofStencils = ptr.get()[0];
+
+ int globalI = this->gridView_().size(0);
+ for(const auto& stencil : fullFaceDofStencils)
+ {
+ for(auto&& globalJ : stencil)
+ this->matrix().addindex(globalI, globalJ);
+ ++globalI;
+ }
+
+ this->matrix().endindices();
+ }
+
+};
+
+} // namespace Dumux
+
+#endif
diff --git a/dumux/implicit/staggered/localjacobian.hh b/dumux/implicit/staggered/localjacobian.hh
new file mode 100644
index 0000000000..1791521d0e
--- /dev/null
+++ b/dumux/implicit/staggered/localjacobian.hh
@@ -0,0 +1,404 @@
+// -*- 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 Caculates the Jacobian of the local residual for fully-implicit models
+ */
+#ifndef DUMUX_STAGGERED_LOCAL_JACOBIAN_HH
+#define DUMUX_STAGGERED_LOCAL_JACOBIAN_HH
+
+#include <dune/istl/io.hh>
+#include <dune/istl/matrix.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/common/valgrind.hh>
+
+#include <dumux/implicit/properties.hh>
+#include <dumux/implicit/localjacobian.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup ImplicitLocalJacobian
+ * \brief Calculates the Jacobian of the local residual for fully-implicit models
+ *
+ * The default behavior is to use numeric differentiation, i.e.
+ * forward or backward differences (2nd order), or central
+ * differences (3rd order). The method used is determined by the
+ * "NumericDifferenceMethod" property:
+ *
+ * - if the value of this property is smaller than 0, backward
+ * differences are used, i.e.:
+ * \f[
+ \frac{\partial f(x)}{\partial x} \approx \frac{f(x) - f(x - \epsilon)}{\epsilon}
+ * \f]
+ *
+ * - if the value of this property is 0, central
+ * differences are used, i.e.:
+ * \f[
+ \frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x - \epsilon)}{2 \epsilon}
+ * \f]
+ *
+ * - if the value of this property is larger than 0, forward
+ * differences are used, i.e.:
+ * \f[
+ \frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x)}{\epsilon}
+ * \f]
+ *
+ * Here, \f$ f \f$ is the residual function for all equations, \f$x\f$
+ * is the value of a sub-control volume's primary variable at the
+ * evaluation point and \f$\epsilon\f$ is a small value larger than 0.
+ */
+template<class TypeTag>
+class StaggeredLocalJacobian : public ImplicitLocalJacobian<TypeTag>
+{
+ using ParentType = ImplicitLocalJacobian<TypeTag>;
+ using Implementation = typename GET_PROP_TYPE(TypeTag, LocalJacobian);
+ using LocalResidual = typename GET_PROP_TYPE(TypeTag, LocalResidual);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
+
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
+ using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using IndexType = typename GridView::IndexSet::IndexType;
+
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+
+ using AssemblyMap = std::vector<std::vector<std::vector<IndexType>>>;
+
+public:
+
+ StaggeredLocalJacobian() : ParentType()
+ {
+ numericDifferenceMethod_ = GET_PARAM_FROM_GROUP(TypeTag, int, Implicit, NumericDifferenceMethod);
+ }
+
+ /*!
+ * \brief Initialize the local Jacobian object.
+ *
+ * At this point we can assume that everything has been allocated,
+ * although some objects may not yet be completely initialized.
+ *
+ * \param problem The problem which we want to simulate.
+ */
+ void init(Problem &problem)
+ {
+ ParentType::init(problem);
+
+// assemblyMap_.resize(problem.gridView().size(0));
+// for (const auto& element : elements(problem.gridView()))
+// {
+// // get a local finite volume geometry object that is bindable
+// auto fvGeometryJ = localView(problem.model().globalFvGeometry());
+//
+// auto globalI = problem.elementMapper().index(element);
+// const auto& neighborStencil = this->model_().stencils(element).neighborStencil();
+//
+// // assemblyMap_[globalI].reserve(neighborStencil.size());
+// for (auto globalJ : neighborStencil)
+// {
+// const auto& elementJ = fvGeometryJ.globalFvGeometry().element(globalJ);
+//
+// // find the flux vars needed for the calculation of the flux into element
+// std::vector<IndexType> fluxVarIndices;
+//
+// // only non-ghost neighbors (J) have to be considered, derivatives from non-ghost to ghost dofs
+// // are assembled when assembling the ghost element (I)
+// if (elementJ.partitionType() != Dune::GhostEntity)
+// {
+// fvGeometryJ.bindElement(elementJ);
+// for (auto&& scvFaceJ : scvfs(fvGeometryJ))
+// {
+// auto fluxVarsIdx = scvFaceJ.index();
+//
+// // if globalI is in flux var stencil, add to list
+// FluxVariables fluxVars;
+// const auto fluxStencil = fluxVars.computeStencil(problem, elementJ, fvGeometryJ, scvFaceJ);
+//
+// for (auto globalIdx : fluxStencil)
+// if (globalIdx == globalI)
+// fluxVarIndices.push_back(fluxVarsIdx);
+// }
+// }
+// // assemblyMap_[globalI].emplace_back(std::move(fluxVarIndices));
+// }
+// }
+ }
+
+ /*!
+ * \brief Assemble an element's local Jacobian matrix of the
+ * defect.
+ *
+ * \param element The DUNE Codim<0> entity which we look at.
+ */
+ void assemble(const Element& element, JacobianMatrix& matrix, SolutionVector& residual)
+ {
+ const bool isGhost = (element.partitionType() == Dune::GhostEntity);
+
+ // prepare the volvars/fvGeometries in case caching is disabled
+ auto fvGeometry = localView(this->model_().globalFvGeometry());
+ fvGeometry.bind(element);
+
+ auto curElemVolVars = localView(this->model_().curGlobalVolVars());
+ curElemVolVars.bind(element, fvGeometry, this->model_().curSol());
+
+ auto prevElemVolVars = localView(this->model_().prevGlobalVolVars());
+ prevElemVolVars.bindElement(element, fvGeometry, this->model_().prevSol());
+
+ auto elemFluxVarsCache = localView(this->model_().globalFluxVarsCache());
+ elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
+
+ // set the actual dof index
+ globalI_ = this->problem_().elementMapper().index(element);
+
+ ElementBoundaryTypes elemBcTypes;
+ elemBcTypes.update(this->problem_(), element, fvGeometry);
+
+ // calculate the local residual
+ if (isGhost)
+ {
+ this->residual_ = 0.0;
+ residual[globalI_] = 0.0;
+ }
+ else
+ {
+ this->localResidual().eval(element, fvGeometry, prevElemVolVars, curElemVolVars, elemBcTypes, elemFluxVarsCache);
+ this->residual_ = this->localResidual().residual();
+ // store residual in global container as well
+ residual[globalI_] = this->localResidual().residual(0);
+ }
+
+ this->model_().updatePVWeights(fvGeometry);
+
+ // calculate derivatives of all dofs in stencil with respect to the dofs in the element
+ evalPartialDerivatives_(element, fvGeometry, prevElemVolVars, curElemVolVars, elemFluxVarsCache, elemBcTypes, matrix, residual, isGhost);
+
+ // TODO: calculate derivatives in the case of an extended source stencil
+ // const auto& extendedSourceStencil = model_().stencils(element).extendedSourceStencil();
+ // for (auto&& globalJ : extendedSourceStencil)
+ // {
+ // for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
+ // {
+ // evalPartialDerivativeSource_(partialDeriv, globalJ, pvIdx, neighborToFluxVars[globalJ]);
+
+ // // update the local stiffness matrix with the partial derivatives
+ // updateLocalJacobian_(j, pvIdx, partialDeriv);
+ // }
+ // ++j;
+ // }
+ }
+
+ const AssemblyMap& assemblyMap() const
+ { return assemblyMap_; }
+
+private:
+ void evalPartialDerivatives_(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& prevElemVolVars,
+ ElementVolumeVariables& curElemVolVars,
+ ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes& elemBcTypes,
+ JacobianMatrix& matrix,
+ SolutionVector& residual,
+ const bool isGhost)
+ {
+ // get stencil informations
+ const auto& neighborStencil = this->model_().stencils(element).neighborStencil();
+ const auto numNeighbors = neighborStencil.size();
+
+ // the localresidual class used for the flux calculations
+ LocalResidual localRes;
+ localRes.init(this->problem_());
+
+ // container to store the neighboring elements
+ std::vector<Element> neighborElements;
+ neighborElements.reserve(numNeighbors);
+
+ // get the elements and calculate the flux into the element in the undeflected state
+ Dune::BlockVector<PrimaryVariables> origFlux(numNeighbors);
+ origFlux = 0.0;
+
+ unsigned int j = 0;
+ for (auto globalJ : neighborStencil)
+ {
+ auto elementJ = fvGeometry.globalFvGeometry().element(globalJ);
+ neighborElements.push_back(elementJ);
+
+ for (auto fluxVarIdx : assemblyMap_[globalI_][j])
+ {
+ auto&& scvf = fvGeometry.scvf(fluxVarIdx);
+ origFlux[j] += localRes.evalFlux_(elementJ, fvGeometry, curElemVolVars, scvf, elemFluxVarsCache[scvf]);
+ }
+
+ ++j;
+ }
+
+ auto&& scv = fvGeometry.scv(globalI_);
+ auto& curVolVars = getCurVolVars(curElemVolVars, scv);
+ // save a copy of the original vol vars
+ VolumeVariables origVolVars(curVolVars);
+
+ // derivatives in the neighbors with repect to the current elements
+ Dune::BlockVector<PrimaryVariables> neighborDeriv(numNeighbors);
+ for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
+ {
+ // derivatives of element dof with respect to itself
+ PrimaryVariables partialDeriv(0.0);
+
+ if (isGhost)
+ partialDeriv[pvIdx] = 1.0;
+
+ neighborDeriv = 0.0;
+ PrimaryVariables priVars(this->model_().curSol()[globalI_]);
+
+ Scalar eps = this->numericEpsilon(scv, curVolVars, pvIdx);
+ Scalar delta = 0;
+
+ if (numericDifferenceMethod_ >= 0)
+ {
+ // we are not using backward differences, i.e. we need to
+ // calculate f(x + \epsilon)
+
+ // deflect primary variables
+ priVars[pvIdx] += eps;
+ delta += eps;
+
+ // update the volume variables and bind the flux var cache again
+ curVolVars.update(priVars, this->problem_(), element, scv);
+ elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
+
+ if (!isGhost)
+ {
+ // calculate the residual with the deflected primary variables
+ this->localResidual().eval(element, fvGeometry, prevElemVolVars, curElemVolVars, elemBcTypes, elemFluxVarsCache);
+
+ // store the residual and the storage term
+ partialDeriv = this->localResidual().residual(0);
+ }
+
+ // calculate the fluxes into element with the deflected primary variables
+ for (std::size_t k = 0; k < numNeighbors; ++k)
+ {
+ for (auto fluxVarIdx : assemblyMap_[globalI_][k])
+ {
+ auto&& scvf = fvGeometry.scvf(fluxVarIdx);
+ neighborDeriv[k] += localRes.evalFlux_(neighborElements[k], fvGeometry, curElemVolVars, scvf, elemFluxVarsCache[scvf]);
+ }
+ }
+ }
+ else
+ {
+ // we are using backward differences, i.e. we don't need
+ // to calculate f(x + \epsilon) and we can recycle the
+ // (already calculated) residual f(x)
+ if (!isGhost)
+ partialDeriv = this->residual(0);
+ neighborDeriv = origFlux;
+ }
+
+ if (numericDifferenceMethod_ <= 0)
+ {
+ // we are not using forward differences, i.e. we
+ // need to calculate f(x - \epsilon)
+
+ // deflect the primary variables
+ priVars[pvIdx] -= delta + eps;
+ delta += eps;
+
+ // update the volume variables and bind the flux var cache again
+ curVolVars.update(priVars, this->problem_(), element, scv);
+ elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
+
+ if (!isGhost)
+ {
+ // calculate the residual with the deflected primary variables
+ this->localResidual().eval(element, fvGeometry, prevElemVolVars, curElemVolVars, elemBcTypes, elemFluxVarsCache);
+
+ // subtract the residual from the derivative storage
+ partialDeriv -= this->localResidual().residual(0);
+ }
+
+ // calculate the fluxes into element with the deflected primary variables
+ for (std::size_t k = 0; k < numNeighbors; ++k)
+ {
+ for (auto fluxVarIdx : assemblyMap_[globalI_][k])
+ {
+ auto&& scvf = fvGeometry.scvf(fluxVarIdx);
+ neighborDeriv[k] -= localRes.evalFlux_(neighborElements[k], fvGeometry, curElemVolVars, scvf, elemFluxVarsCache[scvf]);
+ }
+ }
+ }
+ else
+ {
+ // we are using forward differences, i.e. we don't need to
+ // calculate f(x - \epsilon) and we can recycle the
+ // (already calculated) residual f(x)
+ if (!isGhost)
+ partialDeriv -= this->residual(0);
+ neighborDeriv -= origFlux;
+ }
+
+ // divide difference in residuals by the magnitude of the
+ // deflections between the two function evaluation
+ if (!isGhost)
+ partialDeriv /= delta;
+ neighborDeriv /= delta;
+
+ // restore the original state of the scv's volume variables
+ curVolVars = origVolVars;
+
+ // update the global jacobian matrix with the current partial derivatives
+ this->updateGlobalJacobian_(matrix, globalI_, globalI_, pvIdx, partialDeriv);
+
+ j = 0;
+ for (auto globalJ : neighborStencil)
+ this->updateGlobalJacobian_(matrix, globalJ, globalI_, pvIdx, neighborDeriv[j++]);
+ }
+ }
+
+ //! If the global vol vars caching is enabled we have to modify the global volvar object
+ template<class T = TypeTag>
+ typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), VolumeVariables>::type&
+ getCurVolVars(ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
+ { return this->model_().nonConstCurGlobalVolVars().volVars(scv.index()); }
+
+ //! When global volume variables caching is disabled, return the local volvar object
+ template<class T = TypeTag>
+ typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), VolumeVariables>::type&
+ getCurVolVars(ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
+ { return elemVolVars[scv]; }
+
+ IndexType globalI_;
+ int numericDifferenceMethod_;
+ AssemblyMap assemblyMap_;
+};
+
+}
+
+#endif
diff --git a/dumux/implicit/staggered/localresidual.hh b/dumux/implicit/staggered/localresidual.hh
new file mode 100644
index 0000000000..84e1fd1ff8
--- /dev/null
+++ b/dumux/implicit/staggered/localresidual.hh
@@ -0,0 +1,306 @@
+// -*- 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_LOCAL_RESIDUAL_HH
+#define DUMUX_STAGGERED_LOCAL_RESIDUAL_HH
+
+#include <dune/istl/matrix.hh>
+
+#include <dumux/common/valgrind.hh>
+#include <dumux/implicit/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 StaggeredLocalResidual : public ImplicitLocalResidual<TypeTag>
+{
+ using ParentType = ImplicitLocalResidual<TypeTag>;
+ friend class ImplicitLocalResidual<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+
+ 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);
+
+public:
+ // copying the local residual class is not a good idea
+ StaggeredLocalResidual(const StaggeredLocalResidual &) = delete;
+
+ StaggeredLocalResidual() : ParentType() {}
+
+protected:
+
+ void evalFluxes_(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ 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]);
+ }
+ }
+
+ PrimaryVariables computeFlux_(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace &scvf,
+ const ElementBoundaryTypes& bcTypes,
+ const FluxVariablesCache& fluxVarsCache)
+ {
+ if (!scvf.boundary() /*TODO: || GET_PROP_VALUE(TypeTag, BoundaryReconstruction)*/)
+ return this->asImp_().computeFlux(element, fvGeometry, elemVolVars, scvf, fluxVarsCache);
+ else
+ return PrimaryVariables(0.0);
+
+ }
+
+ void evalBoundary_(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementBoundaryTypes& elemBcTypes,
+ const ElementFluxVariablesCache& elemFluxVarsCache)
+ {
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ if (scvf.boundary())
+ {
+ auto bcTypes = this->problem().boundaryTypes(element, scvf);
+ this->residual_[0] += evalBoundaryFluxes_(element, fvGeometry, elemVolVars, scvf, bcTypes, elemFluxVarsCache[scvf]);
+ }
+ }
+
+ // additionally treat mixed D/N conditions in a strong sense
+ if (elemBcTypes.hasDirichlet())
+ {
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ if (scvf.boundary())
+ this->asImp_().evalDirichlet_(element, fvGeometry, elemVolVars, scvf);
+ }
+ }
+ }
+
+ /*!
+ * \brief Add all fluxes resulting from Neumann, outflow and pure Dirichlet
+ * boundary conditions to the local residual.
+ */
+ PrimaryVariables evalBoundaryFluxes_(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace &scvf,
+ const BoundaryTypes& bcTypes,
+ const FluxVariablesCache& fluxVarsCache)
+ {
+ // evaluate the Neumann conditions at the boundary face
+ PrimaryVariables flux(0);
+ if (bcTypes.hasNeumann() /*TODO: && !GET_PROP_VALUE(TypeTag, BoundaryReconstruction)*/)
+ flux += this->asImp_().evalNeumannSegment_(element, fvGeometry, elemVolVars, scvf, bcTypes);
+
+ // TODO: evaluate the outflow conditions at the boundary face
+ //if (bcTypes.hasOutflow() /*TODO: && !GET_PROP_VALUE(TypeTag, BoundaryReconstruction)*/)
+ // flux += this->asImp_().evalOutflowSegment_(&intersection, bcTypes);
+
+ // evaluate the pure Dirichlet conditions at the boundary face
+ if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann())
+ flux += this->asImp_().evalDirichletSegment_(element, fvGeometry, elemVolVars, scvf, bcTypes, fluxVarsCache);
+
+ return flux;
+ }
+
+
+ /*!
+ * \brief Evaluate Dirichlet conditions on faces that have mixed
+ * Dirichlet/Neumann boundary conditions.
+ */
+ void evalDirichlet_(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace &scvf)
+ {
+ BoundaryTypes bcTypes = this->problem().boundaryTypes(element, scvf);
+
+ if (bcTypes.hasDirichlet() && bcTypes.hasNeumann())
+ this->asImp_().evalDirichletSegmentMixed_(element, fvGeometry, elemVolVars, scvf, bcTypes);
+ }
+
+ /*!
+ * \brief Add Neumann boundary conditions for a single scv face
+ */
+ PrimaryVariables evalNeumannSegment_(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace &scvf,
+ const BoundaryTypes &bcTypes)
+ {
+ // temporary vector to store the neumann boundary fluxes
+ PrimaryVariables flux(0);
+
+ auto neumannFluxes = this->problem().neumann(element, fvGeometry, elemVolVars, scvf);
+
+ // multiply neumann fluxes with the area and the extrusion factor
+ auto&& scv = fvGeometry.scv(scvf.insideScvIdx());
+ neumannFluxes *= scvf.area()*elemVolVars[scv].extrusionFactor();
+
+ // add fluxes to the residual
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ if (bcTypes.isNeumann(eqIdx))
+ flux[eqIdx] += neumannFluxes[eqIdx];
+
+ return flux;
+ }
+
+ /*!
+ * \brief Treat Dirichlet boundary conditions in a weak sense for a single
+ * intersection that only has Dirichlet boundary conditions
+ */
+ PrimaryVariables evalDirichletSegment_(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace &scvf,
+ const BoundaryTypes &bcTypes,
+ const FluxVariablesCache& fluxVarsCache)
+ {
+ // temporary vector to store the Dirichlet boundary fluxes
+ PrimaryVariables flux(0);
+
+ auto dirichletFlux = this->asImp_().computeFlux(element, fvGeometry, elemVolVars, scvf, fluxVarsCache);
+
+ // add fluxes to the residual
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ if (bcTypes.isDirichlet(eqIdx))
+ flux[eqIdx] += dirichletFlux[eqIdx];
+
+ return flux;
+ }
+
+ /*!
+ * \brief Treat Dirichlet boundary conditions in a strong sense for a
+ * single intersection that has mixed D/N boundary conditions
+ */
+ void evalDirichletSegmentMixed_(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace &scvf,
+ 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];
+ }
+ }
+ }
+
+ /*!
+ * \brief Add outflow boundary conditions for a single intersection
+ */
+ /*template <class IntersectionIterator>
+ void evalOutflowSegment_(const IntersectionIterator &isIt,
+ const BoundaryTypes &bcTypes)
+ {
+ if (this->element_().geometry().type().isCube() == false)
+ DUNE_THROW(Dune::InvalidStateException,
+ "for cell-centered models, outflow BCs only work for cubes.");
+
+ // store pointer to the current FVElementGeometry
+ const FVElementGeometry *oldFVGeometryPtr = this->fvElemGeomPtr_;
+
+ // copy the current FVElementGeometry to a local variable
+ // and set the pointer to this local variable
+ FVElementGeometry fvGeometry = this->fvGeometry_();
+ this->fvElemGeomPtr_ = &fvGeometry;
+
+ // get the index of the boundary face
+ unsigned bfIdx = isIt->indexInInside();
+ unsigned oppositeIdx = bfIdx^1;
+
+ // manipulate the corresponding subcontrolvolume face
+ SCVFace& boundaryFace = fvGeometry.boundaryFace[bfIdx];
+
+ // set the second flux approximation index for the boundary face
+ for (int nIdx = 0; nIdx < fvGeometry.numNeighbors-1; nIdx++)
+ {
+ // check whether the two faces are opposite of each other
+ if (fvGeometry.subContVolFace[nIdx].fIdx == oppositeIdx)
+ {
+ boundaryFace.j = nIdx+1;
+ break;
+ }
+ }
+ boundaryFace.fapIndices[1] = boundaryFace.j;
+ boundaryFace.grad[0] *= -0.5;
+ boundaryFace.grad[1] *= -0.5;
+
+ // temporary vector to store the outflow boundary fluxes
+ PrimaryVariables values;
+ Valgrind::SetUndefined(values);
+
+ this->asImp_().computeFlux(values, bfIdx, true);
+ values *= this->curVolVars_(0).extrusionFactor();
+
+ // add fluxes to the residual
+ Valgrind::CheckDefined(values);
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (bcTypes.isOutflow(eqIdx))
+ this->residual_[0][eqIdx] += values[eqIdx];
+ }
+
+ // restore the pointer to the FVElementGeometry
+ this->fvElemGeomPtr_ = oldFVGeometryPtr;
+ }*/
+};
+
+}
+
+#endif // DUMUX_CC_LOCAL_RESIDUAL_HH
diff --git a/dumux/implicit/staggered/propertydefaults.hh b/dumux/implicit/staggered/propertydefaults.hh
index 6f49cbfaba..1b385410ef 100644
--- a/dumux/implicit/staggered/propertydefaults.hh
+++ b/dumux/implicit/staggered/propertydefaults.hh
@@ -28,13 +28,24 @@
#define DUMUX_STAGGERED_PROPERTY_DEFAULTS_HH
#include <dumux/implicit/propertydefaults.hh>
-#include <dumux/porousmediumflow/implicit/fluxvariablescache.hh>
+// #include <dumux/porousmediumflow/implicit/fluxvariablescache.hh>
#include <dumux/discretization/staggered/globalfvgeometry.hh>
#include <dumux/discretization/staggered/fvelementgeometry.hh>
#include <dumux/discretization/staggered/subcontrolvolumeface.hh>
#include <dumux/implicit/staggered/properties.hh>
#include <dumux/discretization/methods.hh>
-#include <dumux/discretization/cellcentered/stencils.hh>
+#include <dumux/discretization/staggered/stencils.hh>
+
+
+#include <dumux/freeflow/staggered/fluxvariables.hh>
+#include <dumux/freeflow/staggered/fluxvariablescache.hh>
+
+
+
+
+#include "assembler.hh"
+#include "localresidual.hh"
+#include "localjacobian.hh"
namespace Dumux {
@@ -88,13 +99,13 @@ SET_TYPE_PROP(StaggeredModel, GlobalVolumeVariables, Dumux::CCGlobalVolumeVariab
SET_TYPE_PROP(StaggeredModel, GlobalFluxVariablesCache, Dumux::CCGlobalFluxVariablesCache<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFluxVariablesCache)>);
//! The local jacobian operator
-SET_TYPE_PROP(StaggeredModel, LocalJacobian, Dumux::CCLocalJacobian<TypeTag>);
+SET_TYPE_PROP(StaggeredModel, LocalJacobian, Dumux::StaggeredLocalJacobian<TypeTag>);
//! Assembler for the global jacobian matrix
-SET_TYPE_PROP(StaggeredModel, JacobianAssembler, Dumux::CCAssembler<TypeTag>);
+SET_TYPE_PROP(StaggeredModel, JacobianAssembler, Dumux::StaggeredAssembler<TypeTag>);
//! The stencil container
-SET_TYPE_PROP(StaggeredModel, StencilsVector, Dumux::CCStencilsVector<TypeTag>);
+SET_TYPE_PROP(StaggeredModel, StencilsVector, Dumux::StaggeredStencilsVector<TypeTag>);
//! The local flux variables cache vector class
SET_TYPE_PROP(StaggeredModel, ElementFluxVariablesCache, Dumux::CCElementFluxVariablesCache<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFluxVariablesCache)>);
@@ -102,12 +113,19 @@ SET_TYPE_PROP(StaggeredModel, ElementFluxVariablesCache, Dumux::CCElementFluxVar
//! The global previous volume variables vector class
SET_TYPE_PROP(StaggeredModel, ElementVolumeVariables, Dumux::CCElementVolumeVariables<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalVolumeVariablesCache)>);
-//! Set the BaseLocalResidual to CCLocalResidual
-SET_TYPE_PROP(StaggeredModel, BaseLocalResidual, Dumux::CCLocalResidual<TypeTag>);
+//! Set the BaseLocalResidual to StaggeredLocalResidual
+SET_TYPE_PROP(StaggeredModel, BaseLocalResidual, Dumux::StaggeredLocalResidual<TypeTag>);
//! indicate that this is no box discretization
SET_BOOL_PROP(StaggeredModel, ImplicitIsBox, false);
+//! The class that contains the different flux variables (i.e. darcy, diffusion, energy)
+//! by default, we set the flux variables to ones for porous media
+SET_TYPE_PROP(StaggeredModel, FluxVariables, FreeFlowFluxVariables<TypeTag>);
+
+//! The flux variables cache class, by default the one for porous media
+SET_TYPE_PROP(StaggeredModel, FluxVariablesCache, FreeFlowFluxVariablesCache<TypeTag>);
+
} // namespace Properties
} // namespace Dumux
diff --git a/test/freeflow/CMakeLists.txt b/test/freeflow/CMakeLists.txt
index 358d7dd60a..f9eb646820 100644
--- a/test/freeflow/CMakeLists.txt
+++ b/test/freeflow/CMakeLists.txt
@@ -1,4 +1,5 @@
add_subdirectory("navierstokes")
+add_subdirectory("staggered")
add_subdirectory("stokes")
add_subdirectory("stokes2c")
add_subdirectory("stokes2cni")
--
GitLab