diff --git a/dumux/assembly/staggeredfvassembler.hh b/dumux/assembly/staggeredfvassembler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..da1834dd3ef8a4df2e42a8f3e30a3fbc6958f41a
--- /dev/null
+++ b/dumux/assembly/staggeredfvassembler.hh
@@ -0,0 +1,457 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief An assembler for the global linear system
+ * for fully implicit models and vertex-centered discretization schemes.
+ */
+#ifndef DUMUX_STAGGERED_FV_ASSEMBLER_HH
+#define DUMUX_STAGGERED_FV_ASSEMBLER_HH
+
+#include
+
+#include
+
+#include
+#include
+#include
+
+#include "diffmethod.hh"
+#include "staggeredlocalassembler.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief An assembler for the global linear system
+ * for fully implicit models and cell-centered discretization schemes.
+ */
+template
+class StaggeredFVAssembler
+{
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using LocalResidual = typename GET_PROP_TYPE(TypeTag, LocalResidual);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using TimeLoop = TimeLoopBase;
+
+ static constexpr int dim = GridView::dimension;
+ static constexpr bool isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox);
+ static constexpr int dofCodim = isBox ? dim : 0;
+
+ using LocalAssembler =StaggeredLocalAssembler;
+
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+ using CCToCCMatrixBlock = typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockCCToCC;
+ using CCToFaceMatrixBlock = typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockCCToFace;
+ using FaceToFaceMatrixBlock = typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockFaceToFace;
+ using FaceToCCMatrixBlock = typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockFaceToCC;
+
+public:
+ using ResidualType = SolutionVector;
+
+ //! The constructor for stationary problems
+ StaggeredFVAssembler(std::shared_ptr problem,
+ std::shared_ptr fvGridGeometry,
+ std::shared_ptr gridVariables)
+ : problem_(problem)
+ , fvGridGeometry_(fvGridGeometry)
+ , gridVariables_(gridVariables)
+ , stationary_(true)
+ {
+ static_assert(isImplicit, "Explicit assembler for stationary problem doesn't make sense!");
+ }
+
+ //! The constructor for instationary problems
+ StaggeredFVAssembler(std::shared_ptr problem,
+ std::shared_ptr fvGridGeometry,
+ std::shared_ptr gridVariables,
+ std::shared_ptr timeLoop)
+ : problem_(problem)
+ , fvGridGeometry_(fvGridGeometry)
+ , gridVariables_(gridVariables)
+ , localResidual_(timeLoop)
+ , stationary_(false)
+ {}
+
+ /*!
+ * \brief Assembles the global Jacobian of the residual
+ * and the residual for the current solution.
+ */
+ void assembleJacobianAndResidual(const SolutionVector& curSol)
+ {
+ if (!stationary_ && localResidual_.isStationary())
+ DUNE_THROW(Dune::InvalidStateException, "Assembling instationary problem but previous solution was not set!");
+
+ if(!jacobian_)
+ {
+ setLinearSystem();
+ }
+
+ if(!residual_)
+ {
+ residual_ = std::make_shared();
+ setResidualSize();
+ }
+
+ resetJacobian_();
+ resetResidual_();
+
+ bool succeeded;
+ // try assembling the global linear system
+ try
+ {
+ // let the local assembler add the element contributions
+ for (const auto element : elements(gridView()))
+ LocalAssembler::assembleJacobianAndResidual(*this, *jacobian_, *residual_, element, curSol);
+
+ // if we get here, everything worked well
+ succeeded = true;
+ if (gridView().comm().size() > 1)
+ succeeded = gridView().comm().min(succeeded);
+
+ // printmatrix(std::cout, (*jacobian_)[cellCenterIdx][cellCenterIdx], "A11", "");
+ // printmatrix(std::cout, (*jacobian_)[cellCenterIdx][faceIdx], "A12", "");
+ // printmatrix(std::cout, (*jacobian_)[faceIdx][cellCenterIdx], "A21", "");
+ // printmatrix(std::cout, (*jacobian_)[faceIdx][faceIdx], "A22", "");
+ }
+ // throw exception if a problem ocurred
+ catch (NumericalProblem &e)
+ {
+ std::cout << "rank " << gridView().comm().rank()
+ << " caught an exception while assembling:" << e.what()
+ << "\n";
+ succeeded = false;
+ if (gridView().comm().size() > 1)
+ succeeded = gridView().comm().min(succeeded);
+ }
+ if (!succeeded)
+ DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system");
+ }
+
+ /*!
+ * \brief Assembles only the global Jacobian of the residual.
+ */
+ void assembleJacobian(const SolutionVector& curSol)
+ {
+ if (!stationary_ && localResidual_.isStationary())
+ DUNE_THROW(Dune::InvalidStateException, "Assembling instationary problem but previous solution was not set!");
+
+ if(!jacobian_)
+ {
+ jacobian_ = std::make_shared();
+ jacobian_->setBuildMode(JacobianMatrix::random);
+ setJacobianPattern();
+ }
+
+ resetJacobian_();
+
+ bool succeeded;
+ // try assembling the global linear system
+ try
+ {
+ // let the local assembler add the element contributions
+ for (const auto element : elements(gridView()))
+ LocalAssembler::assemble(*this, *jacobian_, element, curSol);
+
+ // if we get here, everything worked well
+ succeeded = true;
+ if (gridView().comm().size() > 1)
+ succeeded = gridView().comm().min(succeeded);
+ }
+ // throw exception if a problem ocurred
+ catch (NumericalProblem &e)
+ {
+ std::cout << "rank " << gridView().comm().rank()
+ << " caught an exception while assembling:" << e.what()
+ << "\n";
+ succeeded = false;
+ if (gridView().comm().size() > 1)
+ succeeded = gridView().comm().min(succeeded);
+ }
+ if (!succeeded)
+ DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system");
+ }
+
+ //! compute the residuals
+ void assembleResidual(const SolutionVector& curSol)
+ {
+ if(!residual_)
+ {
+ residual_ = std::make_shared();
+ setResidualSize();
+ }
+
+ assembleResidual(*residual_, curSol);
+ }
+
+ //! compute the residuals
+ void assembleResidual(ResidualType& r, const SolutionVector& curSol) const
+ {
+ if (!stationary_ && localResidual_.isStationary())
+ DUNE_THROW(Dune::InvalidStateException, "Assembling instationary problem but previous solution was not set!");
+
+ // let the local assembler add the element contributions
+ for (const auto element : elements(gridView()))
+ LocalAssembler::assemble(*this, r, element, curSol);
+ }
+
+ //! computes the residual norm
+ Scalar residualNorm(const SolutionVector& curSol) const
+ {
+ ResidualType residual;
+ residual[cellCenterIdx].resize(fvGridGeometry().numCellCenterDofs());
+ residual[faceIdx].resize(fvGridGeometry().numFaceDofs());
+ assembleResidual(residual, curSol);
+
+ // calculate the square norm of the residual
+ Scalar result2 = residual.two_norm2();
+ if (gridView().comm().size() > 1)
+ result2 = gridView().comm().sum(result2);
+
+ using std::sqrt;
+ return sqrt(result2);
+ }
+
+ /*!
+ * \brief Tells the assembler which jacobian and residual to use.
+ * This also resizes the containers to the required sizes and sets the
+ * sparsity pattern of the jacobian matrix.
+ */
+ void setLinearSystem(std::shared_ptr A,
+ std::shared_ptr r)
+ {
+ jacobian_ = A;
+ residual_ = r;
+ //
+ // check and/or set the BCRS matrix's build mode
+ // convenience references
+ CCToCCMatrixBlock& A11 = (*jacobian_)[cellCenterIdx][cellCenterIdx];
+ CCToFaceMatrixBlock& A12 = (*jacobian_)[cellCenterIdx][faceIdx];
+ FaceToCCMatrixBlock& A21 = (*jacobian_)[faceIdx][cellCenterIdx];
+ FaceToFaceMatrixBlock& A22 = (*jacobian_)[faceIdx][faceIdx];
+
+ A11.setBuildMode(CCToCCMatrixBlock::random);
+ A12.setBuildMode(CCToFaceMatrixBlock::random);
+ A21.setBuildMode(FaceToCCMatrixBlock::random);
+ A22.setBuildMode(FaceToFaceMatrixBlock::random);
+
+ setJacobianPattern();
+ setResidualSize();
+ // if (jacobian_->buildMode() == JacobianMatrix::BuildMode::unknown)
+ // jacobian_->setBuildMode(JacobianMatrix::random);
+ // else if (jacobian_->buildMode() != JacobianMatrix::BuildMode::random)
+ // DUNE_THROW(Dune::NotImplemented, "Only BCRS matrices with random build mode are supported at the moment");
+ //
+ // setJacobianPattern();
+ // setResidualSize();
+ }
+
+ /*!
+ * \brief The version without arguments uses the default constructor to create
+ * the jacobian and residual objects in this assembler.
+ */
+ void setLinearSystem()
+ {
+ jacobian_ = std::make_shared();
+
+ // convenience references
+ CCToCCMatrixBlock& A11 = (*jacobian_)[cellCenterIdx][cellCenterIdx];
+ CCToFaceMatrixBlock& A12 = (*jacobian_)[cellCenterIdx][faceIdx];
+ FaceToCCMatrixBlock& A21 = (*jacobian_)[faceIdx][cellCenterIdx];
+ FaceToFaceMatrixBlock& A22 = (*jacobian_)[faceIdx][faceIdx];
+
+ A11.setBuildMode(CCToCCMatrixBlock::random);
+ A12.setBuildMode(CCToFaceMatrixBlock::random);
+ A21.setBuildMode(FaceToCCMatrixBlock::random);
+ A22.setBuildMode(FaceToFaceMatrixBlock::random);
+
+ residual_ = std::make_shared();
+
+ setJacobianPattern();
+ setResidualSize();
+ }
+
+ /*!
+ * \brief Sets the solution from which to start the time integration. Has to be
+ * called prior to assembly for time-dependent problems.
+ */
+ void setPreviousSolution(const SolutionVector& u)
+ { localResidual_.setPreviousSolution(u); }
+
+ /*!
+ * \brief Return the solution that has been set as the previous one.
+ */
+ const SolutionVector& prevSol() const
+ { return localResidual_.prevSol(); }
+
+ /*!
+ * \brief Resizes the jacobian and sets the jacobian' sparsity pattern.
+ */
+
+ void setJacobianPattern()
+ {
+ // resize the jacobian and the residual
+ const auto numDofsCC = fvGridGeometry().numCellCenterDofs();
+ const auto numDofsFace = fvGridGeometry().numFaceDofs();
+
+ // convenience references
+ CCToCCMatrixBlock& A11 = (*jacobian_)[cellCenterIdx][cellCenterIdx];
+ CCToFaceMatrixBlock& A12 = (*jacobian_)[cellCenterIdx][faceIdx];
+ FaceToCCMatrixBlock& A21 = (*jacobian_)[faceIdx][cellCenterIdx];
+ FaceToFaceMatrixBlock& A22 = (*jacobian_)[faceIdx][faceIdx];
+
+ // set the size of the sub-matrizes
+ A11.setSize(numDofsCC, numDofsCC);
+ A12.setSize(numDofsCC, numDofsFace);
+ A21.setSize(numDofsFace, numDofsCC);
+ A22.setSize(numDofsFace, numDofsFace);
+
+
+ // set occupation pattern of the jacobian
+ Dune::MatrixIndexSet occupationPatternA11;
+ Dune::MatrixIndexSet occupationPatternA12;
+ Dune::MatrixIndexSet occupationPatternA21;
+ Dune::MatrixIndexSet occupationPatternA22;
+ occupationPatternA11.resize(numDofsCC, numDofsCC);
+ occupationPatternA12.resize(numDofsCC, numDofsFace);
+ occupationPatternA21.resize(numDofsFace, numDofsCC);
+ occupationPatternA22.resize(numDofsFace, numDofsFace);
+
+ const auto& connectivityMap = fvGridGeometry().connectivityMap();
+
+ // evaluate the acutal pattern
+ for (const auto& element : elements(fvGridGeometry().gridView()))
+ {
+ // the global index of the element at hand
+ const auto ccGlobalI = fvGridGeometry().elementMapper().index(element);
+
+ for (auto&& ccGlobalJ : connectivityMap(cellCenterIdx, cellCenterIdx, ccGlobalI))
+ occupationPatternA11.add(ccGlobalI, ccGlobalJ);
+ for (auto&& faceGlobalJ : connectivityMap(cellCenterIdx, faceIdx, ccGlobalI))
+ occupationPatternA12.add(ccGlobalI, faceGlobalJ);
+
+ auto fvGeometry = localView(fvGridGeometry());
+ fvGeometry.bindElement(element);
+
+ // loop over sub control faces
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ const auto faceGlobalI = scvf.dofIndex();
+ for (auto&& ccGlobalJ : connectivityMap(faceIdx, cellCenterIdx, scvf.index()))
+ occupationPatternA21.add(faceGlobalI, ccGlobalJ);
+ for (auto&& faceGlobalJ : connectivityMap(faceIdx, faceIdx, scvf.index()))
+ occupationPatternA22.add(faceGlobalI, faceGlobalJ);
+ }
+ }
+
+ occupationPatternA11.exportIdx(A11);
+ occupationPatternA12.exportIdx(A12);
+ occupationPatternA21.exportIdx(A21);
+ occupationPatternA22.exportIdx(A22);
+ }
+
+ /*!
+ * \brief Resizes the residual
+ */
+ void setResidualSize()
+ {
+ (*residual_)[cellCenterIdx].resize(fvGridGeometry().numCellCenterDofs());
+ (*residual_)[faceIdx].resize(fvGridGeometry().numFaceDofs());
+ }
+
+ const Problem& problem() const
+ { return *problem_; }
+
+ const FVGridGeometry& fvGridGeometry() const
+ { return *fvGridGeometry_; }
+
+ const GridView& gridView() const
+ { return fvGridGeometry().gridView(); }
+
+ GridVariables& gridVariables()
+ { return *gridVariables_; }
+
+ const GridVariables& gridVariables() const
+ { return *gridVariables_; }
+
+ JacobianMatrix& jacobian()
+ {
+ if (!residual_)
+ DUNE_THROW(Dune::InvalidStateException, "No jacobian was set.");
+ return *jacobian_;
+ }
+
+ SolutionVector& residual()
+ {
+ if (!residual_)
+ DUNE_THROW(Dune::InvalidStateException, "No residual was set.");
+ return *residual_;
+ }
+
+ const LocalResidual& localResidual() const
+ { return localResidual_; }
+
+private:
+ // reset the residual to 0.0
+ void resetResidual_()
+ {
+ (*residual_)[cellCenterIdx] = 0.0;
+ (*residual_)[faceIdx] = 0.0;
+ }
+
+ // reset the jacobian to 0.0
+ void resetJacobian_()
+ {
+ (*jacobian_)[cellCenterIdx][cellCenterIdx] = 0.0;
+ (*jacobian_)[cellCenterIdx][faceIdx] = 0.0;
+ (*jacobian_)[faceIdx][cellCenterIdx] = 0.0;
+ (*jacobian_)[faceIdx][faceIdx] = 0.0;
+ }
+
+ // pointer to the problem to be solved
+ std::shared_ptr problem_;
+
+ // the finite volume geometry of the grid
+ std::shared_ptr fvGridGeometry_;
+
+ // the variables container for the grid
+ std::shared_ptr gridVariables_;
+
+ // shared pointers to the jacobian matrix and residual
+ std::shared_ptr jacobian_;
+ std::shared_ptr residual_;
+
+ // class computing the residual of an element
+ LocalResidual localResidual_;
+
+ // if this assembler is assembling a time dependent problem
+ bool stationary_;
+};
+
+} // namespace Dumux
+
+#endif
diff --git a/dumux/assembly/staggeredlocalassembler.hh b/dumux/assembly/staggeredlocalassembler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..65dad662104a2c5ff5995191aa2d799d0d383709
--- /dev/null
+++ b/dumux/assembly/staggeredlocalassembler.hh
@@ -0,0 +1,625 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief An assembler for the global linear system for fully implicit models
+ * and cell-centered discretization schemes using Newton's method.
+ */
+#ifndef DUMUX_CC_LOCAL_ASSEMBLER_HH
+#define DUMUX_CC_LOCAL_ASSEMBLER_HH
+
+#include
+#include
+
+#include
+#include
+
+#include
+#include
+
+namespace Dumux {
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief An assembler for the local contributions (per element) to the global
+ * linear system for fully implicit models and cell-centered discretization schemes.
+ */
+template
+class StaggeredLocalAssembler;
+
+
+template
+class StaggeredLocalAssembler
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
+ using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GlobalVolumeVariables);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
+
+ using NumCellCenterEqVector = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+ using NumFaceEqVector = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+
+ using FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
+ using FaceSolution = typename GET_PROP_TYPE(TypeTag, StaggeredFaceSolution);
+
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using PriVarIndices = typename Dumux::PriVarIndices;
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+ using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+ using FaceVariables = typename GET_PROP_TYPE(TypeTag, FaceVariables);
+ using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+
+ static constexpr bool enableGlobalFluxVarsCache = GET_PROP_VALUE(TypeTag, EnableGlobalFluxVariablesCache);
+
+public:
+
+ /*!
+ * \brief Computes the derivatives with respect to the given element and adds them
+ * to the global matrix. The element residual is written into the right hand side.
+ */
+ template
+ static void assembleJacobianAndResidual(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
+ const Element& element, const SolutionVector& curSol)
+ {
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+
+ // get some references for convenience
+ const auto& problem = assembler.problem();
+ auto& localResidual = assembler.localResidual();
+ auto& gridVariables = assembler.gridVariables();
+
+ // prepare the local views
+ auto fvGeometry = localView(assembler.fvGridGeometry());
+ fvGeometry.bind(element);
+
+ auto curElemVolVars = localView(gridVariables.curGridVolVars());
+ curElemVolVars.bind(element, fvGeometry, curSol);
+
+ auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
+ elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
+
+ auto curElemFaceVars = localView(gridVariables.curGridFaceVars());
+ curElemFaceVars.bind(element, fvGeometry, curSol);
+
+ const bool isStationary = localResidual.isStationary();
+ auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
+ auto prevElemFaceVars = localView(gridVariables.prevGridFaceVars());
+ if (!isStationary)
+ {
+ prevElemVolVars.bindElement(element, fvGeometry, localResidual.prevSol());
+ prevElemFaceVars.bindElement(element, fvGeometry, localResidual.prevSol());
+ }
+
+ // for compatibility with box models
+ ElementBoundaryTypes elemBcTypes;
+
+ const auto cellCenterGlobalI = assembler.fvGridGeometry().elementMapper().index(element);
+ res[cellCenterIdx][cellCenterGlobalI] = localResidual.evalCellCenter(problem,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ curElemVolVars,
+ prevElemFaceVars,
+ curElemFaceVars,
+ elemBcTypes,
+ elemFluxVarsCache);
+
+ // treat the local residua of the face dofs:
+ // create a cache to reuse some results for the calculation of the derivatives
+ FaceSolutionVector faceResidualCache;
+ faceResidualCache.resize(fvGeometry.numScvf());
+ faceResidualCache = 0.0;
+
+ for(auto&& scvf : scvfs(fvGeometry))
+ {
+ faceResidualCache[scvf.localFaceIdx()] = localResidual.evalFace(problem,
+ element,
+ fvGeometry,
+ scvf,
+ prevElemVolVars,
+ curElemVolVars,
+ prevElemFaceVars,
+ curElemFaceVars,
+ elemBcTypes,
+ elemFluxVarsCache);
+
+ res[faceIdx][scvf.dofIndex()] += faceResidualCache[scvf.localFaceIdx()] ;
+ }
+
+ // calculate derivatives of all dofs in stencil with respect to the dofs in the element
+ evalPartialDerivatives_(assembler,
+ element,
+ fvGeometry,
+ curSol,
+ prevElemVolVars,
+ curElemVolVars,
+ prevElemFaceVars,
+ curElemFaceVars,
+ elemFluxVarsCache,
+ elemBcTypes,
+ jac,
+ res[cellCenterIdx][cellCenterGlobalI],
+ faceResidualCache);
+
+ }
+
+ /*!
+ * \brief Computes the derivatives with respect to the given element and adds them
+ * to the global matrix.
+ */
+ template
+ static void assemble(Assembler& assembler, JacobianMatrix& jac,
+ const Element& element, const SolutionVector& curSol)
+ {
+ std::cout << "calling wrong \n";
+ assemble_(assembler, jac, element, curSol);
+ }
+
+ /*!
+ * \brief Assemble the residual only
+ */
+ template
+ static void assemble(Assembler& assembler, SolutionVector& res,
+ const Element& element, const SolutionVector& curSol)
+ {
+ std::cout << "calling wrong \n";
+ // using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ // typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ // typename DofTypeIndices::FaceIdx faceIdx;
+
+ // const auto cellCenterGlobalI = assembler.fvGridGeometry().elementMapper().index(element);
+ // res[cellCenterIdx][cellCenterGlobalI] = localResidual.evalCellCenter(problem,
+ // element,
+ // fvGeometry,
+ // prevElemVolVars,
+ // curElemVolVars,
+ // curGlobalFaceVars,
+ // prevGlobalFaceVars,
+ // elemBcTypes,
+ // elemFluxVarsCache)[0];
+
+
+
+ // treat the local residua of the face dofs:
+ // create a cache to reuse some results for the calculation of the derivatives
+ // FaceSolutionVector faceResidualCache;
+ // faceResidualCache.resize(assembler.fvGridGeometry().numScvf());
+ // faceResidualCache = 0.0;
+ //
+ // auto fvGeometry = localView(assembler.fvGridGeometry());
+ // fvGeometry.bind(element);
+ // auto faceResiduals = assembleFace_(assembler, element, curSol);
+ //
+ // for(auto&& scvf : scvfs(fvGeometry))
+ // {
+ // res[faceIdx][scvf.dofIndex()] += faceResiduals[scvf.localFaceIdx()];
+ // faceResidualCache[scvf.localFaceIdx()] = faceResiduals[scvf.localFaceIdx()];
+ // }
+ }
+
+ /*!
+ * \brief Computes the epsilon used for numeric differentiation
+ * for a given value of a primary variable.
+ *
+ * \param priVar The value of the primary variable
+ */
+ static Scalar numericEpsilon(const Scalar priVar)
+ {
+ // define the base epsilon as the geometric mean of 1 and the
+ // resolution of the scalar type. E.g. for standard 64 bit
+ // floating point values, the resolution is about 10^-16 and
+ // the base epsilon is thus approximately 10^-8.
+ /*
+ static const Scalar baseEps
+ = Dumux::geometricMean(std::numeric_limits::epsilon(), 1.0);
+ */
+ static const Scalar baseEps = 1e-10;
+ assert(std::numeric_limits::epsilon()*1e4 < baseEps);
+ // the epsilon value used for the numeric differentiation is
+ // now scaled by the absolute value of the primary variable...
+ return baseEps*(std::abs(priVar) + 1.0);
+ }
+
+protected:
+ template
+ static void evalPartialDerivatives_(Assembler& assembler,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SolutionVector& curSol,
+ const ElementVolumeVariables& prevElemVolVars,
+ ElementVolumeVariables& curElemVolVars,
+ const ElementFaceVariables& prevElemFaceVars,
+ ElementFaceVariables& curElemFaceVars,
+ ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes& elemBcTypes,
+ JacobianMatrix& matrix,
+ const NumCellCenterEqVector& ccResidual,
+ const FaceSolutionVector& faceResidualCache)
+{
+ // compute the derivatives of the cell center residuals with respect to cell center dofs
+ dCCdCC_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, elemFluxVarsCache, elemBcTypes, matrix, ccResidual);
+
+ // compute the derivatives of the cell center residuals with respect to face dofs
+ dCCdFace_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, elemFluxVarsCache, elemBcTypes, matrix, ccResidual);
+
+ // compute the derivatives of the face residuals with respect to cell center dofs
+ dFacedCC_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, elemFluxVarsCache, elemBcTypes, matrix, faceResidualCache);
+
+ // compute the derivatives of the face residuals with respect to face dofs
+ dFacedFace_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, elemFluxVarsCache, elemBcTypes, matrix, faceResidualCache);
+}
+
+/*!
+* \brief Computes the derivatives of the cell center residuals with respect to cell center dofs
+*/
+template
+static void dCCdCC_(Assembler& assembler,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SolutionVector& curSol,
+ const ElementVolumeVariables& prevElemVolVars,
+ ElementVolumeVariables& curElemVolVars,
+ const ElementFaceVariables& prevElemFaceVars,
+ const ElementFaceVariables& curElemFaceVars,
+ ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes& elemBcTypes,
+ JacobianMatrix& matrix,
+ const NumCellCenterEqVector& ccResidual)
+{
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+
+ const auto& problem = assembler.problem();
+ auto& localResidual = assembler.localResidual();
+ auto& gridVariables = assembler.gridVariables();
+
+ // build derivatives with for cell center dofs w.r.t. cell center dofs
+ const auto cellCenterGlobalI = assembler.fvGridGeometry().elementMapper().index(element);
+
+ const auto& connectivityMap = assembler.fvGridGeometry().connectivityMap();
+
+ for(const auto& globalJ : connectivityMap(cellCenterIdx, cellCenterIdx, cellCenterGlobalI))
+ {
+ // get the volVars of the element with respect to which we are going to build the derivative
+ auto&& scvJ = fvGeometry.scv(globalJ);
+ const auto elementJ = fvGeometry.fvGridGeometry().element(globalJ);
+ auto& curVolVars = getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scvJ);
+ VolumeVariables origVolVars(curVolVars);
+
+ for(auto pvIdx : PriVarIndices(cellCenterIdx))
+ {
+ PrimaryVariables priVars(CellCenterPrimaryVariables(curSol[cellCenterIdx][globalJ]),
+ FacePrimaryVariables(0.0));
+
+ const Scalar eps = numericEpsilon(priVars[pvIdx], cellCenterIdx, cellCenterIdx);
+ priVars[pvIdx] += eps;
+ ElementSolutionVector elemSol{std::move(priVars)};
+ curVolVars.update(elemSol, problem, elementJ, scvJ);
+
+ const auto deflectedResidual = localResidual.evalCellCenter(problem, element, fvGeometry, prevElemVolVars, curElemVolVars,
+ prevElemFaceVars, curElemFaceVars,
+ elemBcTypes, elemFluxVarsCache);
+
+ auto partialDeriv = (deflectedResidual - ccResidual);
+ partialDeriv /= eps;
+
+ // update the global jacobian matrix with the current partial derivatives
+ updateGlobalJacobian_(matrix[cellCenterIdx][cellCenterIdx], cellCenterGlobalI, globalJ, pvIdx, partialDeriv);
+
+ // restore the original volVars
+ curVolVars = origVolVars;
+ }
+ }
+}
+
+/*!
+* \brief Computes the derivatives of the cell center residuals with respect to face dofs
+*/
+template
+static void dCCdFace_(Assembler& assembler,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SolutionVector& curSol,
+ const ElementVolumeVariables& prevElemVolVars,
+ const ElementVolumeVariables& curElemVolVars,
+ const ElementFaceVariables& prevElemFaceVars,
+ ElementFaceVariables& curElemFaceVars,
+ ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes& elemBcTypes,
+ JacobianMatrix& matrix,
+ const NumCellCenterEqVector& ccResidual)
+{
+ // build derivatives with for cell center dofs w.r.t. face dofs
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+ const auto& problem = assembler.problem();
+ auto& localResidual = assembler.localResidual();
+ auto& gridVariables = assembler.gridVariables();
+
+ // build derivatives with for cell center dofs w.r.t. cell center dofs
+ const auto cellCenterGlobalI = assembler.fvGridGeometry().elementMapper().index(element);
+
+ for(const auto& scvfJ : scvfs(fvGeometry))
+ {
+ const auto globalJ = scvfJ.dofIndex();
+
+ // get the faceVars of the face with respect to which we are going to build the derivative
+ auto& faceVars = getFaceVarAccess(gridVariables.curGridFaceVars(), curElemFaceVars, scvfJ);
+ const auto origFaceVars(faceVars);
+
+ for(auto pvIdx : PriVarIndices(faceIdx))
+ {
+ FacePrimaryVariables facePriVars(curSol[faceIdx][globalJ]);
+ const Scalar eps = numericEpsilon(facePriVars[pvIdx], cellCenterIdx, faceIdx);
+ facePriVars[pvIdx] += eps;
+
+ faceVars.updateOwnFaceOnly(facePriVars);
+
+ const auto deflectedResidual = localResidual.evalCellCenter(problem, element, fvGeometry,
+ prevElemVolVars, curElemVolVars,
+ prevElemFaceVars, curElemFaceVars,
+ elemBcTypes, elemFluxVarsCache);
+
+ auto partialDeriv = (deflectedResidual - ccResidual);
+ partialDeriv /= eps;
+
+ // update the global jacobian matrix with the current partial derivatives
+ updateGlobalJacobian_(matrix[cellCenterIdx][faceIdx], cellCenterGlobalI, globalJ, pvIdx - Indices::faceOffset, partialDeriv);
+
+ // restore the original faceVars
+ faceVars = origFaceVars;
+ }
+ }
+}
+
+/*!
+* \brief Computes the derivatives of the face residuals with respect to cell center dofs
+*/
+template
+static void dFacedCC_(Assembler& assembler,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SolutionVector& curSol,
+ const ElementVolumeVariables& prevElemVolVars,
+ ElementVolumeVariables& curElemVolVars,
+ const ElementFaceVariables& prevElemFaceVars,
+ const ElementFaceVariables& curElemFaceVars,
+ ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes& elemBcTypes,
+ JacobianMatrix& matrix,
+ const FaceSolutionVector& cachedResidual)
+{
+ for(auto&& scvf : scvfs(fvGeometry))
+ {
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+ const auto& problem = assembler.problem();
+ auto& localResidual = assembler.localResidual();
+ auto& gridVariables = assembler.gridVariables();
+ const auto& connectivityMap = assembler.fvGridGeometry().connectivityMap();
+
+ // set the actual dof index
+ const auto faceGlobalI = scvf.dofIndex();
+
+ // build derivatives with for face dofs w.r.t. cell center dofs
+ for(const auto& globalJ : connectivityMap(faceIdx, cellCenterIdx, scvf.index()))
+ {
+ // get the volVars of the element with respect to which we are going to build the derivative
+ auto&& scvJ = fvGeometry.scv(globalJ);
+ const auto elementJ = fvGeometry.fvGridGeometry().element(globalJ);
+ auto& curVolVars = getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scvJ);
+ VolumeVariables origVolVars(curVolVars);
+
+ for(auto pvIdx : PriVarIndices(cellCenterIdx))
+ {
+ PrimaryVariables priVars(CellCenterPrimaryVariables(curSol[cellCenterIdx][globalJ]),
+ FacePrimaryVariables(0.0));
+
+ const Scalar eps = numericEpsilon(priVars[pvIdx], faceIdx, cellCenterIdx);
+ priVars[pvIdx] += eps;
+ ElementSolutionVector elemSol{std::move(priVars)};
+ curVolVars.update(elemSol, problem, elementJ, scvJ);
+
+ const auto deflectedResidual = localResidual.evalFace(problem, element, fvGeometry, scvf,
+ prevElemVolVars, curElemVolVars,
+ prevElemFaceVars, curElemFaceVars,
+ elemBcTypes, elemFluxVarsCache);
+
+ auto partialDeriv = (deflectedResidual - cachedResidual[scvf.localFaceIdx()]);
+ partialDeriv /= eps;
+
+ // update the global jacobian matrix with the current partial derivatives
+ updateGlobalJacobian_(matrix[faceIdx][cellCenterIdx], faceGlobalI, globalJ, pvIdx, partialDeriv);
+
+ // restore the original volVars
+ curVolVars = origVolVars;
+ }
+ }
+ }
+}
+
+/*!
+* \brief Computes the derivatives of the face residuals with respect to cell center dofs
+*/
+template
+static void dFacedFace_(Assembler& assembler,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SolutionVector& curSol,
+ const ElementVolumeVariables& prevElemVolVars,
+ const ElementVolumeVariables& curElemVolVars,
+ const ElementFaceVariables& prevElemFaceVars,
+ ElementFaceVariables& curElemFaceVars,
+ ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes& elemBcTypes,
+ JacobianMatrix& matrix,
+ const FaceSolutionVector& cachedResidual)
+{
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+ const auto& problem = assembler.problem();
+ auto& localResidual = assembler.localResidual();
+ const auto& connectivityMap = assembler.fvGridGeometry().connectivityMap();
+ auto& gridVariables = assembler.gridVariables();
+
+ for(auto&& scvf : scvfs(fvGeometry))
+ {
+ // set the actual dof index
+ const auto faceGlobalI = scvf.dofIndex();
+
+ // build derivatives with for face dofs w.r.t. cell center dofs
+ for(const auto& globalJ : connectivityMap(faceIdx, faceIdx, scvf.index()))
+ {
+ // get the faceVars of the face with respect to which we are going to build the derivative
+ auto& faceVars = getFaceVarAccess(gridVariables.curGridFaceVars(), curElemFaceVars, scvf);
+ const auto origFaceVars(faceVars);
+
+ for(auto pvIdx : PriVarIndices(faceIdx))
+ {
+ auto faceSolution = FaceSolution(scvf, curSol[faceIdx], assembler.fvGridGeometry());
+
+ const Scalar eps = numericEpsilon(faceSolution[globalJ][pvIdx], faceIdx, faceIdx);
+
+ faceSolution[globalJ][pvIdx] += eps;
+ faceVars.update(faceSolution, problem, element, fvGeometry, scvf);
+
+ const auto deflectedResidual = localResidual.evalFace(problem, element, fvGeometry, scvf,
+ prevElemVolVars, curElemVolVars,
+ prevElemFaceVars, curElemFaceVars,
+ elemBcTypes, elemFluxVarsCache);
+
+ auto partialDeriv = (deflectedResidual - cachedResidual[scvf.localFaceIdx()]);
+ partialDeriv /= eps;
+
+ // update the global jacobian matrix with the current partial derivatives
+ updateGlobalJacobian_(matrix[faceIdx][faceIdx], faceGlobalI, globalJ, pvIdx - Indices::faceOffset, partialDeriv);
+
+ // restore the original faceVars
+ faceVars = origFaceVars;
+ }
+ }
+ }
+}
+
+
+
+static Scalar numericEpsilon(const Scalar priVar, const int idx1, const int idx2)
+{
+ // define the base epsilon as the geometric mean of 1 and the
+ // resolution of the scalar type. E.g. for standard 64 bit
+ // floating point values, the resolution is about 10^-16 and
+ // the base epsilon is thus approximately 10^-8.
+ /*
+ static const Scalar baseEps
+ = Dumux::geometricMean(std::numeric_limits::epsilon(), 1.0);
+ */
+ using BaseEpsilon = typename GET_PROP(TypeTag, BaseEpsilon);
+ const std::array, 2> baseEps_ = BaseEpsilon::getEps();
+
+
+ static const Scalar baseEps = baseEps_[idx1][idx2];
+ assert(std::numeric_limits::epsilon()*1e4 < baseEps);
+ // the epsilon value used for the numeric differentiation is
+ // now scaled by the absolute value of the primary variable...
+ return baseEps*(std::abs(priVar) + 1.0);
+}
+
+
+/*!
+ * \brief Updates the current global Jacobian matrix with the
+ * partial derivatives of all equations in regard to the
+ * primary variable 'pvIdx' at dof 'col'. Specialization for cc methods.
+ */
+template
+static void updateGlobalJacobian_(SubMatrix& matrix,
+ const int globalI,
+ const int globalJ,
+ const int pvIdx,
+ const CCOrFacePrimaryVariables &partialDeriv)
+{
+ for (int eqIdx = 0; eqIdx < partialDeriv.size(); eqIdx++)
+ {
+ // A[i][col][eqIdx][pvIdx] is the rate of change of
+ // the residual of equation 'eqIdx' at dof 'i'
+ // depending on the primary variable 'pvIdx' at dof
+ // 'col'.
+
+ assert(pvIdx >= 0);
+ assert(eqIdx < matrix[globalI][globalJ].size());
+ assert(pvIdx < matrix[globalI][globalJ][eqIdx].size());
+ matrix[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ Valgrind::CheckDefined(matrix[globalI][globalJ][eqIdx][pvIdx]);
+ }
+}
+
+
+
+
+
+private:
+ template
+ static typename std::enable_if::type
+ getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
+ { return elemVolVars[scv]; }
+
+ template
+ static typename std::enable_if::type
+ getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
+ { return gridVolVars.volVars(scv); }
+
+ template
+ static typename std::enable_if::type
+ getFaceVarAccess(GlobalFaceVars& gridFaceVars, ElementFaceVariables& elemFaceVars, const SubControlVolumeFace& scvf)
+ { return elemFaceVars[scvf]; }
+
+ template
+ static typename std::enable_if::type
+ getFaceVarAccess(GlobalFaceVars& gridFaceVars, ElementFaceVariables& elemFaceVars, const SubControlVolumeFace& scvf)
+ { return gridFaceVars.faceVars(scvf.index()); }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh
index 5b79cae1e2e6bbd18de694e8887539754146c77a..cb03e8348de6ba7a5988adbebbd6b9a33283ba84 100644
--- a/dumux/common/properties.hh
+++ b/dumux/common/properties.hh
@@ -157,5 +157,23 @@ NEW_PROP_TAG(EvaluatePermeabilityAtScvfIP);
NEW_PROP_TAG(Chemistry); //!< The chemistry class with which solves equlibrium reactions
NEW_PROP_TAG(NumMajorComponents); //!< Number of major fluid components which are considered in the calculation of the phase density
NEW_PROP_TAG(SetMoleFractionsForWettingPhase); //!< Set the mole fraction in the wetting or non-wetting phase
+
+/////////////////////////////////////////////////////////////
+// Properties used by the staggered-grid discretization method
+/////////////////////////////////////////////////////////////
+NEW_PROP_TAG(NumEqCellCenter); //! The number of equations for cell-centered dofs
+NEW_PROP_TAG(NumEqFace); //! The number of equations for face dofs
+NEW_PROP_TAG(CellCenterSolutionVector); //! The solution vector type for cell-centered dofs
+NEW_PROP_TAG(FaceSolutionVector); //! The solution vector type for face dofs
+NEW_PROP_TAG(GlobalFaceVars); //! Class containing face-related data
+NEW_PROP_TAG(CellCenterPrimaryVariables); //! The primary variables container type for cell-centered dofs
+NEW_PROP_TAG(FacePrimaryVariables); //! The primary variables container type for face dofs
+NEW_PROP_TAG(IntersectionMapper); //! Specifies the intersection mapper
+NEW_PROP_TAG(DofTypeIndices); //! Specifies index types for accessing the multi type block vectors/matrices
+NEW_PROP_TAG(StaggeredGeometryHelper); //! Specifies a helper class for the staggered grid geometry
+NEW_PROP_TAG(StaggeredPrimaryVariables); //! The hybrid primary variables container type
+NEW_PROP_TAG(BaseEpsilon); //! A base epsilon for numerical differentiation, can contain multiple values
+NEW_PROP_TAG(FaceVariables); //! Class containing local face-related data
+NEW_PROP_TAG(BoundaryValues); //! Class containing local boundary data
}
}
diff --git a/dumux/common/staggeredfvproblem.hh b/dumux/common/staggeredfvproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..384bbcd3cc62d210a4c12ef3400a86185bc7d51d
--- /dev/null
+++ b/dumux/common/staggeredfvproblem.hh
@@ -0,0 +1,159 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for all problems
+ */
+#ifndef DUMUX_STAGGERD_FV_PROBLEM_HH
+#define DUMUX_STAGGERD_FV_PROBLEM_HH
+
+#include
+#include
+
+namespace Dumux
+{
+/*!
+ * \ingroup Problems
+ * \brief Base class for all finite-volume problems
+ *
+ * \note All quantities (regarding the units) are specified assuming a
+ * three-dimensional world. Problems discretized using 2D grids
+ * are assumed to be extruded by \f$1 m\f$ and 1D grids are assumed
+ * to have a cross section of \f$1m \times 1m\f$.
+ */
+template
+class StaggeredFVProblem : public FVProblem
+{
+ using ParentType = FVProblem;
+ using Implementation = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using VertexMapper = typename GET_PROP_TYPE(TypeTag, VertexMapper);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using PointSource = typename GET_PROP_TYPE(TypeTag, PointSource);
+ using PointSourceHelper = typename GET_PROP_TYPE(TypeTag, PointSourceHelper);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using InitialValues = typename GET_PROP_TYPE(TypeTag, BoundaryValues);
+
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Vertex = typename GridView::template Codim::Entity;
+ using CoordScalar = typename GridView::ctype;
+ using GlobalPosition = Dune::FieldVector;
+
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+
+ // using GridAdaptModel = ImplicitGridAdapt;
+
+public:
+ /*!
+ * \brief Constructor
+ *
+ * \param gridView The simulation's idea about physical space
+ */
+ StaggeredFVProblem(std::shared_ptr fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ { }
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ // Throw an exception (there is no reasonable default value
+ // for initial values)
+ DUNE_THROW(Dune::InvalidStateException,
+ "The problem does not provide "
+ "an initial() or an initialAtPos() method.");
+ }
+
+ /*!
+ * \brief Evaluate the initial value for
+ * an element (for cell-centered models)
+ * or vertex (for box / vertex-centered models)
+ *
+ * \param entity The dof entity (element or vertex)
+ */
+ template
+ InitialValues initial(const Entity& entity) const
+ {
+ // static_assert(int(Entity::codimension) == 0 || int(Entity::codimension) == dim, "Entity must be element or vertex");
+ return asImp_().initialAtPos(entity.center());
+ }
+
+ /*!
+ * \brief Applies the initial solution for all degrees of freedom of the grid.
+ *
+ */
+ void applyInitialSolution(SolutionVector& sol) const
+ {
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+
+ // loop over sub control volumes
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ // let the problem do the dirty work of nailing down
+ // the initial solution.
+ auto initPriVars = asImp_().initial(scv)[cellCenterIdx];
+ auto dofIdxGlobal = scv.dofIndex();
+ sol[cellCenterIdx][dofIdxGlobal] += initPriVars;
+ }
+
+ // loop over faces
+ for(auto&& scvf : scvfs(fvGeometry))
+ {
+ auto initPriVars = asImp_().initial(scvf)[faceIdx][scvf.directionIndex()];
+ sol[faceIdx][scvf.dofIndex()] = initPriVars;
+ }
+ }
+ }
+
+protected:
+ //! Returns the implementation of the problem (i.e. static polymorphism)
+ Implementation &asImp_()
+ { return *static_cast(this); }
+
+ //! \copydoc asImp_()
+ const Implementation &asImp_() const
+ { return *static_cast(this); }
+
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/implicit/staggered/assemblymap.hh b/dumux/discretization/staggered/connectivitymap.hh
similarity index 83%
rename from dumux/implicit/staggered/assemblymap.hh
rename to dumux/discretization/staggered/connectivitymap.hh
index 45331e432b4a90986a1eb1e197e625b90f069f33..c71587679a5aa483b30114db32352ed029a7cb98 100644
--- a/dumux/implicit/staggered/assemblymap.hh
+++ b/dumux/discretization/staggered/connectivitymap.hh
@@ -22,22 +22,22 @@
* that contribute to the derivative calculation. This is used for
* finite-volume schemes with symmetric sparsity pattern in the global matrix.
*/
-#ifndef DUMUX_STAGGERED_ASSEMBLY_MAP_HH
-#define DUMUX_STAGGERED_ASSEMBLY_MAP_HH
+#ifndef DUMUX_STAGGERED_CONNECTIVITY_MAP_HH
+#define DUMUX_STAGGERED_CONNECTIVITY_MAP_HH
-#include
-
-#include
+#include
+#include
namespace Dumux
{
template
-class StaggeredAssemblyMap
+class StaggeredConnectivityMap
{
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 FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using IndexType = typename GridView::IndexSet::IndexType;
@@ -62,11 +62,11 @@ public:
*
* \param problem The problem which we want to simulate.
*/
- void init(const Problem& problem)
+ void update(const FVGridGeometry& fvGridGeometry)
{
- const auto numDofsCC = problem.model().numCellCenterDofs();
- const auto numDofsFace = problem.model().numFaceDofs();
- const auto numBoundaryFacets = problem.model().fvGridGeometry().numBoundaryScvf();
+ const auto numDofsCC = fvGridGeometry.gridView().size(0);
+ const auto numDofsFace = fvGridGeometry.gridView().size(1);
+ const auto numBoundaryFacets = fvGridGeometry.numBoundaryScvf();
cellCenterToCellCenterMap_.resize(numDofsCC);
cellCenterToFaceMap_.resize(numDofsCC);
faceToCellCenterMap_.resize(2*numDofsFace - numBoundaryFacets);
@@ -78,22 +78,22 @@ public:
fullfaceToFaceStencils.resize(numDofsFace);
FluxVariables fluxVars;
- for(auto&& element: elements(problem.gridView()))
+ for(auto&& element: elements(fvGridGeometry.gridView()))
{
// restrict the FvGeometry locally and bind to the element
- auto fvGeometry = localView(problem.model().fvGridGeometry());
+ auto fvGeometry = localView(fvGridGeometry);
fvGeometry.bindElement(element);
// loop over sub control faces
for (auto&& scvf : scvfs(fvGeometry))
{
- const auto dofIdxCellCenter = problem.elementMapper().index(element);
- fluxVars.computeCellCenterToCellCenterStencil(cellCenterToCellCenterMap_[dofIdxCellCenter], problem, element, fvGeometry, scvf);
- fluxVars.computeCellCenterToFaceStencil(cellCenterToFaceMap_[dofIdxCellCenter], problem, element, fvGeometry, scvf);
+ const auto dofIdxCellCenter = fvGridGeometry.elementMapper().index(element);
+ fluxVars.computeCellCenterToCellCenterStencil(cellCenterToCellCenterMap_[dofIdxCellCenter], element, fvGeometry, scvf);
+ fluxVars.computeCellCenterToFaceStencil(cellCenterToFaceMap_[dofIdxCellCenter], element, fvGeometry, scvf);
const auto scvfIdx = scvf.index();
- fluxVars.computeFaceToCellCenterStencil(faceToCellCenterMap_[scvfIdx],problem, fvGeometry, scvf);
- fluxVars.computeFaceToFaceStencil(faceToFaceMap_[scvfIdx],problem, fvGeometry, scvf);
+ fluxVars.computeFaceToCellCenterStencil(faceToCellCenterMap_[scvfIdx], fvGeometry, scvf);
+ fluxVars.computeFaceToFaceStencil(faceToFaceMap_[scvfIdx], fvGeometry, scvf);
}
}
}
diff --git a/dumux/discretization/staggered/elementfacevariables.hh b/dumux/discretization/staggered/elementfacevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b9ab9e87f191d906fdb7a6a6063c4d335c54a1c9
--- /dev/null
+++ b/dumux/discretization/staggered/elementfacevariables.hh
@@ -0,0 +1,175 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief The face variables class for free flow staggered grid models
+ */
+#ifndef DUMUX_DISCRETIZATION_STAGGERED_ELEMENTFACEVARIABLES_HH
+#define DUMUX_DISCRETIZATION_STAGGERED_ELEMENTFACEVARIABLES_HH
+
+#include
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief Base class for the face variables vector
+ */
+template
+class StaggeredElementFaceVariables
+{};
+
+
+template
+class StaggeredElementFaceVariables
+{
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
+ using FaceVariables = typename GET_PROP_TYPE(TypeTag, FaceVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using IndexType = typename GridView::IndexSet::IndexType;
+
+public:
+
+ StaggeredElementFaceVariables(const GlobalFaceVars& globalFacesVars) : globalFaceVarsPtr_(&globalFacesVars) {}
+
+ const FaceVariables& operator [](const SubControlVolumeFace& scvf) const
+ { return globalFaceVars().faceVars(scvf.index()); }
+
+ // operator for the access with an index
+ // needed for cc methods for the access to the boundary volume variables
+ const FaceVariables& operator [](const IndexType scvfIdx) const
+ { return globalFaceVars().faceVars(scvfIdx); }
+
+
+ //! For compatibility reasons with the case of not storing the face vars.
+ //! function to be called before assembling an element, preparing the vol vars within the stencil
+ void bind(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SolutionVector& sol)
+ {}
+
+ // Binding of an element, prepares only the face variables of the element
+ // specialization for Staggered models
+ void bindElement(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SolutionVector& sol)
+ {}
+
+
+ //! The global volume variables object we are a restriction of
+ const GlobalFaceVars& globalFaceVars() const
+ { return *globalFaceVarsPtr_; }
+
+
+private:
+ const GlobalFaceVars* globalFaceVarsPtr_;
+};
+
+template
+class StaggeredElementFaceVariables
+{
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
+ using FaceVariables = typename GET_PROP_TYPE(TypeTag, FaceVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using IndexType = typename GridView::IndexSet::IndexType;
+
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+public:
+
+ StaggeredElementFaceVariables(const GlobalFaceVars& globalFacesVars) : globalFaceVarsPtr_(&globalFacesVars) {}
+
+ const FaceVariables& operator [](const SubControlVolumeFace& scvf) const
+ { return faceVariables_[scvf.localFaceIdx()]; }
+
+ // operator for the access with an index
+ const FaceVariables& operator [](const IndexType scvfIdx) const
+ { return faceVariables_[getLocalIdx_(scvfIdx)]; }
+
+ FaceVariables& operator [](const SubControlVolumeFace& scvf)
+ { return faceVariables_[scvf.localFaceIdx()]; }
+
+ // operator for the access with an index
+ FaceVariables& operator [](const IndexType scvfIdx)
+ { return faceVariables_[getLocalIdx_(scvfIdx)]; }
+
+ //! For compatibility reasons with the case of not storing the vol vars.
+ //! function to be called before assembling an element, preparing the vol vars within the stencil
+ void bind(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SolutionVector& sol)
+ {
+ faceVariables_.resize(fvGeometry.numScvf());
+ faceVarIndices_.resize(fvGeometry.numScvf());
+
+ for(auto&& scvf : scvfs(fvGeometry))
+ {
+ faceVariables_[scvf.localFaceIdx()].update(sol[faceIdx], globalFaceVars().problem(), element, fvGeometry, scvf);
+ faceVarIndices_[scvf.localFaceIdx()] = scvf.index();
+ }
+ }
+
+ // Binding of an element, prepares only the face variables of the element
+ // specialization for Staggered models
+ void bindElement(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SolutionVector& sol)
+ {
+ faceVariables_.resize(fvGeometry.numScvf());
+ faceVarIndices_.resize(fvGeometry.numScvf());
+
+ for(auto&& scvf : scvfs(fvGeometry))
+ {
+ faceVariables_[scvf.localFaceIdx()].updateOwnFaceOnly(sol[faceIdx][scvf.dofIndex()]);
+ faceVarIndices_[scvf.localFaceIdx()] = scvf.index();
+ }
+ }
+
+ //! The global volume variables object we are a restriction of
+ const GlobalFaceVars& globalFaceVars() const
+ { return *globalFaceVarsPtr_; }
+
+private:
+
+ const int getLocalIdx_(const int scvfIdx) const
+ {
+ auto it = std::find(faceVarIndices_.begin(), faceVarIndices_.end(), scvfIdx);
+ assert(it != faceVarIndices_.end() && "Could not find the current face variables for scvfIdx!");
+ return std::distance(faceVarIndices_.begin(), it);
+ }
+
+ const GlobalFaceVars* globalFaceVarsPtr_;
+ std::vector faceVarIndices_;
+ std::vector faceVariables_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/staggered/elementfluxvariablescache.hh b/dumux/discretization/staggered/elementfluxvariablescache.hh
index 5feb420cea47a244087e541603edce95dcf22c73..adf2e895cefcf40fee06d28c5d556ee3f1ca30c3 100644
--- a/dumux/discretization/staggered/elementfluxvariablescache.hh
+++ b/dumux/discretization/staggered/elementfluxvariablescache.hh
@@ -23,7 +23,7 @@
#ifndef DUMUX_DISCRETIZATION_STAGGERED_ELEMENT_FLUXVARSCACHE_HH
#define DUMUX_DISCRETIZATION_STAGGERED_ELEMENT_FLUXVARSCACHE_HH
-#include
+#include
namespace Dumux
{
diff --git a/dumux/discretization/staggered/elementvolumevariables.hh b/dumux/discretization/staggered/elementvolumevariables.hh
index 0ada499689039043105ab4dbc8ba6542af3f6e9b..9e2cc021ceed2da95fcc04fa9cc5385fa7c8c1ea 100644
--- a/dumux/discretization/staggered/elementvolumevariables.hh
+++ b/dumux/discretization/staggered/elementvolumevariables.hh
@@ -23,7 +23,7 @@
#ifndef DUMUX_DISCRETIZATION_STAGGERED_ELEMENT_VOLUMEVARIABLES_HH
#define DUMUX_DISCRETIZATION_STAGGERED_ELEMENT_VOLUMEVARIABLES_HH
-#include
+#include
namespace Dumux
{
@@ -93,8 +93,11 @@ class StaggeredElementVolumeVariables::Entity;
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+ static constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
+
public:
//! Constructor
@@ -115,30 +124,35 @@ public:
const FVElementGeometry& fvGeometry,
const SolutionVector& sol)
{
- auto eIdx = globalVolVars().problem_().elementMapper().index(element);
+ clear();
- // stencil information
- const auto& neighborStencil = globalVolVars().problem_().model().stencils(element).neighborStencil();
- const auto numDofs = neighborStencil.size() + 1;
+ const auto& problem = globalVolVars().problem();
+ const auto& fvGridGeometry = fvGeometry.fvGridGeometry();
+ const auto globalI = fvGridGeometry.elementMapper().index(element);
+ const auto map = fvGridGeometry.connectivityMap();
+ const auto& connectivityMapI = map(cellCenterIdx, cellCenterIdx, globalI);
+ const auto numDofs = connectivityMapI.size();
+
+ auto&& scvI = fvGeometry.scv(globalI);
// resize local containers to the required size (for internal elements)
volumeVariables_.resize(numDofs);
volVarIndices_.resize(numDofs);
int localIdx = 0;
- // update the volume variables of the element at hand
- auto&& scvI = fvGeometry.scv(eIdx);
- volumeVariables_[localIdx].update(sol[eIdx], globalVolVars().problem_(), element, scvI);
- volVarIndices_[localIdx] = scvI.index();
- ++localIdx;
-
- // Update the volume variables of the neighboring elements
- for (auto globalJ : neighborStencil)
+ // Update the volume variables of the element at hand and the neighboring elements
+ for (auto globalJ : connectivityMapI)
{
- const auto& elementJ = fvGeometry.fvGridGeometry().element(globalJ);
+ const auto& elementJ = fvGridGeometry.element(globalJ);
auto&& scvJ = fvGeometry.scv(globalJ);
- volumeVariables_[localIdx].update(sol[globalJ], globalVolVars().problem_(), elementJ, scvJ);
- volVarIndices_[localIdx] = scvJ.index();
+ PrimaryVariables priVars(0.0);
+ priVars[cellCenterIdx] = sol[cellCenterIdx][globalJ];
+ ElementSolution elemSol{std::move(priVars)};
+ volumeVariables_[localIdx].update(elemSol,
+ problem,
+ elementJ,
+ scvJ);
+ volVarIndices_[localIdx] = scvJ.dofIndex();
++localIdx;
}
@@ -149,18 +163,33 @@ public:
if (!scvf.boundary())
continue;
- // check if boundary is a pure dirichlet boundary
- const auto bcTypes = globalVolVars().problem_().boundaryTypes(element, scvf);
- if (bcTypes.hasOnlyDirichlet())
- {
- const auto dirichletPriVars = globalVolVars().problem_().dirichlet(element, scvf);
+ const auto bcTypes = problem.boundaryTypes(element, scvf);
+
+ PrimaryVariables boundaryPriVars(0.0);
- volumeVariables_.resize(localIdx+1);
- volVarIndices_.resize(localIdx+1);
- volumeVariables_[localIdx].update(dirichletPriVars, globalVolVars().problem_(), element, scvI);
- volVarIndices_[localIdx] = scvf.outsideScvIdx();
- ++localIdx;
+ for(int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
+ {
+ if(bcTypes.isDirichlet(eqIdx) || bcTypes.isDirichletCell(eqIdx))
+ boundaryPriVars[cellCenterIdx][eqIdx] = problem.dirichlet(element, scvf)[cellCenterIdx][eqIdx];
+ else if(bcTypes.isNeumann(eqIdx) || bcTypes.isOutflow(eqIdx) || bcTypes.isSymmetry())
+ boundaryPriVars[cellCenterIdx][eqIdx] = sol[cellCenterIdx][scvf.insideScvIdx()][eqIdx];
+ //TODO: this assumes a zero-gradient for e.g. the pressure on the boundary
+ // could be made more general by allowing a non-zero-gradient, provided in problem file
+ else
+ if(eqIdx == Indices::pressureIdx)
+ DUNE_THROW(Dune::InvalidStateException, "Face at: " << scvf.center() << " has neither Dirichlet nor Neumann BC.");
}
+
+ volumeVariables_.resize(localIdx+1);
+ volVarIndices_.resize(localIdx+1);
+
+ ElementSolution elemSol{std::move(boundaryPriVars)};
+ volumeVariables_[localIdx].update(elemSol,
+ problem,
+ element,
+ scvI);
+ volVarIndices_[localIdx] = scvf.outsideScvIdx();
+ ++localIdx;
}
}
@@ -170,13 +199,21 @@ public:
const FVElementGeometry& fvGeometry,
const SolutionVector& sol)
{
- auto eIdx = globalVolVars().problem_().elementMapper().index(element);
+ clear();
+
+ const auto eIdx = fvGeometry.fvGridGeometry().elementMapper().index(element);
volumeVariables_.resize(1);
volVarIndices_.resize(1);
// update the volume variables of the element
auto&& scv = fvGeometry.scv(eIdx);
- volumeVariables_[0].update(sol[eIdx], globalVolVars().problem_(), element, scv);
+ PrimaryVariables priVars(0.0);
+ priVars[cellCenterIdx] = sol[cellCenterIdx][eIdx];
+ ElementSolution elemSol{std::move(priVars)};
+ volumeVariables_[0].update(elemSol,
+ globalVolVars().problem(),
+ element,
+ scv);
volVarIndices_[0] = scv.dofIndex();
}
@@ -196,6 +233,13 @@ public:
const GlobalVolumeVariables& globalVolVars() const
{ return *globalVolVarsPtr_; }
+ //! Clear all local storage
+ void clear()
+ {
+ volVarIndices_.clear();
+ volumeVariables_.clear();
+ }
+
private:
const GlobalVolumeVariables* globalVolVarsPtr_;
diff --git a/dumux/discretization/staggered/facesolution.hh b/dumux/discretization/staggered/facesolution.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d36ea605aa5dd38704e9f01a6c16e3e82e6f8ccd
--- /dev/null
+++ b/dumux/discretization/staggered/facesolution.hh
@@ -0,0 +1,92 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief The global volume variables class for cell centered models
+ */
+#ifndef DUMUX_DISCRETIZATION_STAGGERED_FACE_SOLUTION_HH
+#define DUMUX_DISCRETIZATION_STAGGERED_FACE_SOLUTION_HH
+
+#include
+#include
+
+namespace Dumux
+{
+
+template
+class StaggeredFaceSolution
+{
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
+ using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+public:
+
+ StaggeredFaceSolution(const SubControlVolumeFace& scvf, const FaceSolutionVector& sol,
+ const FVGridGeometry& fvGridGeometry)
+ {
+ const auto& connectivityMap = fvGridGeometry.connectivityMap();
+ const auto& stencil = connectivityMap(faceIdx, faceIdx, scvf.index());
+
+ facePriVars_.reserve(stencil.size());
+ map_.reserve(stencil.size());
+
+ for(const auto dofJ : stencil)
+ {
+ map_.push_back(dofJ);
+ facePriVars_.push_back(sol[dofJ]);
+ }
+ }
+
+ //! bracket operator const access
+ template
+ const FacePrimaryVariables& operator [](IndexType globalFaceDofIdx) const
+ {
+ const auto pos = std::find(map_.begin(), map_.end(), globalFaceDofIdx);
+ assert (pos != map_.end());
+ return facePriVars_[pos - map_.begin()];
+ }
+
+ //! bracket operator
+ template
+ FacePrimaryVariables& operator [](IndexType globalFaceDofIdx)
+ {
+ const auto pos = std::find(map_.begin(), map_.end(), globalFaceDofIdx);
+ assert (pos != map_.end());
+ return facePriVars_[pos - map_.begin()];
+ }
+
+
+private:
+
+ std::vector facePriVars_;
+ std::vector map_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/staggered/freeflow/facevariables.hh b/dumux/discretization/staggered/freeflow/facevariables.hh
index f9591f5101453f9afe0099efc4b7e13a72187293..337b6e419d801d6615d80ff09232be05117b7f4b 100644
--- a/dumux/discretization/staggered/freeflow/facevariables.hh
+++ b/dumux/discretization/staggered/freeflow/facevariables.hh
@@ -18,34 +18,171 @@
*****************************************************************************/
/*!
* \file
- * \brief The face variables class for free flow staggered grid models
+ * \brief The face variables class for free flow staggered grid models.
+ * Contains all relevant velocities for the assembly of the momentum balance.
*/
#ifndef DUMUX_DISCRETIZATION_STAGGERED_FREEFLOW_FACEVARIABLES_HH
#define DUMUX_DISCRETIZATION_STAGGERED_FREEFLOW_FACEVARIABLES_HH
-#include
+#include
namespace Dumux
{
+
+namespace Properties
+{
+ NEW_PROP_TAG(StaggeredFaceSolution);
+}
+
template
class StaggeredFaceVariables
{
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 FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+
+ static constexpr int dimWorld = GridView::dimensionworld;
+ static constexpr int numPairs = (dimWorld == 2) ? 2 : 4;
+
+ using GlobalPosition = Dune::FieldVector;
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
public:
- void update(const FacePrimaryVariables &facePrivars)
+
+ /*!
+ * \brief Partial update of the face variables. Only the face itself is considered.
+ *
+ * \param priVars The face-specific primary variales
+ */
+ void updateOwnFaceOnly(const FacePrimaryVariables& priVars)
+ {
+ velocitySelf_ = priVars[0];
+ }
+
+ /*!
+ * \brief Complete update of the face variables (i.e. velocities for free flow)
+ * for a given face
+ *
+ * \param faceSol The face-specific solution vector
+ * \param problem The problem
+ * \param element The element
+ * \param fvGeometry The finite-volume geometry
+ * \param scvf The sub-control volume face of interest
+ */
+ template
+ void update(const SolVector& faceSol,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf)
+ {
+ velocitySelf_ = faceSol[scvf.dofIndex()];
+ velocityOpposite_ = faceSol[scvf.dofIndexOpposingFace()];
+
+ // lambda to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
+ auto makeGhostFace = [](const auto& pos)
+ {
+ return SubControlVolumeFace(pos, std::vector{0,0});
+ };
+
+ // lambda to check whether there is a parallel face neighbor
+ auto hasParallelNeighbor = [](const auto& subFaceData)
+ {
+ return subFaceData.outerParallelFaceDofIdx >= 0;
+ };
+
+ // lambda to check whether there is a normal face neighbor
+ auto hasNormalNeighbor = [](const auto& subFaceData)
+ {
+ return subFaceData.normalPair.second >= 0;
+ };
+
+ // handle all sub faces
+ for(int i = 0; i < scvf.pairData().size(); ++i)
+ {
+ const auto& subFaceData = scvf.pairData(i);
+
+ // treat the velocities normal to the face
+ velocityNormalInside_[i] = faceSol[subFaceData.normalPair.first];
+
+ if(hasNormalNeighbor(subFaceData))
+ {
+ velocityNormalOutside_[i] = faceSol[subFaceData.normalPair.second];
+ }
+ else
+ {
+ const auto& normalFace = fvGeometry.scvf(scvf.insideScvIdx(), subFaceData.localNormalFaceIdx);
+ const auto normalDirIdx = normalFace.directionIndex();
+ velocityNormalOutside_[i] = problem.dirichlet(element, makeGhostFace(subFaceData.virtualOuterNormalFaceDofPos))[faceIdx][normalDirIdx];
+ }
+
+ // treat the velocity parallel to the face
+ velocityParallel_[i] = hasParallelNeighbor(subFaceData) ?
+ velocityParallel_[i] = faceSol[subFaceData.outerParallelFaceDofIdx] :
+ problem.dirichlet(element, makeGhostFace(subFaceData.virtualOuterParallelFaceDofPos))[faceIdx][scvf.directionIndex()];
+ }
+ }
+
+ /*!
+ * \brief Returns the velocity at the face itself
+ */
+ Scalar velocitySelf() const
{
- velocity_ = facePrivars[0];
+ return velocitySelf_;
}
- Scalar velocity() const
+ /*!
+ * \brief Returns the velocity at the opposing face
+ */
+ Scalar velocityOpposite() const
{
- return velocity_;
+ return velocityOpposite_;
}
+ /*!
+ * \brief Returns the velocity at the parallel face
+ *
+ * \param localSubFaceIdx The local index of the subface
+ */
+ Scalar velocityParallel(const int localSubFaceIdx) const
+ {
+ return velocityParallel_[localSubFaceIdx];
+ }
+
+ /*!
+ * \brief Returns the velocity at the inner normal face
+ *
+ * \param localSubFaceIdx The local index of the subface
+ */
+ Scalar velocityNormalInside(const int localSubFaceIdx) const
+ {
+ return velocityNormalInside_[localSubFaceIdx];
+ }
+
+ /*!
+ * \brief Returns the velocity at the outer normal face
+ *
+ * \param localSubFaceIdx The local index of the subface
+ */
+ Scalar velocityNormalOutside(const int localSubFaceIdx) const
+ {
+ return velocityNormalOutside_[localSubFaceIdx];
+ }
private:
- Scalar velocity_;
+
+ Scalar velocitySelf_;
+ Scalar velocityOpposite_;
+ std::array velocityParallel_;
+ std::array velocityNormalInside_;
+ std::array velocityNormalOutside_;
};
} // end namespace
diff --git a/dumux/discretization/staggered/fvelementgeometry.hh b/dumux/discretization/staggered/fvelementgeometry.hh
index 1806d66e596e7dbd4a06c773229daf0cebf70624..ef34bf41c2da24b92ad7bb0891fc67ad7db5afb5 100644
--- a/dumux/discretization/staggered/fvelementgeometry.hh
+++ b/dumux/discretization/staggered/fvelementgeometry.hh
@@ -28,7 +28,6 @@
#include
#include
-#include
namespace Dumux
{
@@ -135,7 +134,7 @@ public:
void bindElement(const Element& element)
{
elementPtr_ = &element;
- scvIndices_ = std::vector({fvGridGeometry().problem_().elementMapper().index(*elementPtr_)});
+ scvIndices_ = std::vector({fvGridGeometry().elementMapper().index(*elementPtr_)});
}
//! The global finite volume geometry we are a restriction of
diff --git a/dumux/discretization/staggered/fvgridgeometry.hh b/dumux/discretization/staggered/fvgridgeometry.hh
index 8fd8b53b95d5c96c2d172ca38b4420143fa82da5..5dce791e7350f3791322887ed8071b0b15658762 100644
--- a/dumux/discretization/staggered/fvgridgeometry.hh
+++ b/dumux/discretization/staggered/fvgridgeometry.hh
@@ -26,7 +26,9 @@
#define DUMUX_DISCRETIZATION_STAGGERED_GLOBAL_FVGEOMETRY_HH
#include
-#include
+#include
+#include
+#include
namespace Dumux
{
@@ -43,9 +45,9 @@ class StaggeredFVGridGeometry
// specialization in case the FVElementGeometries are stored globally
template
-class StaggeredFVGridGeometry
+class StaggeredFVGridGeometry : public BaseFVGridGeometry
{
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using ParentType = BaseFVGridGeometry;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
@@ -64,11 +66,12 @@ class StaggeredFVGridGeometry
};
using GeometryHelper = typename GET_PROP_TYPE(TypeTag, StaggeredGeometryHelper);
+ using ConnectivityMap = StaggeredConnectivityMap;
public:
//! Constructor
StaggeredFVGridGeometry(const GridView& gridView)
- : gridView_(gridView), elementMap_(gridView), intersectionMapper_(gridView) {}
+ : ParentType(gridView), elementMap_(gridView), intersectionMapper_(gridView) {}
//! The total number of sub control volumes
std::size_t numScv() const
@@ -94,6 +97,16 @@ public:
return intersectionMapper_.numIntersections();
}
+ //! the total number of dofs
+ std::size_t numDofs() const
+ { return numCellCenterDofs() + numFaceDofs(); }
+
+ std::size_t numCellCenterDofs() const
+ { return this->gridView().size(0); }
+
+ std::size_t numFaceDofs() const
+ { return this->gridView().size(1); }
+
// Get an element from a sub control volume contained in it
Element element(const SubControlVolume& scv) const
{ return elementMap_.element(scv.elementIndex()); }
@@ -102,15 +115,9 @@ public:
Element element(IndexType eIdx) const
{ return elementMap_.element(eIdx); }
- //! Return the gridView this global object lives on
- const GridView& gridView() const
- { return gridView_; }
-
//! update all fvElementGeometries (do this again after grid adaption)
- void update(const Problem& problem)
+ void update()
{
- problemPtr_ = &problem;
-
// clear containers (necessary after grid refinement)
scvs_.clear();
scvfs_.clear();
@@ -119,9 +126,9 @@ public:
intersectionMapper_.update();
// determine size of containers
- IndexType numScvs = gridView_.size(0);
+ IndexType numScvs = this->gridView().size(0);
IndexType numScvf = 0;
- for (const auto& element : elements(gridView_))
+ for (const auto& element : elements(this->gridView()))
numScvf += element.subEntities(1);
// reserve memory
@@ -134,9 +141,9 @@ public:
// Build the scvs and scv faces
IndexType scvfIdx = 0;
numBoundaryScvf_ = 0;
- for (const auto& element : elements(gridView_))
+ for (const auto& element : elements(this->gridView()))
{
- auto eIdx = problem.elementMapper().index(element);
+ auto eIdx = this->elementMapper().index(element);
// reserve memory for the localToGlobalScvfIdx map
auto numLocalFaces = intersectionMapper_.numFaces(element);
@@ -151,9 +158,9 @@ public:
std::vector scvfsIndexSet;
scvfsIndexSet.reserve(numLocalFaces);
- GeometryHelper geometryHelper(element, gridView_);
+ GeometryHelper geometryHelper(element, this->gridView());
- for (const auto& intersection : intersections(gridView_, element))
+ for (const auto& intersection : intersections(this->gridView(), element))
{
geometryHelper.updateLocalFace(intersectionMapper_, intersection);
const int localFaceIndex = geometryHelper.localFaceIndex();
@@ -161,7 +168,7 @@ public:
// inner sub control volume faces
if (intersection.neighbor())
{
- auto nIdx = problem.elementMapper().index(intersection.outside());
+ auto nIdx = this->elementMapper().index(intersection.outside());
scvfs_.emplace_back(intersection,
intersection.geometry(),
scvfIdx,
@@ -177,7 +184,7 @@ public:
scvfs_.emplace_back(intersection,
intersection.geometry(),
scvfIdx,
- std::vector({eIdx, gridView_.size(0) + numBoundaryScvf_++}),
+ std::vector({eIdx, this->gridView().size(0) + numBoundaryScvf_++}),
geometryHelper
);
localToGlobalScvfIndices_[eIdx][localFaceIndex] = scvfIdx;
@@ -188,15 +195,10 @@ public:
// Save the scvf indices belonging to this scv to build up fv element geometries fast
scvfIndicesOfScv_[eIdx] = scvfsIndexSet;
}
- }
- /*!
- * \brief Return a local restriction of this global object
- * The local object is only functional after calling its bind/bindElement method
- * This is a free function that will be found by means of ADL
- */
- friend inline FVElementGeometry localView(const StaggeredFVGridGeometry& global)
- { return FVElementGeometry(global); }
+ // build the connectivity map for an effecient assembly
+ connectivityMap_.update(*this);
+ }
//private:
@@ -228,16 +230,21 @@ public:
return scvf(localToGlobalScvfIndex(eIdx, localScvfIdx));
}
+ /*!
+ * \brief Returns the connectivity map of which dofs have derivatives with respect
+ * to a given dof.
+ */
+ const ConnectivityMap &connectivityMap() const
+ { return connectivityMap_; }
-private:
- const Problem& problem_() const
- { return *problemPtr_; }
- const Problem* problemPtr_;
+private:
- const GridView gridView_;
+ // mappers
+ ConnectivityMap connectivityMap_;
Dumux::ElementMap elementMap_;
IntersectionMapper intersectionMapper_;
+
std::vector scvs_;
std::vector scvfs_;
std::vector> scvfIndicesOfScv_;
@@ -249,132 +256,7 @@ private:
template
class StaggeredFVGridGeometry
{
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using IndexType = typename GridView::IndexSet::IndexType;
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using Element = typename GridView::template Codim<0>::Entity;
- //! The local fvGeometry needs access to the problem
- friend typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
-
-public:
- //! Constructor
- StaggeredFVGridGeometry(const GridView& gridView)
- : gridView_(gridView), elementMap_(gridView) {}
-
- //! The total number of sub control volumes
- std::size_t numScv() const
- {
- return numScvs_;
- }
-
- //! The total number of sub control volume faces
- std::size_t numScvf() const
- {
- return numScvf_;
- }
-
- //! The total number of boundary sub control volume faces
- std::size_t numBoundaryScvf() const
- {
- return numBoundaryScvf_;
- }
-
- // Get an element from a sub control volume contained in it
- Element element(const SubControlVolume& scv) const
- { return elementMap_.element(scv.elementIndex()); }
-
- // Get an element from a global element index
- Element element(IndexType eIdx) const
- { return elementMap_.element(eIdx); }
-
- //! Return the gridView this global object lives on
- const GridView& gridView() const
- { return gridView_; }
-
- //! update all fvElementGeometries (do this again after grid adaption)
- void update(const Problem& problem)
- {
- problemPtr_ = &problem;
- elementMap_.clear();
-
- // reserve memory or resize the containers
- numScvs_ = gridView_.size(0);
- numScvf_ = 0;
- numBoundaryScvf_ = 0;
- elementMap_.resize(numScvs_);
- scvfIndicesOfScv_.resize(numScvs_);
- neighborVolVarIndices_.resize(numScvs_);
-
- // Build the SCV and SCV face
- for (const auto& element : elements(gridView_))
- {
- auto eIdx = problem.elementMapper().index(element);
-
- // fill the element map with seeds
- elementMap_[eIdx] = element.seed();
-
- // the element-wise index sets for finite volume geometry
- auto numLocalFaces = element.subEntities(1);
- std::vector scvfsIndexSet;
- std::vector neighborVolVarIndexSet;
- scvfsIndexSet.reserve(numLocalFaces);
- neighborVolVarIndexSet.reserve(numLocalFaces);
- for (const auto& intersection : intersections(gridView_, element))
- {
- // inner sub control volume faces
- if (intersection.neighbor())
- {
- scvfsIndexSet.push_back(numScvf_++);
- neighborVolVarIndexSet.push_back(problem.elementMapper().index(intersection.outside()));
- }
- // boundary sub control volume faces
- else if (intersection.boundary())
- {
- scvfsIndexSet.push_back(numScvf_++);
- neighborVolVarIndexSet.push_back(numScvs_ + numBoundaryScvf_++);
- }
- }
-
- // store the sets of indices in the data container
- scvfIndicesOfScv_[eIdx] = scvfsIndexSet;
- neighborVolVarIndices_[eIdx] = neighborVolVarIndexSet;
- }
- }
-
- const std::vector& scvfIndicesOfScv(IndexType scvIdx) const
- { return scvfIndicesOfScv_[scvIdx]; }
-
- const std::vector& neighborVolVarIndices(IndexType scvIdx) const
- { return neighborVolVarIndices_[scvIdx]; }
-
- /*!
- * \brief Return a local restriction of this global object
- * The local object is only functional after calling its bind/bindElement method
- * This is a free function that will be found by means of ADL
- */
- friend inline FVElementGeometry localView(const StaggeredFVGridGeometry& global)
- { return FVElementGeometry(global); }
-
-private:
- const Problem& problem_() const
- { return *problemPtr_; }
-
- const Problem* problemPtr_;
-
- const GridView gridView_;
-
- // Information on the global number of geometries
- IndexType numScvs_;
- IndexType numScvf_;
- IndexType numBoundaryScvf_;
-
- // vectors that store the global data
- Dumux::ElementMap elementMap_;
- std::vector> scvfIndicesOfScv_;
- std::vector> neighborVolVarIndices_;
+ // TODO: implement without caching
};
} // end namespace
diff --git a/dumux/discretization/staggered/globalfacevariables.hh b/dumux/discretization/staggered/globalfacevariables.hh
index 1e63de429179b8050d417206313404084c7870ed..f368f250da120fc90bad067e4be0eea2fba402e6 100644
--- a/dumux/discretization/staggered/globalfacevariables.hh
+++ b/dumux/discretization/staggered/globalfacevariables.hh
@@ -23,32 +23,54 @@
#ifndef DUMUX_DISCRETIZATION_STAGGERED_GLOBAL_FACEVARIABLES_HH
#define DUMUX_DISCRETIZATION_STAGGERED_GLOBAL_FACEVARIABLES_HH
-#include
+#include
+#include
namespace Dumux
{
-template
+namespace Properties
+{
+ NEW_PROP_TAG(ElementFaceVariables);
+}
+
+template
class StaggeredGlobalFaceVariables
+{};
+
+template
+class StaggeredGlobalFaceVariables
{
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 FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using FaceVariables = typename GET_PROP_TYPE(TypeTag, FaceVariables);
+ using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using IndexType = typename GridView::IndexSet::IndexType;
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
public:
- void update(Problem& problem, const FaceSolutionVector& sol)
- {
- problemPtr_ = &problem;
+ StaggeredGlobalFaceVariables(const Problem& problem) : problemPtr_(&problem) {}
- faceVariables_.resize(problem.model().numFaceDofs());
- assert(faceVariables_.size() == sol.size());
+ void update(const FVGridGeometry& fvGridGeometry, const SolutionVector& sol)
+ {
+ const auto& faceSol = sol[faceIdx];
+ faceVariables_.resize(fvGridGeometry.numScvf());
- for(int i = 0; i < problem.model().numFaceDofs(); ++i)
+ for(auto&& element : elements(fvGridGeometry.gridView()))
{
- faceVariables_[i].update(sol[i]);
+ auto fvGeometry = localView(fvGridGeometry);
+ fvGeometry.bindElement(element);
+
+ for(auto&& scvf : scvfs(fvGeometry))
+ {
+ faceVariables_[scvf.index()].update(faceSol, problem(), element, fvGeometry, scvf);
+ }
}
}
@@ -57,15 +79,65 @@ public:
FaceVariables& faceVars(const IndexType facetIdx)
{ return faceVariables_[facetIdx]; }
-private:
- const Problem& problem_() const
+
+
+ /*!
+ * \brief Return a local restriction of this global object
+ * The local object is only functional after calling its bind/bindElement method
+ * This is a free function that will be found by means of ADL
+ */
+ friend inline ElementFaceVariables localView(const StaggeredGlobalFaceVariables& global)
+ { return ElementFaceVariables(global); }
+
+ const Problem& problem() const
{ return *problemPtr_; }
+private:
+
const Problem* problemPtr_;
std::vector faceVariables_;
};
+template
+class StaggeredGlobalFaceVariables
+{
+ 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 SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using FaceVariables = typename GET_PROP_TYPE(TypeTag, FaceVariables);
+ using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using IndexType = typename GridView::IndexSet::IndexType;
+
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+public:
+ StaggeredGlobalFaceVariables(const Problem& problem) : problemPtr_(&problem) {}
+
+ void update(const FVGridGeometry& fvGridGeometry, const SolutionVector& sol)
+ { }
+
+ /*!
+ * \brief Return a local restriction of this global object
+ * The local object is only functional after calling its bind/bindElement method
+ * This is a free function that will be found by means of ADL
+ */
+ friend inline ElementFaceVariables localView(const StaggeredGlobalFaceVariables& global)
+ { return ElementFaceVariables(global); }
+
+ const Problem& problem() const
+ { return *problemPtr_; }
+
+
+private:
+
+ const Problem* problemPtr_;
+};
+
} // end namespace
diff --git a/dumux/discretization/staggered/globalfluxvariablescache.hh b/dumux/discretization/staggered/globalfluxvariablescache.hh
index 5e4dfd38cf687086c75fab3990c21e668ffdd51f..432b148255209f516fe418e4dc5eca446a5dbffb 100644
--- a/dumux/discretization/staggered/globalfluxvariablescache.hh
+++ b/dumux/discretization/staggered/globalfluxvariablescache.hh
@@ -23,7 +23,7 @@
#ifndef DUMUX_DISCRETIZATION_STAGGERED_GLOBAL_FLUXVARSCACHE_HH
#define DUMUX_DISCRETIZATION_STAGGERED_GLOBAL_FLUXVARSCACHE_HH
-#include
+#include
#include
namespace Dumux
@@ -47,32 +47,38 @@ class StaggeredGlobalFluxVariablesCache
friend StaggeredElementFluxVariablesCache;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GlobalVolumeVariables);
using IndexType = typename GridView::IndexSet::IndexType;
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
public:
+ StaggeredGlobalFluxVariablesCache(const Problem& problem) : problemPtr_(&problem) {}
+
// When global caching is enabled, precompute transmissibilities and stencils for all the scv faces
- void update(Problem& problem)
+ void update(const FVGridGeometry& fvGridGeometry,
+ const GridVolumeVariables& gridVolVars,
+ const SolutionVector& sol,
+ bool forceUpdate = false)
{
- problemPtr_ = &problem;
- const auto& fvGridGeometry = problem.model().fvGridGeometry();
- fluxVarsCache_.resize(fvGridGeometry.numScvf());
- for (const auto& element : elements(problem.gridView()))
- {
- // Prepare the geometries within the elements of the stencil
- auto fvGeometry = localView(fvGridGeometry);
- fvGeometry.bind(element);
-
- auto elemVolVars = localView(problem.model().curGlobalVolVars());
- elemVolVars.bind(element, fvGeometry, problem.model().curSol());
-
- for (auto&& scvf : scvfs(fvGeometry))
- {
- fluxVarsCache_[scvf.index()].update(problem, element, fvGeometry, elemVolVars, scvf);
- }
- }
+ // fluxVarsCache_.resize(fvGridGeometry.numScvf());
+ // for (const auto& element : elements(fvGridGeometry.gridView()))
+ // {
+ // // Prepare the geometries within the elements of the stencil
+ // auto fvGeometry = localView(fvGridGeometry);
+ // fvGeometry.bind(element);
+ //
+ // auto elemVolVars = localView(gridVolVars);
+ // elemVolVars.bind(element, fvGeometry, sol);
+ //
+ // for (auto&& scvf : scvfs(fvGeometry))
+ // {
+ // fluxVarsCache_[scvf.index()].update(problem, element, fvGeometry, elemVolVars, scvf);
+ // }
+ // }
}
/*!
@@ -83,6 +89,9 @@ public:
friend inline ElementFluxVariablesCache localView(const StaggeredGlobalFluxVariablesCache& global)
{ return ElementFluxVariablesCache(global); }
+ const Problem& problem() const
+ { return *problemPtr_; }
+
private:
// access operators in the case of caching
const FluxVariablesCache& operator [](IndexType scvfIdx) const
@@ -91,9 +100,6 @@ private:
FluxVariablesCache& operator [](IndexType scvfIdx)
{ return fluxVarsCache_[scvfIdx]; }
- const Problem& problem_() const
- { return *problemPtr_; }
-
const Problem* problemPtr_;
std::vector fluxVarsCache_;
diff --git a/dumux/discretization/staggered/globalvolumevariables.hh b/dumux/discretization/staggered/globalvolumevariables.hh
index 04477eb7bd04e9b3da4cd60026d604e2f077634e..3514b8c35e38dbb0566edc8c1a6dc38fbf56358b 100644
--- a/dumux/discretization/staggered/globalvolumevariables.hh
+++ b/dumux/discretization/staggered/globalvolumevariables.hh
@@ -23,9 +23,6 @@
#ifndef DUMUX_DISCRETIZATION_STAGGERED_GLOBAL_VOLUMEVARIABLES_HH
#define DUMUX_DISCRETIZATION_STAGGERED_GLOBAL_VOLUMEVARIABLES_HH
-#include
-#include
-
namespace Dumux
{
@@ -41,17 +38,10 @@ class StaggeredGlobalVolumeVariables
template
class StaggeredGlobalVolumeVariables
{
- // The local class needs to access and change volVars
- friend StaggeredElementVolumeVariables;
- // The local jacobian needs to access and change volVars for derivative calculation
- friend typename GET_PROP_TYPE(TypeTag, LocalJacobian);
- // as does the primary variable switch
- friend class PrimaryVariableSwitch;
- friend typename GET_PROP_TYPE(TypeTag, PrimaryVariableSwitch);
-
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
@@ -71,17 +61,17 @@ class StaggeredGlobalVolumeVariables
enum { numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter) };
public:
- void update(Problem& problem, const SolutionVector& sol)
- {
- problemPtr_ = &problem;
+ StaggeredGlobalVolumeVariables(const Problem& problem) : problemPtr_(&problem) {}
- auto numScv = problem.model().fvGridGeometry().numScv();
- auto numBoundaryScvf = problem.model().fvGridGeometry().numBoundaryScvf();
+ void update(const FVGridGeometry& fvGridGeometry, const SolutionVector& sol)
+ {
+ auto numScv = fvGridGeometry.numScv();
+ auto numBoundaryScvf = fvGridGeometry.numBoundaryScvf();
volumeVariables_.resize(numScv + numBoundaryScvf);
- for (const auto& element : elements(problem.gridView()))
+ for (const auto& element : elements(fvGridGeometry.gridView()))
{
- auto fvGeometry = localView(problem.model().fvGridGeometry());
+ auto fvGeometry = localView(fvGridGeometry);
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
@@ -89,7 +79,7 @@ public:
PrimaryVariables priVars(0.0);
priVars[cellCenterIdx] = sol[cellCenterIdx][scv.dofIndex()];
ElementSolutionVector elemSol{std::move(priVars)};
- volumeVariables_[scv.dofIndex()].update(elemSol, problem, element, scv);
+ volumeVariables_[scv.dofIndex()].update(elemSol, problem(), element, scv);
}
// handle the boundary volume variables
@@ -99,7 +89,7 @@ public:
if (!scvf.boundary())
continue;
- const auto bcTypes = problem.boundaryTypes(element, scvf);
+ const auto bcTypes = problem().boundaryTypes(element, scvf);
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideScv = fvGeometry.scv(insideScvIdx);
@@ -108,7 +98,7 @@ public:
for(int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
{
if(bcTypes.isDirichlet(eqIdx) || bcTypes.isDirichletCell(eqIdx))
- boundaryPriVars[cellCenterIdx][eqIdx] = problem.dirichlet(element, scvf)[cellCenterIdx][eqIdx];
+ boundaryPriVars[cellCenterIdx][eqIdx] = problem().dirichlet(element, scvf)[cellCenterIdx][eqIdx];
else if(bcTypes.isNeumann(eqIdx) || bcTypes.isOutflow(eqIdx) || bcTypes.isSymmetry())
boundaryPriVars[cellCenterIdx][eqIdx] = sol[cellCenterIdx][scvf.insideScvIdx()][eqIdx];
//TODO: this assumes a zero-gradient for e.g. the pressure on the boundary
@@ -118,7 +108,7 @@ public:
DUNE_THROW(Dune::InvalidStateException, "Face at: " << scvf.center() << " has neither Dirichlet nor Neumann BC.");
}
ElementSolutionVector elemSol{std::move(boundaryPriVars)};
- volumeVariables_[scvf.outsideScvIdx()].update(elemSol, problem, element, insideScv);
+ volumeVariables_[scvf.outsideScvIdx()].update(elemSol, problem(), element, insideScv);
}
}
}
@@ -143,10 +133,11 @@ public:
VolumeVariables& volVars(const SubControlVolume scv)
{ return volumeVariables_[scv.dofIndex()]; }
-private:
- const Problem& problem_() const
+ const Problem& problem() const
{ return *problemPtr_; }
+private:
+
const Problem* problemPtr_;
std::vector volumeVariables_;
@@ -157,15 +148,15 @@ private:
template
class StaggeredGlobalVolumeVariables
{
- // local class needs access to the problem
- friend StaggeredElementVolumeVariables;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
public:
- void update(Problem& problem, const SolutionVector& sol)
- { problemPtr_ = &problem; }
+ StaggeredGlobalVolumeVariables(const Problem& problem) : problemPtr_(&problem) {}
+
+ void update(const FVGridGeometry& fvGridGeometry, const SolutionVector& sol) {}
/*!
* \brief Return a local restriction of this global object
@@ -175,11 +166,12 @@ public:
friend inline ElementVolumeVariables localView(const StaggeredGlobalVolumeVariables& global)
{ return ElementVolumeVariables(global); }
-private:
- Problem& problem_() const
+ const Problem& problem() const
{ return *problemPtr_;}
- Problem* problemPtr_;
+private:
+
+ const Problem* problemPtr_;
};
} // end namespace
diff --git a/dumux/implicit/staggered/propertydefaults.hh b/dumux/discretization/staggered/properties.hh
similarity index 70%
rename from dumux/implicit/staggered/propertydefaults.hh
rename to dumux/discretization/staggered/properties.hh
index b1dc5044a9c3fc1800217f696484d76b2e3b794e..c54fbf52041ec073d39c76d9c2c5c78f2a734be0 100644
--- a/dumux/implicit/staggered/propertydefaults.hh
+++ b/dumux/discretization/staggered/properties.hh
@@ -18,69 +18,69 @@
*****************************************************************************/
/*!
* \ingroup Properties
- * \ingroup CCTpfaProperties
- * \ingroup StaggeredModel
* \file
*
- * \brief Default properties for cell centered models
+ * \brief Defines a type tag and some properties for models using the staggered scheme.
*/
-#ifndef DUMUX_STAGGERED_PROPERTY_DEFAULTS_HH
-#define DUMUX_STAGGERED_PROPERTY_DEFAULTS_HH
-#include
-#include
-#include
-#include
+#ifndef DUMUX_STAGGERDs_PROPERTIES_HH
+#define DUMUX_STAGGERDs_PROPERTIES_HH
+
+#include
+
#include
+#include
+#include
+#include
+#include
+#include
+
+#include
#include
#include
-
-#include
#include
+#include
+#include
+#include
#include
+#include
+#include
-#include
-#include
-#include
-#include
+#include
#include
#include
-#include
-
-#include
-#include
-
-#include "assembler.hh"
-#include "localresidual.hh"
-#include "localjacobian.hh"
-#include "properties.hh"
-#include "newtoncontroller.hh"
-#include "newtonconvergencewriter.hh"
-#include "model.hh"
-#include "primaryvariables.hh"
+#include
-namespace Dumux {
+namespace Dumux
+{
// forward declarations
template class CCElementBoundaryTypes;
-namespace Properties {
+namespace Properties
+{
+
+NEW_PROP_TAG(CellCenterSolutionVector);
+NEW_PROP_TAG(FaceSolutionVector);
+NEW_PROP_TAG(StaggeredFaceSolution);
+NEW_PROP_TAG(ElementFaceVariables);
+NEW_PROP_TAG(EnableGlobalFaceVariablesCache);
+
+//! Type tag for the box scheme.
+NEW_TYPE_TAG(StaggeredModel, INHERITS_FROM(FiniteVolumeModel, NumericModel, LinearSolverTypeTag));
+
//! Set the corresponding discretization method property
SET_PROP(StaggeredModel, DiscretizationMethod)
{
static const DiscretizationMethods value = DiscretizationMethods::Staggered;
};
-
-SET_TYPE_PROP(StaggeredModel, BaseModel, StaggeredBaseModel);
-
-
-//! Set the default for the global finite volume geometry
+//! Set the default for the FVElementGeometry vector
SET_TYPE_PROP(StaggeredModel, FVGridGeometry, StaggeredFVGridGeometry);
-//! Set the default for the local finite volume geometry
+//! Set the default for the FVElementGeometry vector
SET_TYPE_PROP(StaggeredModel, FVElementGeometry, StaggeredFVElementGeometry);
//! The sub control volume
@@ -94,53 +94,39 @@ public:
typedef Dumux::CCSubControlVolume type;
};
-SET_TYPE_PROP(StaggeredModel, GlobalFaceVars, Dumux::StaggeredGlobalFaceVariables);
+SET_TYPE_PROP(StaggeredModel, GlobalFaceVars, Dumux::StaggeredGlobalFaceVariables);
//! Set the default for the ElementBoundaryTypes
SET_TYPE_PROP(StaggeredModel, ElementBoundaryTypes, Dumux::CCElementBoundaryTypes);
-//! Mapper for the degrees of freedoms.
-SET_TYPE_PROP(StaggeredModel, DofMapper, typename GET_PROP_TYPE(TypeTag, ElementMapper));
+//! The global volume variables vector class
+SET_TYPE_PROP(StaggeredModel, GlobalVolumeVariables, StaggeredGlobalVolumeVariables);
-//! The global current volume variables vector class
-SET_TYPE_PROP(StaggeredModel, GlobalVolumeVariables, Dumux::StaggeredGlobalVolumeVariables);
+//! The element volume variables vector class
+SET_TYPE_PROP(StaggeredModel, ElementVolumeVariables, StaggeredElementVolumeVariables);
//! The global flux variables cache vector class
-SET_TYPE_PROP(StaggeredModel, GlobalFluxVariablesCache, Dumux::StaggeredGlobalFluxVariablesCache);
-
-//! The local jacobian operator
-SET_TYPE_PROP(StaggeredModel, LocalJacobian, Dumux::StaggeredLocalJacobian);
-
-//! Assembler for the global jacobian matrix
-SET_TYPE_PROP(StaggeredModel, JacobianAssembler, Dumux::StaggeredAssembler);
+SET_TYPE_PROP(StaggeredModel, GlobalFluxVariablesCache, StaggeredGlobalFluxVariablesCache);
//! The local flux variables cache vector class
-SET_TYPE_PROP(StaggeredModel, ElementFluxVariablesCache, Dumux::StaggeredElementFluxVariablesCache);
-
-//! The global previous volume variables vector class
-SET_TYPE_PROP(StaggeredModel, ElementVolumeVariables, Dumux::StaggeredElementVolumeVariables);
+SET_TYPE_PROP(StaggeredModel, ElementFluxVariablesCache, StaggeredElementFluxVariablesCache);
//! Set the BaseLocalResidual to StaggeredLocalResidual
-SET_TYPE_PROP(StaggeredModel, BaseLocalResidual, Dumux::StaggeredLocalResidual);
+SET_TYPE_PROP(StaggeredModel, BaseLocalResidual, StaggeredLocalResidual);
-//! Set the BaseLocalResidual to StaggeredLocalResidual
SET_TYPE_PROP(StaggeredModel, IntersectionMapper, Dumux::ConformingGridIntersectionMapper);
-//! indicate that this is no box discretization
-SET_BOOL_PROP(StaggeredModel, ImplicitIsBox, false);
+SET_TYPE_PROP(StaggeredModel, StaggeredFaceSolution, StaggeredFaceSolution);
-SET_TYPE_PROP(StaggeredModel, NewtonController, StaggeredNewtonController);
+SET_TYPE_PROP(StaggeredModel, ElementFaceVariables, StaggeredElementFaceVariables);
-SET_INT_PROP(StaggeredModel, NumEqCellCenter, 1);
-SET_INT_PROP(StaggeredModel, NumEqFace, 1);
+SET_BOOL_PROP(StaggeredModel, EnableGlobalFaceVariablesCache, true);
-SET_PROP(StaggeredModel, NumEq)
+//! Definition of the indices for cell center and face dofs in the global solution vector
+SET_PROP(StaggeredModel, DofTypeIndices)
{
-private:
- static constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
- static constexpr auto numEqFace = GET_PROP_VALUE(TypeTag, NumEqFace);
-public:
- static constexpr auto value = numEqCellCenter + numEqFace;
+ using CellCenterIdx = Dune::index_constant<0>;
+ using FaceIdx = Dune::index_constant<1>;
};
//! A vector of primary variables
@@ -208,19 +194,23 @@ public:
using type = typename Dune::MultiTypeBlockMatrix;
};
-//! Definition of the indices for cell center and face dofs in the global solution vector
-SET_PROP(StaggeredModel, DofTypeIndices)
+SET_PROP(StaggeredModel, NumEq)
{
- using CellCenterIdx = Dune::index_constant<0>;
- using FaceIdx = Dune::index_constant<1>;
+private:
+ static constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
+ static constexpr auto numEqFace = GET_PROP_VALUE(TypeTag, NumEqFace);
+public:
+ static constexpr auto value = numEqCellCenter + numEqFace;
};
-//! set default solver
-// SET_TYPE_PROP(StaggeredModel, LinearSolver, Dumux::GSBiCGSTABBackend);
-SET_TYPE_PROP(StaggeredModel, LinearSolver, Dumux::UMFPackBackend);
+SET_PROP(StaggeredModel, LinearSolverBlockSize)
+{
+ // LinearSolverAcceptsMultiTypeMatrix::value
+ // TODO: make somehow dependend? or only relevant for direct solvers?
+public:
+ static constexpr auto value = 1;
+};
-//! set the block level to 2, suitable for e.g. the Dune::MultiTypeBlockMatrix
-SET_INT_PROP(StaggeredModel, LinearSolverPreconditionerBlockLevel, 2);
//! Boundary types at a single degree of freedom
SET_PROP(StaggeredModel, BoundaryTypes)
@@ -243,17 +233,7 @@ public:
using type = StaggeredPrimaryVariables;
};
-//! use the plain newton convergence writer by default
-SET_TYPE_PROP(StaggeredModel, NewtonConvergenceWriter, StaggeredNewtonConvergenceWriter);
-
-//! Write separate vtp files for face variables by default
-SET_BOOL_PROP(StaggeredModel, VtkWriteFaceData, true);
-
-//! For compatibility
-SET_BOOL_PROP(StaggeredModel, EnableInteriorBoundaries, false);
-
-//! Set staggered vtk output module
-SET_TYPE_PROP(StaggeredModel, VtkOutputModule, StaggeredVtkOutputModule);
+SET_TYPE_PROP(StaggeredModel, GridVariables, StaggeredGridVariables);
//! Set one or different base epsilons for the calculations of the localJacobian's derivatives
SET_PROP(StaggeredModel, BaseEpsilon)
@@ -273,8 +253,9 @@ public:
}
};
-} // namespace Properties
+
+} // namespace Properties
} // namespace Dumux
#endif
diff --git a/dumux/freeflow/properties.hh b/dumux/freeflow/properties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b959a7503797421f4baf2ea3c5279557ac1ffcd7
--- /dev/null
+++ b/dumux/freeflow/properties.hh
@@ -0,0 +1,88 @@
+// -*- 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 . *
+ *****************************************************************************/
+/*!
+ * \ingroup Properties
+ * \file
+ *
+ * \brief Defines a type tag and some properties for free flow models.
+ */
+
+#ifndef DUMUX_FREE_FLOW_PROPERTIES_HH
+#define DUMUX_FREE_FLOW_PROPERTIES_HH
+
+#include
+#include
+#include
+#include
+
+#include "./staggered/boundarytypes.hh"
+
+namespace Dumux
+{
+namespace Properties
+{
+//! Type tag for models involving flow in porous media
+NEW_TYPE_TAG(FreeFlow);
+
+SET_INT_PROP(FreeFlow, NumEqFace, 1); //!< set the number of equations to 1
+
+//! The sub-controlvolume face
+SET_PROP(FreeFlow, SubControlVolumeFace)
+{
+private:
+ using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
+ using ScvfGeometry = typename Grid::template Codim<1>::Geometry;
+ using IndexType = typename Grid::LeafGridView::IndexSet::IndexType;
+public:
+ typedef Dumux::StaggeredSubControlVolumeFace type;
+};
+
+//! The geometry helper required for the stencils, etc.
+SET_PROP(FreeFlow, StaggeredGeometryHelper)
+{
+private:
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+public:
+ using type = StaggeredGeometryHelper;
+};
+
+//! The variables living on the faces
+SET_TYPE_PROP(FreeFlow, FaceVariables, StaggeredFaceVariables);
+
+//! A container class used to specify values for boundary conditions
+SET_PROP(FreeFlow, BoundaryValues)
+{
+private:
+ using CellCenterBoundaryValues = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using FaceBoundaryValues = Dune::FieldVector;
+public:
+ using type = StaggeredPrimaryVariables;
+};
+
+//! Boundary types at a single degree of freedom
+SET_TYPE_PROP(FreeFlow,
+ BoundaryTypes,
+ StaggeredFreeFlowBoundaryTypes);
+
+} // namespace Properties
+} // namespace Dumux
+
+ #endif
diff --git a/dumux/freeflow/staggered/fluxvariables.hh b/dumux/freeflow/staggered/fluxvariables.hh
index 0d941b9df6fbdebd7bcee224c9b121bca47ed2ac..a85555cdc4279fbe4201b78ae96dbff8f178247b 100644
--- a/dumux/freeflow/staggered/fluxvariables.hh
+++ b/dumux/freeflow/staggered/fluxvariables.hh
@@ -23,7 +23,7 @@
#ifndef DUMUX_FREELOW_IMPLICIT_FLUXVARIABLES_HH
#define DUMUX_FREELOW_IMPLICIT_FLUXVARIABLES_HH
-#include
+#include
#include
namespace Dumux
@@ -66,7 +66,6 @@ class FreeFlowFluxVariablesImpl
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
@@ -74,6 +73,9 @@ class FreeFlowFluxVariablesImpl
using IndexType = typename GridView::IndexSet::IndexType;
using Stencil = std::vector;
+ using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
+ using FaceVariables = typename GET_PROP_TYPE(TypeTag, FaceVariables);
+
static constexpr bool navierStokes = GET_PROP_VALUE(TypeTag, EnableInertiaTerms);
enum {
@@ -100,7 +102,7 @@ public:
const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
- const GlobalFaceVars& globalFaceVars,
+ const ElementFaceVariables& elemFaceVars,
const SubControlVolumeFace &scvf,
const FluxVariablesCache& fluxVarsCache)
{
@@ -111,13 +113,13 @@ public:
const auto& outsideVolVars = scvf.boundary() ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
CellCenterPrimaryVariables flux(0.0);
- const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndex()).velocity();
+ const Scalar velocity = elemFaceVars[scvf].velocitySelf();
const bool insideIsUpstream = sign(scvf.outerNormalScalar()) == sign(velocity) ? true : false;
const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
const auto& downstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
- const Scalar upWindWeight = GET_PROP_VALUE(TypeTag, ImplicitUpwindWeight);
+ static const Scalar upWindWeight = getParamFromGroup(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Implicit.UpwindWeight");
flux = (upWindWeight * upstreamVolVars.density() +
(1.0 - upWindWeight) * downstreamVolVars.density()) * velocity;
@@ -128,7 +130,6 @@ public:
}
void computeCellCenterToCellCenterStencil(Stencil& stencil,
- const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
@@ -141,7 +142,6 @@ public:
}
void computeCellCenterToFaceStencil(Stencil& stencil,
- const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
@@ -150,7 +150,6 @@ public:
}
void computeFaceToCellCenterStencil(Stencil& stencil,
- const Problem& problem,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
{
@@ -167,7 +166,6 @@ public:
}
void computeFaceToFaceStencil(Stencil& stencil,
- const Problem& problem,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
{
@@ -194,19 +192,19 @@ public:
* \param scvf The sub control volume face
* \param fvGeometry The finite-volume geometry
* \param elemVolVars All volume variables for the element
- * \param globalFaceVars The face variables
+ * \param elementFaceVars The face variables
*/
FacePrimaryVariables computeNormalMomentumFlux(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
- const GlobalFaceVars& globalFaceVars)
+ const ElementFaceVariables& elementFaceVars)
{
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideVolVars = elemVolVars[insideScvIdx];
- const Scalar velocitySelf = globalFaceVars.faceVars(scvf.dofIndex()).velocity() ;
- const Scalar velocityOpposite = globalFaceVars.faceVars(scvf.dofIndexOpposingFace()).velocity();
+ const Scalar velocitySelf = elementFaceVars[scvf].velocitySelf() ;
+ const Scalar velocityOpposite = elementFaceVars[scvf].velocityOpposite();
FacePrimaryVariables normalFlux(0.0);
if(navierStokes)
@@ -246,31 +244,27 @@ public:
* \param scvf The sub control volume face
* \param fvGeometry The finite-volume geometry
* \param elemVolVars All volume variables for the element
- * \param globalFaceVars The face variables
+ * \param elementFaceVars The face variables
*/
FacePrimaryVariables computeTangetialMomentumFlux(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
- const GlobalFaceVars& globalFaceVars)
+ const ElementFaceVariables& elementFaceVars)
{
FacePrimaryVariables tangentialFlux(0.0);
-
- // convenience function to get the velocity on a face
- auto velocity = [&globalFaceVars](const int dofIdx)
- {
- return globalFaceVars.faceVars(dofIdx).velocity();
- };
+ auto& faceVars = elementFaceVars[scvf];
+ const int numSubFaces = scvf.pairData().size();
// account for all sub-faces
- for(const auto& subFaceData : scvf.pairData())
+ for(int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx)
{
const auto eIdx = scvf.insideScvIdx();
- const auto& normalFace = fvGeometry.scvf(eIdx, subFaceData.localNormalFaceIdx);
+ const auto& normalFace = fvGeometry.scvf(eIdx, scvf.pairData()[localSubFaceIdx].localNormalFaceIdx);
// Check if we have a symmetry boundary condition. If yes, the tangental part of the momentum flux can be neglected.
- if(subFaceData.outerParallelFaceDofIdx < 0)
+ if(scvf.pairData()[localSubFaceIdx].outerParallelFaceDofIdx < 0)
{
// lambda to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
auto makeGhostFace = [eIdx] (const GlobalPosition& pos)
@@ -279,58 +273,41 @@ public:
};
// use the ghost face to check if there is a symmetry boundary condition and skip any further steps if yes
- const auto bcTypes = problem.boundaryTypes(element, makeGhostFace(subFaceData.virtualOuterParallelFaceDofPos));
+ const auto bcTypes = problem.boundaryTypes(element, makeGhostFace(scvf.pairData()[localSubFaceIdx].virtualOuterParallelFaceDofPos));
if(bcTypes.isSymmetry())
continue;
}
// if there is no symmetry boundary condition, proceed to calculate the tangential momentum flux
if(navierStokes)
- tangentialFlux += computeAdvectivePartOfTangentialMomentumFlux_(problem, element, scvf, normalFace, subFaceData, elemVolVars, velocity);
+ tangentialFlux += computeAdvectivePartOfTangentialMomentumFlux_(problem, element, scvf, normalFace, elemVolVars, faceVars, localSubFaceIdx);
- tangentialFlux += computeDiffusivePartOfTangentialMomentumFlux_(problem, element, scvf, normalFace, subFaceData, elemVolVars, velocity);
+ tangentialFlux += computeDiffusivePartOfTangentialMomentumFlux_(problem, element, scvf, normalFace, elemVolVars, faceVars, localSubFaceIdx);
}
return tangentialFlux;
}
private:
- template
FacePrimaryVariables computeAdvectivePartOfTangentialMomentumFlux_(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf,
const SubControlVolumeFace& normalFace,
- const SubFaceData& subFaceData,
const ElementVolumeVariables& elemVolVars,
- VelocityHelper velocity)
+ const FaceVariables& faceVars,
+ const int localSubFaceIdx)
{
- const Scalar transportingVelocity = velocity(subFaceData.normalPair.first);
+ const Scalar transportingVelocity = faceVars.velocityNormalInside(localSubFaceIdx);
const auto insideScvIdx = normalFace.insideScvIdx();
const auto outsideScvIdx = normalFace.outsideScvIdx();
- // lambda to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
- // auto makeGhostFace = [insideScvIdx] (const GlobalPosition& pos)
- // {
- // return SubControlVolumeFace(pos, std::vector{insideScvIdx,insideScvIdx});
- // };
-
const bool innerElementIsUpstream = ( sign(normalFace.outerNormalScalar()) == sign(transportingVelocity) );
const auto& upVolVars = innerElementIsUpstream ? elemVolVars[insideScvIdx] : elemVolVars[outsideScvIdx];
- Scalar transportedVelocity(0.0);
-
- if(innerElementIsUpstream)
- transportedVelocity = velocity(scvf.dofIndex());
- else
- {
- const int outerDofIdx = subFaceData.outerParallelFaceDofIdx;
- if(outerDofIdx >= 0)
- transportedVelocity = velocity(outerDofIdx);
- else // this is the case when the outer parallal dof would lie outside the domain TODO: discuss which one is better
- // transportedVelocity = problem.dirichlet(makeGhostFace(subFaceData.virtualOuterParallelFaceDofPos))[faceIdx][scvf.directionIndex()];
- transportedVelocity = problem.dirichlet(element, scvf)[faceIdx][scvf.directionIndex()];
- }
+ const Scalar transportedVelocity = innerElementIsUpstream ?
+ faceVars.velocitySelf() :
+ faceVars.velocityParallel(localSubFaceIdx);
const Scalar momentum = upVolVars.density() * transportedVelocity;
const int sgn = sign(normalFace.outerNormalScalar());
@@ -338,63 +315,46 @@ private:
return transportingVelocity * momentum * sgn * normalFace.area() * 0.5;
}
- template
FacePrimaryVariables computeDiffusivePartOfTangentialMomentumFlux_(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf,
const SubControlVolumeFace& normalFace,
- const SubFaceData& subFaceData,
const ElementVolumeVariables& elemVolVars,
- VelocityHelper velocity)
+ const FaceVariables& faceVars,
+ const int localSubFaceIdx)
{
FacePrimaryVariables tangentialDiffusiveFlux(0.0);
- const auto normalDirIdx = normalFace.directionIndex();
const auto insideScvIdx = normalFace.insideScvIdx();
const auto outsideScvIdx = normalFace.outsideScvIdx();
const auto& insideVolVars = elemVolVars[insideScvIdx];
const auto& outsideVolVars = elemVolVars[outsideScvIdx];
- // lambda to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
- auto makeGhostFace = [insideScvIdx] (const GlobalPosition& pos)
- {
- return SubControlVolumeFace(pos, std::vector{insideScvIdx,insideScvIdx});
- };
-
// the averaged viscosity at the face normal to our face of interest (where we assemble the face residual)
const Scalar muAvg = (insideVolVars.viscosity() + outsideVolVars.viscosity()) * 0.5;
// the normal derivative
- const int innerNormalVelocityIdx = subFaceData.normalPair.first;
- const int outerNormalVelocityIdx = subFaceData.normalPair.second;
-
- const Scalar innerNormalVelocity = velocity(innerNormalVelocityIdx);
-
- const Scalar outerNormalVelocity = outerNormalVelocityIdx >= 0 ?
- velocity(outerNormalVelocityIdx) :
- problem.dirichlet(element, makeGhostFace(subFaceData.virtualOuterNormalFaceDofPos))[faceIdx][normalDirIdx];
+ const Scalar innerNormalVelocity = faceVars.velocityNormalInside(localSubFaceIdx);
+ const Scalar outerNormalVelocity = faceVars.velocityNormalOutside(localSubFaceIdx);
const Scalar normalDeltaV = scvf.normalInPosCoordDir() ?
(outerNormalVelocity - innerNormalVelocity) :
(innerNormalVelocity - outerNormalVelocity);
- const Scalar normalDerivative = normalDeltaV / subFaceData.normalDistance;
+ const Scalar normalDerivative = normalDeltaV / scvf.pairData(localSubFaceIdx).normalDistance;
tangentialDiffusiveFlux -= muAvg * normalDerivative;
// the parallel derivative
- const Scalar innerParallelVelocity = velocity(scvf.dofIndex());
+ const Scalar innerParallelVelocity = faceVars.velocitySelf();
- const int outerParallelFaceDofIdx = subFaceData.outerParallelFaceDofIdx;
- const Scalar outerParallelVelocity = outerParallelFaceDofIdx >= 0 ?
- velocity(outerParallelFaceDofIdx) :
- problem.dirichlet(element, makeGhostFace(subFaceData.virtualOuterParallelFaceDofPos))[faceIdx][scvf.directionIndex()];
+ const Scalar outerParallelVelocity = faceVars.velocityParallel(localSubFaceIdx);
const Scalar parallelDeltaV = normalFace.normalInPosCoordDir() ?
(outerParallelVelocity - innerParallelVelocity) :
(innerParallelVelocity - outerParallelVelocity);
- const Scalar parallelDerivative = parallelDeltaV / subFaceData.parallelDistance;
+ const Scalar parallelDerivative = parallelDeltaV / scvf.pairData(localSubFaceIdx).parallelDistance;
tangentialDiffusiveFlux -= muAvg * parallelDerivative;
const Scalar sgn = sign(normalFace.outerNormalScalar());
diff --git a/dumux/freeflow/staggered/fluxvariablescache.hh b/dumux/freeflow/staggered/fluxvariablescache.hh
index e84095a7f12495816c0d8cef5de0836362b87b4c..10a568425b87c1889ae6b5617eca65425824af3c 100644
--- a/dumux/freeflow/staggered/fluxvariablescache.hh
+++ b/dumux/freeflow/staggered/fluxvariablescache.hh
@@ -23,7 +23,7 @@
#ifndef DUMUX_FREEFLOW_IMPLICIT_FLUXVARIABLESCACHE_HH
#define DUMUX_FREEFLOW_IMPLICIT_FLUXVARIABLESCACHE_HH
-#include
+#include
#include
#include