From 25e61c9fd16b565ea17256f4b834d79555d81273 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Fri, 27 Oct 2017 17:35:01 +0200
Subject: [PATCH] [staggered] First attempt to port to new structure of next
---
dumux/assembly/staggeredfvassembler.hh | 451 ++++++
dumux/assembly/staggeredlocalassembler.hh | 1397 +++++++++++++++++
dumux/common/staggeredfvproblem.hh | 159 ++
.../staggered/connectivitymap.hh} | 32 +-
.../staggered/fvelementgeometry.hh | 2 +-
.../staggered/fvgridgeometry.hh | 187 +--
.../staggered/globalfacevariables.hh | 24 +-
.../staggered/globalfluxvariablescache.hh | 48 +-
.../staggered/globalvolumevariables.hh | 31 +-
dumux/freeflow/staggered/fluxvariables.hh | 4 -
dumux/freeflow/staggered/localresidual.hh | 4 +-
dumux/freeflow/staggered/problem.hh | 33 +-
dumux/freeflow/staggered/propertydefaults.hh | 3 +
dumux/freeflow/staggered/velocityoutput.hh | 19 +-
dumux/freeflow/staggered/vtkoutputfields.hh | 61 +
dumux/implicit/staggered/assembler.hh | 174 --
dumux/implicit/staggered/gridvariables.hh | 115 ++
dumux/implicit/staggered/localjacobian.hh | 594 -------
dumux/implicit/staggered/localresidual.hh | 232 +--
dumux/implicit/staggered/newtoncontroller.hh | 4 +-
dumux/implicit/staggered/propertydefaults.hh | 15 +-
dumux/io/pointcloudvtkwriter.hh | 2 +-
dumux/io/staggeredvtkoutputmodule.hh | 427 ++---
test/freeflow/staggered/doneatestproblem.hh | 48 +-
test/freeflow/staggered/test_donea.cc | 115 +-
test/freeflow/staggered/test_donea.input | 9 +-
26 files changed, 2825 insertions(+), 1365 deletions(-)
create mode 100644 dumux/assembly/staggeredfvassembler.hh
create mode 100644 dumux/assembly/staggeredlocalassembler.hh
create mode 100644 dumux/common/staggeredfvproblem.hh
rename dumux/{implicit/staggered/assemblymap.hh => discretization/staggered/connectivitymap.hh} (84%)
create mode 100644 dumux/freeflow/staggered/vtkoutputfields.hh
delete mode 100644 dumux/implicit/staggered/assembler.hh
create mode 100644 dumux/implicit/staggered/gridvariables.hh
delete mode 100644 dumux/implicit/staggered/localjacobian.hh
diff --git a/dumux/assembly/staggeredfvassembler.hh b/dumux/assembly/staggeredfvassembler.hh
new file mode 100644
index 0000000000..28020a6bed
--- /dev/null
+++ b/dumux/assembly/staggeredfvassembler.hh
@@ -0,0 +1,451 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief An assembler for the global linear system
+ * for fully implicit models and vertex-centered discretization schemes.
+ */
+#ifndef DUMUX_STAGGERED_FV_ASSEMBLER_HH
+#define DUMUX_STAGGERED_FV_ASSEMBLER_HH
+
+#include <type_traits>
+
+#include <dune/istl/matrixindexset.hh>
+
+#include <dumux/common/timeloop.hh>
+#include <dumux/implicit/staggered/properties.hh>
+#include <dumux/implicit/staggered/localresidual.hh>
+#include <dumux/discretization/methods.hh>
+
+#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 TypeTag, DiffMethod diffMethod, bool isImplicit = true>
+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<Scalar>;
+
+ 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<TypeTag, diffMethod, isImplicit>;
+
+ 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<const Problem> problem,
+ std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+ std::shared_ptr<GridVariables> 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<const Problem> problem,
+ std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+ std::shared_ptr<GridVariables> gridVariables,
+ std::shared_ptr<TimeLoop> 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_)
+ {
+ jacobian_ = std::make_shared<JacobianMatrix>();
+ jacobian_->setBuildMode(JacobianMatrix::random);
+ setJacobianPattern();
+ }
+
+ if(!residual_)
+ {
+ residual_ = std::make_shared<SolutionVector>();
+ 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::assemble(*this, *jacobian_, *residual_, 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");
+ }
+
+ /*!
+ * \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<JacobianMatrix>();
+ 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<SolutionVector>();
+ 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(numDofs());
+ 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<JacobianMatrix> A,
+ std::shared_ptr<SolutionVector> r)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Setting linear system with existing matrix not implemented yet.");
+ // jacobian_ = A;
+ // residual_ = r;
+ //
+ // // check and/or set the BCRS matrix's build mode
+ // 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<JacobianMatrix>();
+
+ // convenience references
+ CCToCCMatrixBlock& A11 = (*jacobian_)[cellCenterIdx][cellCenterIdx];
+ CCToFaceMatrixBlock& A12 = (*jacobian_)[cellCenterIdx][faceIdx];
+ FaceToFaceMatrixBlock& A21 = (*jacobian_)[faceIdx][cellCenterIdx];
+ FaceToCCMatrixBlock& A22 = (*jacobian_)[faceIdx][faceIdx];
+
+ A11.setBuildMode(CCToCCMatrixBlock::random);
+ A12.setBuildMode(CCToFaceMatrixBlock::random);
+ A21.setBuildMode(FaceToFaceMatrixBlock::random);
+ A22.setBuildMode(FaceToCCMatrixBlock::random);
+
+ residual_ = std::make_shared<SolutionVector>();
+
+ 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 = numCellCenterDofs();
+ const auto numDofsFace = numFaceDofs();
+
+ // convenience references
+ CCToCCMatrixBlock& A11 = (*jacobian_)[cellCenterIdx][cellCenterIdx];
+ CCToFaceMatrixBlock& A12 = (*jacobian_)[cellCenterIdx][faceIdx];
+ FaceToFaceMatrixBlock& A21 = (*jacobian_)[faceIdx][cellCenterIdx];
+ FaceToCCMatrixBlock& 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(numCellCenterDofs());
+ (*residual_)[faceIdx].resize(numFaceDofs());
+ }
+
+ //! cell-centered schemes have one dof per cell
+ std::size_t numDofs() const
+ { return numCellCenterDofs() + numFaceDofs(); }
+
+ std::size_t numCellCenterDofs() const
+ { return gridView().size(0); }
+
+ std::size_t numFaceDofs() const
+ { return gridView().size(1); }
+
+ 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<const Problem> problem_;
+
+ // the finite volume geometry of the grid
+ std::shared_ptr<const FVGridGeometry> fvGridGeometry_;
+
+ // the variables container for the grid
+ std::shared_ptr<GridVariables> gridVariables_;
+
+ // shared pointers to the jacobian matrix and residual
+ std::shared_ptr<JacobianMatrix> jacobian_;
+ std::shared_ptr<SolutionVector> 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 0000000000..272d4f2d76
--- /dev/null
+++ b/dumux/assembly/staggeredlocalassembler.hh
@@ -0,0 +1,1397 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief An assembler for the global 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 <dune/istl/matrixindexset.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dumux/implicit/properties.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+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 TypeTag,
+ DiffMethod DM = DiffMethod::numeric,
+ bool implicit = true>
+class StaggeredLocalAssembler;
+
+
+template<class TypeTag>
+class StaggeredLocalAssembler<TypeTag,
+ DiffMethod::numeric,
+ /*implicit=*/true>
+{
+ 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 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);
+
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+
+ 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<class Assembler>
+ static void assemble(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
+ const Element& element, const SolutionVector& curSol)
+ {
+ const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
+ res[globalI] = assemble_(assembler, jac, element, curSol);
+ }
+
+ /*!
+ * \brief Computes the derivatives with respect to the given element and adds them
+ * to the global matrix.
+ */
+ template<class Assembler>
+ static void assemble(Assembler& assembler, JacobianMatrix& jac,
+ const Element& element, const SolutionVector& curSol)
+ {
+ assemble_(assembler, jac, element, curSol);
+ }
+
+ /*!
+ * \brief Assemble the residual only
+ */
+ template<class Assembler>
+ static void assemble(Assembler& assembler, SolutionVector& res,
+ const Element& element, const SolutionVector& curSol)
+ {
+ const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
+ res[globalI] = assemble_(assembler, element, curSol);
+ }
+
+ /*!
+ * \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<Scalar>(std::numeric_limits<Scalar>::epsilon(), 1.0);
+ */
+ static const Scalar baseEps = 1e-10;
+ assert(std::numeric_limits<Scalar>::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);
+ }
+
+private:
+ /*!
+ * \brief Computes the residual
+ *
+ * \return The element residual at the current solution.
+ */
+ template<class Assembler>
+ static NumEqVector assemble_(Assembler& assembler,
+ const Element& element, const SolutionVector& curSol)
+ {
+ // is the actual element a ghost element?
+ const bool isGhost = (element.partitionType() == Dune::GhostEntity);
+ if (isGhost) return NumEqVector(0.0);
+
+ // 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);
+
+ const bool isStationary = localResidual.isStationary();
+ auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
+ if (!isStationary)
+ prevElemVolVars.bindElement(element, fvGeometry, localResidual.prevSol());
+
+ // for compatibility with box models
+ ElementBoundaryTypes elemBcTypes;
+
+ // the actual element's current residual
+ NumEqVector residual(0.0);
+ if (isStationary)
+ {
+ residual = localResidual.eval(problem,
+ element,
+ fvGeometry,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+ }
+ else
+ {
+ residual = localResidual.eval(problem,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+ }
+
+ return residual;
+ }
+
+ /*!
+ * \brief Computes the derivatives with respect to the given element and adds them
+ * to the global matrix.
+ *
+ * \return The element residual at the current solution.
+ */
+ template<class Assembler>
+ static NumEqVector assemble_(Assembler& assembler, JacobianMatrix& A,
+ const Element& element, const SolutionVector& curSol)
+ {
+ // get some references for convenience
+ const auto& problem = assembler.problem();
+ const auto& fvGridGeometry = assembler.fvGridGeometry();
+ const auto& connectivityMap = fvGridGeometry.connectivityMap();
+ 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);
+
+ const bool isStationary = localResidual.isStationary();
+ auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
+ if (!isStationary)
+ prevElemVolVars.bindElement(element, fvGeometry, localResidual.prevSol());
+
+ // the global dof of the actual element
+ const auto globalI = fvGridGeometry.elementMapper().index(element);
+
+ // check for boundaries on the element
+ // TODO Do we need them for cell-centered models?
+ ElementBoundaryTypes elemBcTypes;
+ elemBcTypes.update(problem, element, fvGeometry);
+
+ // is the actual element a ghost element?
+ const bool isGhost = (element.partitionType() == Dune::GhostEntity);
+
+ // the actual element's current residual
+ NumEqVector residual(0.0);
+ if (!isGhost)
+ {
+ if (isStationary)
+ {
+ residual = localResidual.eval(problem,
+ element,
+ fvGeometry,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+ }
+ else
+ {
+ residual = localResidual.eval(problem,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+ }
+ }
+
+
+ // TODO Do we really need this??????????
+ // this->model_().updatePVWeights(fvGeometry);
+
+ //////////////////////////////////////////////////////////////////////////////////////////////////
+ // //
+ // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
+ // neighboring elements we do so by computing the derivatives of the fluxes which depend on the //
+ // actual element. In the actual element we evaluate the derivative of the entire residual. //
+ // //
+ //////////////////////////////////////////////////////////////////////////////////////////////////
+
+ static const std::string group = GET_PROP_VALUE(TypeTag, ModelParameterGroup);
+ static const int numericDifferenceMethod = getParamFromGroup<int>(group, "Implicit.NumericDifferenceMethod");
+
+ // get stencil informations
+ const auto numNeighbors = connectivityMap[globalI].size();
+
+ // container to store the neighboring elements
+ std::vector<Element> neighborElements;
+ neighborElements.reserve(numNeighbors);
+
+ // get the elements in which we need to evaluate the fluxes
+ // and calculate these in the undeflected state
+ Dune::BlockVector<NumEqVector> origFlux(numNeighbors);
+ origFlux = 0.0;
+ unsigned int j = 0;
+ for (const auto& dataJ : connectivityMap[globalI])
+ {
+ neighborElements.emplace_back(fvGridGeometry.element(dataJ.globalJ));
+ for (const auto scvfIdx : dataJ.scvfsJ)
+ {
+ origFlux[j] += localResidual.evalFlux(problem,
+ neighborElements.back(),
+ fvGeometry,
+ curElemVolVars,
+ elemFluxVarsCache,
+ fvGeometry.scvf(scvfIdx));
+ }
+ // increment neighbor counter
+ ++j;
+ }
+
+ // reference to the element's scv (needed later) and corresponding vol vars
+ const auto& scv = fvGeometry.scv(globalI);
+ auto& curVolVars = getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
+
+ // save a copy of the original privars and vol vars in order
+ // to restore the original solution after deflection
+ const auto origPriVars = curSol[globalI];
+ const auto origVolVars = curVolVars;
+
+ // element solution container to be deflected
+ ElementSolutionVector elemSol({origPriVars});
+
+ // derivatives in the neighbors with repect to the current elements
+ Dune::BlockVector<NumEqVector> neighborDeriv(numNeighbors);
+ for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
+ {
+ // reset derivatives of element dof with respect to itself
+ // as well as neighbor derivatives
+ NumEqVector partialDeriv(0.0);
+ neighborDeriv = 0.0;
+
+ if (isGhost)
+ partialDeriv[pvIdx] = 1.0;
+
+ Scalar eps = numericEpsilon(curVolVars.priVar(pvIdx));
+ Scalar delta = 0;
+
+ if (numericDifferenceMethod >= 0)
+ {
+ // we are not using backward differences, i.e. we need to
+ // calculate f(x + \epsilon)
+
+ // deflect primary variables
+ elemSol[0][pvIdx] += eps;
+ delta += eps;
+
+ // update the volume variables and the flux var cache
+ curVolVars.update(elemSol, problem, element, scv);
+ if (enableGlobalFluxVarsCache)
+ gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
+ else
+ elemFluxVarsCache.update(element, fvGeometry, curElemVolVars);
+
+ // calculate the residual with the deflected primary variables
+ if (!isGhost)
+ {
+ if (isStationary)
+ {
+ partialDeriv = localResidual.eval(problem,
+ element,
+ fvGeometry,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+ }
+ else
+ {
+ partialDeriv = localResidual.eval(problem,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+ }
+ }
+
+ // calculate the fluxes in the neighbors with the deflected primary variables
+ for (std::size_t k = 0; k < numNeighbors; ++k)
+ for (auto scvfIdx : connectivityMap[globalI][k].scvfsJ)
+ {
+ neighborDeriv[k] += localResidual.evalFlux(problem,
+ neighborElements[k],
+ fvGeometry,
+ curElemVolVars,
+ elemFluxVarsCache,
+ fvGeometry.scvf(scvfIdx));
+ }
+ }
+ else
+ {
+ // we are using backward differences, i.e. we don't need
+ // to calculate f(x + \epsilon) and we can recycle the
+ // (already calculated) residual f(x)
+ if (!isGhost)
+ partialDeriv = residual;
+ neighborDeriv = origFlux;
+ }
+
+ if (numericDifferenceMethod <= 0)
+ {
+ // we are not using forward differences, i.e. we
+ // need to calculate f(x - \epsilon)
+
+ // deflect the primary variables
+ elemSol[0][pvIdx] -= delta + eps;
+ delta += eps;
+
+ // update the volume variables and the flux var cache
+ curVolVars.update(elemSol, problem, element, scv);
+ if (enableGlobalFluxVarsCache)
+ gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
+ else
+ elemFluxVarsCache.update(element, fvGeometry, curElemVolVars);
+
+ // calculate the residual with the deflected primary variables and subtract it
+ if (!isGhost)
+ {
+ if (isStationary)
+ {
+ partialDeriv -= localResidual.eval(problem,
+ element,
+ fvGeometry,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+ }
+ else
+ {
+ partialDeriv -= localResidual.eval(problem,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+ }
+ }
+
+ // calculate the fluxes into element with the deflected primary variables
+ for (std::size_t k = 0; k < numNeighbors; ++k)
+ for (auto scvfIdx : connectivityMap[globalI][k].scvfsJ)
+ {
+ neighborDeriv[k] += localResidual.evalFlux(problem,
+ neighborElements[k],
+ fvGeometry,
+ curElemVolVars,
+ elemFluxVarsCache,
+ fvGeometry.scvf(scvfIdx));
+ }
+ }
+ else
+ {
+ // we are using forward differences, i.e. we don't need to
+ // calculate f(x - \epsilon) and we can recycle the
+ // (already calculated) residual f(x)
+ if (!isGhost)
+ partialDeriv -= residual;
+ neighborDeriv -= origFlux;
+ }
+
+ // divide difference in residuals by the magnitude of the
+ // deflections between the two function evaluation
+ if (!isGhost)
+ partialDeriv /= delta;
+ neighborDeriv /= delta;
+
+ // restore the original state of the scv's volume variables
+ curVolVars = origVolVars;
+
+ // restore the current element solution
+ elemSol[0][pvIdx] = origPriVars[pvIdx];
+
+ // add the current partial derivatives to the global jacobian matrix
+ for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ {
+ // the diagonal entries
+ A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
+
+ // off-diagonal entries
+ j = 0;
+ for (const auto& dataJ : connectivityMap[globalI])
+ A[dataJ.globalJ][globalI][eqIdx][pvIdx] += neighborDeriv[j++][eqIdx];
+ }
+ }
+
+ //////////////////////////////////////////////////////////////////////////////////////////////
+ // //
+ // Calculate derivatives of the dofs in the element with respect to user-defined additional //
+ // dof dependencies. We do so by evaluating the change in the source term of the current //
+ // element with respect to the primary variables at the given additional dofs. //
+ // //
+ //////////////////////////////////////////////////////////////////////////////////////////////
+
+ // const auto& additionalDofDepedencies = problem.getAdditionalDofDependencies(globalI);
+ // if (!additionalDofDepedencies.empty() && !isGhost)
+ // {
+ // // compute the source in the undeflected state
+ // auto source = localResidual.computeSource(element, fvGeometry, curElemVolVars, scv);
+ // source *= -scv.volume()*curVolVarsI.extrusionFactor();
+
+ // // deflect solution at given dofs and recalculate the source
+ // for (auto globalJ : additionalDofDependencies)
+ // {
+ // const auto& scvJ = fvGeometry.scv(globalJ);
+ // auto& curVolVarsJ = curElemVolVars[scv];
+ // const auto& elementJ = fvGridGeometry.element(globalJ);
+
+ // // save a copy of the original privars and volvars
+ // // to restore original solution after deflection
+ // const auto origPriVars = curSol[globalJ];
+ // const auto origVolVarsJ = curVolVarsJ;
+
+ // // derivatives with repect to the additional DOF we depend on
+ // for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
+ // {
+ // // derivatives of element dof with respect to itself
+ // NumEqVector partialDeriv(0.0);
+ // const auto eps = numericEpsilon(curVolVarsJ.priVar(pvIdx));
+ // Scalar delta = 0;
+
+ // if (numericDifferenceMethod >= 0)
+ // {
+ // // we are not using backward differences, i.e. we need to
+ // // calculate f(x + \epsilon)
+
+ // // deflect primary variables
+ // curSol[globalJ][pvIdx] += eps;
+ // delta += eps;
+
+ // // update the volume variables and the flux var cache
+ // curVolVarsJ.update(gridVariables.elementSolution(elementJ, curSol), problem, elementJ, scvJ);
+
+ // // calculate the source with the deflected primary variables
+ // auto deflSource = localResidual.computeSource(element, fvGeometry, curElemVolVars, scv);
+ // deflSource *= -scv.volume()*curVolVarsI.extrusionFactor();
+ // partialDeriv = std::move(deflSource);
+ // }
+ // else
+ // {
+ // // we are using backward differences, i.e. we don't need
+ // // to calculate f(x + \epsilon) and we can recycle the
+ // // (already calculated) source f(x)
+ // partialDeriv = source;
+ // }
+
+ // if (numericDifferenceMethod <= 0)
+ // {
+ // // we are not using forward differences, i.e. we
+ // // need to calculate f(x - \epsilon)
+
+ // // deflect the primary variables
+ // curSol[globalJ][pvIdx] -= delta + eps;
+ // delta += eps;
+
+ // // update the volume variables and the flux var cache
+ // curVolVarsJ.update(gridVariables.elementSolution(elementJ, curSol), problem, elementJ, scvJ);
+
+ // // calculate the source with the deflected primary variables and subtract
+ // auto deflSource = localResidual.computeSource(element, fvGeometry, curElemVolVars, scv);
+ // deflSource *= -scv.volume()*curVolVarsI.extrusionFactor();
+ // partialDeriv -= std::move(deflSource);
+ // }
+ // else
+ // {
+ // // we are using forward differences, i.e. we don't need to
+ // // calculate f(x - \epsilon) and we can recycle the
+ // // (already calculated) source f(x)
+ // partialDeriv -= source;
+ // }
+
+ // // divide difference in residuals by the magnitude of the
+ // // deflections between the two function evaluation
+ // partialDeriv /= delta;
+
+ // // restore the original state of the dofs privars and the volume variables
+ // curSol[globalJ] = origPriVars;
+ // curVolVarsJ = origVolVarsJ;
+
+ // // add the current partial derivatives to the global jacobian matrix
+ // for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ // A[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ // }
+ // }
+ // }
+
+ // return the original residual
+ return residual;
+ }
+private:
+ template<class T = TypeTag>
+ static typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), VolumeVariables&>::type
+ getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
+ { return elemVolVars[scv]; }
+
+ template<class T = TypeTag>
+ static typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), VolumeVariables&>::type
+ getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
+ { return gridVolVars.volVars(scv); }
+};
+
+//! Explicit assembler with numeric differentiation
+template<class TypeTag>
+class StaggeredLocalAssembler<TypeTag,
+ DiffMethod::numeric,
+ /*implicit=*/false>
+{
+ 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 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);
+
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+
+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<class Assembler>
+ static void assemble(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
+ const Element& element, const SolutionVector& curSol)
+ {
+ const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
+ res[globalI] = assemble_(assembler, jac, element, curSol);
+ }
+
+ /*!
+ * \brief Computes the derivatives with respect to the given element and adds them
+ * to the global matrix.
+ */
+ template<class Assembler>
+ static void assemble(Assembler& assembler, JacobianMatrix& jac,
+ const Element& element, const SolutionVector& curSol)
+ {
+ assemble_(assembler, jac, element, curSol);
+ }
+
+ /*!
+ * \brief Assemble the residual only
+ */
+ template<class Assembler>
+ static void assemble(Assembler& assembler, SolutionVector& res,
+ const Element& element, const SolutionVector& curSol)
+ {
+ const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
+ res[globalI] = assemble_(assembler, element, curSol);
+ }
+
+ /*!
+ * \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<Scalar>(std::numeric_limits<Scalar>::epsilon(), 1.0);
+ */
+ static const Scalar baseEps = 1e-10;
+ assert(std::numeric_limits<Scalar>::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);
+ }
+
+private:
+
+ /*!
+ * \brief Computes the residual
+ *
+ * \return The element residual at the current solution.
+ */
+ template<class Assembler>
+ static NumEqVector assemble_(Assembler& assembler,
+ const Element& element, const SolutionVector& curSol)
+ {
+ // is the actual element a ghost element?
+ const bool isGhost = (element.partitionType() == Dune::GhostEntity);
+ if (isGhost) return NumEqVector(0.0);
+
+ // get some references for convenience
+ const auto& problem = assembler.problem();
+ auto& localResidual = assembler.localResidual();
+ auto& gridVariables = assembler.gridVariables();
+
+ // using an explicit assembler doesn't make sense for stationary problems
+ if (localResidual.isStationary())
+ DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
+
+ // prepare the local views
+ auto fvGeometry = localView(assembler.fvGridGeometry());
+ fvGeometry.bind(element);
+
+ auto curElemVolVars = localView(gridVariables.curGridVolVars());
+ curElemVolVars.bindElement(element, fvGeometry, curSol);
+
+ auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
+ prevElemVolVars.bind(element, fvGeometry, localResidual.prevSol());
+
+ auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
+ elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
+
+ // compatibility with box method
+ ElementBoundaryTypes elemBcTypes;
+
+ // the actual element's previous time step residual
+ auto residual = localResidual.eval(problem,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+
+ auto storageResidual = localResidual.evalStorage(problem,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+
+ residual += storageResidual;
+
+ return residual;
+ }
+
+ /*!
+ * \brief Computes the derivatives with respect to the given element and adds them
+ * to the global matrix.
+ *
+ * \return The element residual at the current solution.
+ */
+ template<class Assembler>
+ static NumEqVector assemble_(Assembler& assembler, JacobianMatrix& A,
+ const Element& element, const SolutionVector& curSol)
+ {
+ // get some references for convenience
+ const auto& problem = assembler.problem();
+ const auto& fvGridGeometry = assembler.fvGridGeometry();
+ auto& localResidual = assembler.localResidual();
+ auto& gridVariables = assembler.gridVariables();
+
+ // using an explicit assembler doesn't make sense for stationary problems
+ if (localResidual.isStationary())
+ DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
+
+ // prepare the local views
+ auto fvGeometry = localView(assembler.fvGridGeometry());
+ fvGeometry.bind(element);
+
+ auto curElemVolVars = localView(gridVariables.curGridVolVars());
+ curElemVolVars.bindElement(element, fvGeometry, curSol);
+
+ auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
+ prevElemVolVars.bind(element, fvGeometry, localResidual.prevSol());
+
+ auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
+ elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
+
+ // the global dof of the actual element
+ const auto globalI = fvGridGeometry.elementMapper().index(element);
+
+ // check for boundaries on the element
+ // TODO Do we need them for cell-centered models?
+ ElementBoundaryTypes elemBcTypes;
+ elemBcTypes.update(problem, element, fvGeometry);
+
+ // is the actual element a ghost element?
+ const bool isGhost = (element.partitionType() == Dune::GhostEntity);
+
+ // the actual element's previous time step residual
+ NumEqVector residual(0.0), storageResidual(0.0);
+ if (!isGhost)
+ {
+ residual = localResidual.eval(problem,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+
+ storageResidual = localResidual.evalStorage(problem,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+
+ residual += storageResidual;
+ }
+
+ //////////////////////////////////////////////////////////////////////////////////////////////////
+ // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
+ // neighboring elements all derivatives are zero. For the assembled element only the storage //
+ // derivatives are non-zero. //
+ //////////////////////////////////////////////////////////////////////////////////////////////////
+
+ static const std::string group = GET_PROP_VALUE(TypeTag, ModelParameterGroup);
+ static const int numericDifferenceMethod = getParamFromGroup<int>(group, "Implicit.NumericDifferenceMethod");
+
+ // reference to the element's scv (needed later) and corresponding vol vars
+ const auto& scv = fvGeometry.scv(globalI);
+ auto& curVolVars = getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
+
+ // save a copy of the original privars and vol vars in order
+ // to restore the original solution after deflection
+ const auto origPriVars = curSol[globalI];
+ const auto origVolVars = curVolVars;
+
+ // element solution container to be deflected
+ ElementSolutionVector elemSol({origPriVars});
+
+ // derivatives in the neighbors with repect to the current elements
+ for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
+ {
+ // reset derivatives of element dof with respect to itself
+ // as well as neighbor derivatives
+ NumEqVector partialDeriv(0.0);
+
+ if (isGhost)
+ partialDeriv[pvIdx] = 1.0;
+
+ Scalar eps = numericEpsilon(curVolVars.priVar(pvIdx));
+ Scalar delta = 0;
+
+ if (numericDifferenceMethod >= 0)
+ {
+ // we are not using backward differences, i.e. we need to
+ // calculate f(x + \epsilon)
+
+ // deflect primary variables
+ elemSol[0][pvIdx] += eps;
+ delta += eps;
+
+ // update the volume variables and the flux var cache
+ curVolVars.update(elemSol, problem, element, scv);
+
+ // calculate the residual with the deflected primary variables
+ if (!isGhost)
+ {
+ partialDeriv = localResidual.evalStorage(problem,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+ }
+ }
+ else
+ {
+ // we are using backward differences, i.e. we don't need
+ // to calculate f(x + \epsilon) and we can recycle the
+ // (already calculated) residual f(x)
+ if (!isGhost)
+ partialDeriv = storageResidual;
+ }
+
+ if (numericDifferenceMethod <= 0)
+ {
+ // we are not using forward differences, i.e. we
+ // need to calculate f(x - \epsilon)
+
+ // deflect the primary variables
+ elemSol[0][pvIdx] -= delta + eps;
+ delta += eps;
+
+ // update the volume variables and the flux var cache
+ curVolVars.update(elemSol, problem, element, scv);
+
+ // calculate the residual with the deflected primary variables and subtract it
+ if (!isGhost)
+ {
+ partialDeriv -= localResidual.evalStorage(problem,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+ }
+ }
+ else
+ {
+ // we are using forward differences, i.e. we don't need to
+ // calculate f(x - \epsilon) and we can recycle the
+ // (already calculated) residual f(x)
+ if (!isGhost)
+ partialDeriv -= storageResidual;
+ }
+
+ // divide difference in residuals by the magnitude of the
+ // deflections between the two function evaluation
+ if (!isGhost)
+ partialDeriv /= delta;
+
+ // restore the original state of the scv's volume variables
+ curVolVars = origVolVars;
+
+ // restore the current element solution
+ elemSol[0][pvIdx] = origPriVars[pvIdx];
+
+ // add the current partial derivatives to the global jacobian matrix
+ for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ {
+ // the diagonal entries
+ A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ }
+ }
+
+ // return the original residual
+ return residual;
+ }
+private:
+ template<class T = TypeTag>
+ static typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), VolumeVariables&>::type
+ getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
+ { return elemVolVars[scv]; }
+
+ template<class T = TypeTag>
+ static typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), VolumeVariables&>::type
+ getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
+ { return gridVolVars.volVars(scv); }
+};
+
+//! implicit assembler using analytic differentiation
+template<class TypeTag>
+class StaggeredLocalAssembler<TypeTag,
+ DiffMethod::analytic,
+ /*implicit=*/true>
+{
+ 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 Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
+ using IndexType = typename GET_PROP_TYPE(TypeTag, GridView)::IndexSet::IndexType;
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
+ 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);
+
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+
+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<class Assembler>
+ static void assemble(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
+ const Element& element, const SolutionVector& curSol)
+ {
+ const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
+ res[globalI] = assemble_(assembler, jac, element, curSol);
+ }
+
+ /*!
+ * \brief Computes the derivatives with respect to the given element and adds them
+ * to the global matrix.
+ */
+ template<class Assembler>
+ static void assemble(Assembler& assembler, JacobianMatrix& jac,
+ const Element& element, const SolutionVector& curSol)
+ {
+ assemble_(assembler, jac, element, curSol);
+ }
+
+ /*!
+ * \brief Assemble the residual only
+ */
+ template<class Assembler>
+ static void assemble(Assembler& assembler, SolutionVector& res,
+ const Element& element, const SolutionVector& curSol)
+ {
+ const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
+ res[globalI] = assemble_(assembler, element, curSol);
+ }
+
+private:
+
+ /*!
+ * \brief Computes the residual
+ *
+ * \return The element residual at the current solution.
+ */
+ template<class Assembler>
+ static NumEqVector assemble_(Assembler& assembler,
+ const Element& element, const SolutionVector& curSol)
+ {
+ // is the actual element a ghost element?
+ const bool isGhost = (element.partitionType() == Dune::GhostEntity);
+ if (isGhost) return NumEqVector(0.0);
+
+ // 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 prevElemVolVars = localView(gridVariables.prevGridVolVars());
+
+ // check for boundaries on the element
+ // TODO Do we need them for cell-centered models?
+ ElementBoundaryTypes elemBcTypes;
+ elemBcTypes.update(problem, element, fvGeometry);
+
+ NumEqVector residual(0.0);
+ if (!localResidual.isStationary())
+ {
+ prevElemVolVars.bindElement(element, fvGeometry, localResidual.prevSol());
+
+ residual = localResidual.eval(problem,
+ element,
+ fvGeometry,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+ }
+ else
+ {
+ residual = localResidual.eval(problem,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+ }
+
+ return residual;
+ }
+
+ /*!
+ * \brief Computes the derivatives with respect to the given element and adds them
+ * to the global matrix.
+ *
+ * \return The element residual at the current solution.
+ */
+ template<class Assembler>
+ static NumEqVector assemble_(Assembler& assembler, JacobianMatrix& A,
+ const Element& element, const SolutionVector& curSol)
+ {
+ // get some references for convenience
+ const auto& problem = assembler.problem();
+ const auto& fvGridGeometry = assembler.fvGridGeometry();
+ 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);
+
+ const bool isStationary = localResidual.isStationary();
+ auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
+ if (!isStationary)
+ prevElemVolVars.bindElement(element, fvGeometry, localResidual.prevSol());
+
+ // the global dof of the actual element
+ const auto globalI = fvGridGeometry.elementMapper().index(element);
+
+ // check for boundaries on the element
+ // TODO Do we need them for cell-centered models?
+ ElementBoundaryTypes elemBcTypes;
+ elemBcTypes.update(problem, element, fvGeometry);
+
+ // is the actual element a ghost element?
+ const bool isGhost = (element.partitionType() == Dune::GhostEntity);
+
+ // the actual element's current residual (will be returned by this function)
+ NumEqVector residual(0.0);
+ if (!isGhost)
+ {
+ if (isStationary)
+ {
+ residual = localResidual.eval(problem,
+ element,
+ fvGeometry,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+ }
+ else
+ {
+ residual = localResidual.eval(problem,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+ }
+ }
+
+ // get reference to the element's current vol vars
+ const auto& scv = fvGeometry.scv(globalI);
+ const auto& volVars = curElemVolVars[scv];
+
+ // if the problem is instationary, add derivative of storage term
+ if (!isStationary)
+ localResidual.addStorageDerivatives(A[globalI][globalI],
+ problem,
+ element,
+ fvGeometry,
+ volVars,
+ scv);
+
+ // add source term derivatives
+ localResidual.addSourceDerivatives(A[globalI][globalI],
+ problem,
+ element,
+ fvGeometry,
+ volVars,
+ scv);
+
+ // add flux derivatives for each scvf
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+ if (!scvf.boundary())
+ {
+ localResidual.addFluxDerivatives(A[globalI],
+ problem,
+ element,
+ fvGeometry,
+ curElemVolVars,
+ elemFluxVarsCache,
+ scvf);
+ }
+ else
+ {
+ const auto& bcTypes = problem.boundaryTypes(element, scvf);
+
+ // add Dirichlet boundary flux derivatives
+ if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann())
+ {
+ localResidual.addCCDirichletFluxDerivatives(A[globalI],
+ problem,
+ element,
+ fvGeometry,
+ curElemVolVars,
+ elemFluxVarsCache,
+ scvf);
+ }
+ // add Robin ("solution dependent Neumann") boundary flux derivatives
+ else if (bcTypes.hasNeumann() && !bcTypes.hasDirichlet())
+ {
+ localResidual.addRobinFluxDerivatives(A[globalI],
+ problem,
+ element,
+ fvGeometry,
+ curElemVolVars,
+ elemFluxVarsCache,
+ scvf);
+ }
+ else
+ DUNE_THROW(Dune::NotImplemented, "Mixed boundary conditions. Use pure boundary conditions by converting Dirichlet BCs to Robin BCs");
+ }
+ }
+
+ // TODO Do we really need this??????????
+ // this->model_().updatePVWeights(fvGeometry);
+
+ // TODO: Additional dof dependencies???
+
+ // return element residual
+ return residual;
+ }
+};
+
+//! explicit assembler using analytic differentiation
+template<class TypeTag>
+class StaggeredLocalAssembler<TypeTag,
+ DiffMethod::analytic,
+ /*implicit=*/false>
+{
+ 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 Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
+ using IndexType = typename GET_PROP_TYPE(TypeTag, GridView)::IndexSet::IndexType;
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
+ 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);
+
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+
+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<class Assembler>
+ static void assemble(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
+ const Element& element, const SolutionVector& curSol)
+ {
+ const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
+ res[globalI] = assemble_(assembler, jac, element, curSol);
+ }
+
+ /*!
+ * \brief Computes the derivatives with respect to the given element and adds them
+ * to the global matrix.
+ */
+ template<class Assembler>
+ static void assemble(Assembler& assembler, JacobianMatrix& jac,
+ const Element& element, const SolutionVector& curSol)
+ {
+ assemble_(assembler, jac, element, curSol);
+ }
+
+ /*!
+ * \brief Assemble the residual only
+ */
+ template<class Assembler>
+ static void assemble(Assembler& assembler, SolutionVector& res,
+ const Element& element, const SolutionVector& curSol)
+ {
+ const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
+ res[globalI] = assemble_(assembler, element, curSol);
+ }
+
+private:
+
+ /*!
+ * \brief Computes the residual
+ *
+ * \return The element residual at the current solution.
+ */
+ template<class Assembler>
+ static NumEqVector assemble_(Assembler& assembler,
+ const Element& element, const SolutionVector& curSol)
+ {
+ // is the actual element a ghost element?
+ const bool isGhost = (element.partitionType() == Dune::GhostEntity);
+ if (isGhost) return NumEqVector(0.0);
+
+ // get some references for convenience
+ const auto& problem = assembler.problem();
+ auto& localResidual = assembler.localResidual();
+ auto& gridVariables = assembler.gridVariables();
+
+ // using an explicit assembler doesn't make sense for stationary problems
+ if (localResidual.isStationary())
+ DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
+
+ // prepare the local views
+ auto fvGeometry = localView(assembler.fvGridGeometry());
+ fvGeometry.bind(element);
+
+ auto curElemVolVars = localView(gridVariables.curGridVolVars());
+ curElemVolVars.bindElement(element, fvGeometry, curSol);
+
+ auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
+ prevElemVolVars.bind(element, fvGeometry, localResidual.prevSol());
+
+ auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
+ elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
+
+ // check for boundaries on the element
+ // TODO Do we need them for cell-centered models?
+ ElementBoundaryTypes elemBcTypes;
+ elemBcTypes.update(problem, element, fvGeometry);
+
+ // the actual element's previous time step residual
+ auto residual = localResidual.eval(problem,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+
+ auto storageResidual = localResidual.evalStorage(problem,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+
+ residual += storageResidual;
+
+ return residual;
+ }
+
+ /*!
+ * \brief Computes the derivatives with respect to the given element and adds them
+ * to the global matrix.
+ *
+ * \return The element residual at the current solution.
+ */
+ template<class Assembler>
+ static NumEqVector assemble_(Assembler& assembler, JacobianMatrix& A,
+ const Element& element, const SolutionVector& curSol)
+ {
+ // get some references for convenience
+ const auto& problem = assembler.problem();
+ const auto& fvGridGeometry = assembler.fvGridGeometry();
+ auto& localResidual = assembler.localResidual();
+ auto& gridVariables = assembler.gridVariables();
+
+ // using an explicit assembler doesn't make sense for stationary problems
+ if (localResidual.isStationary())
+ DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
+
+ // prepare the local views
+ auto fvGeometry = localView(assembler.fvGridGeometry());
+ fvGeometry.bind(element);
+
+ auto curElemVolVars = localView(gridVariables.curGridVolVars());
+ curElemVolVars.bindElement(element, fvGeometry, curSol);
+
+ auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
+ prevElemVolVars.bind(element, fvGeometry, localResidual.prevSol());
+
+ auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
+ elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
+
+ // the global dof of the actual element
+ const auto globalI = fvGridGeometry.elementMapper().index(element);
+
+ // check for boundaries on the element
+ // TODO Do we need them for cell-centered models?
+ ElementBoundaryTypes elemBcTypes;
+ elemBcTypes.update(problem, element, fvGeometry);
+
+ // is the actual element a ghost element?
+ const bool isGhost = (element.partitionType() == Dune::GhostEntity);
+
+ // the actual element's previous time step residual
+ NumEqVector residual(0.0), storageResidual(0.0);
+ if (!isGhost)
+ {
+ residual = localResidual.eval(problem,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+
+ storageResidual = localResidual.evalStorage(problem,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+
+ residual += storageResidual;
+ }
+
+ // get reference to the element's current vol vars
+ const auto& scv = fvGeometry.scv(globalI);
+ const auto& volVars = curElemVolVars[scv];
+
+ // add derivative of storage term
+ localResidual.addStorageDerivatives(A[globalI][globalI],
+ problem,
+ element,
+ fvGeometry,
+ volVars,
+ scv);
+
+ // return the original residual
+ return residual;
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/common/staggeredfvproblem.hh b/dumux/common/staggeredfvproblem.hh
new file mode 100644
index 0000000000..248c965b50
--- /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 <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for all problems
+ */
+#ifndef DUMUX_STAGGERD_FV_PROBLEM_HH
+#define DUMUX_STAGGERD_FV_PROBLEM_HH
+
+#include <dumux/implicit/staggered/properties.hh>
+#include <dumux/common/fvproblem.hh>
+
+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 TypeTag>
+class StaggeredFVProblem : public FVProblem<TypeTag>
+{
+ using ParentType = FVProblem<TypeTag>;
+ 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<dim>::Entity;
+ using CoordScalar = typename GridView::ctype;
+ using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
+
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+
+ // using GridAdaptModel = ImplicitGridAdapt<TypeTag, adaptiveGrid>;
+
+public:
+ /*!
+ * \brief Constructor
+ *
+ * \param gridView The simulation's idea about physical space
+ */
+ StaggeredFVProblem(std::shared_ptr<const FVGridGeometry> 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<class Entity>
+ 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<Implementation *>(this); }
+
+ //! \copydoc asImp_()
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation *>(this); }
+
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/implicit/staggered/assemblymap.hh b/dumux/discretization/staggered/connectivitymap.hh
similarity index 84%
rename from dumux/implicit/staggered/assemblymap.hh
rename to dumux/discretization/staggered/connectivitymap.hh
index 45331e432b..9863d8eb49 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
-
-#include <dune/istl/bcrsmatrix.hh>
+#ifndef DUMUX_STAGGERED_CONNECTIVITY_MAP_HH
+#define DUMUX_STAGGERED_CONNECTIVITY_MAP_HH
+#include <vector>
#include <dumux/implicit/properties.hh>
namespace Dumux
{
template<class TypeTag>
-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/fvelementgeometry.hh b/dumux/discretization/staggered/fvelementgeometry.hh
index 1806d66e59..8dca20b69e 100644
--- a/dumux/discretization/staggered/fvelementgeometry.hh
+++ b/dumux/discretization/staggered/fvelementgeometry.hh
@@ -135,7 +135,7 @@ public:
void bindElement(const Element& element)
{
elementPtr_ = &element;
- scvIndices_ = std::vector<IndexType>({fvGridGeometry().problem_().elementMapper().index(*elementPtr_)});
+ scvIndices_ = std::vector<IndexType>({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 8fd8b53b95..31883f8f02 100644
--- a/dumux/discretization/staggered/fvgridgeometry.hh
+++ b/dumux/discretization/staggered/fvgridgeometry.hh
@@ -27,6 +27,9 @@
#include <dumux/common/elementmap.hh>
#include <dumux/implicit/staggered/properties.hh>
+#include <dumux/discretization/basefvgridgeometry.hh>
+#include <dumux/common/boundingboxtree.hh>
+#include <dumux/discretization/staggered/connectivitymap.hh>
namespace Dumux
{
@@ -43,9 +46,9 @@ class StaggeredFVGridGeometry
// specialization in case the FVElementGeometries are stored globally
template<class TypeTag>
-class StaggeredFVGridGeometry<TypeTag, true>
+class StaggeredFVGridGeometry<TypeTag, true> : public BaseFVGridGeometry<TypeTag>
{
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using ParentType = BaseFVGridGeometry<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
@@ -64,11 +67,12 @@ class StaggeredFVGridGeometry<TypeTag, true>
};
using GeometryHelper = typename GET_PROP_TYPE(TypeTag, StaggeredGeometryHelper);
+ using ConnectivityMap = StaggeredConnectivityMap<TypeTag>;
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
@@ -102,15 +106,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 +117,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 +132,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 +149,9 @@ public:
std::vector<IndexType> 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 +159,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 +175,7 @@ public:
scvfs_.emplace_back(intersection,
intersection.geometry(),
scvfIdx,
- std::vector<IndexType>({eIdx, gridView_.size(0) + numBoundaryScvf_++}),
+ std::vector<IndexType>({eIdx, this->gridView().size(0) + numBoundaryScvf_++}),
geometryHelper
);
localToGlobalScvfIndices_[eIdx][localFaceIndex] = scvfIdx;
@@ -188,15 +186,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 +221,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<GridView> elementMap_;
IntersectionMapper intersectionMapper_;
+
std::vector<SubControlVolume> scvs_;
std::vector<SubControlVolumeFace> scvfs_;
std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
@@ -249,132 +247,7 @@ private:
template<class TypeTag>
class StaggeredFVGridGeometry<TypeTag, false>
{
- 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<IndexType> scvfsIndexSet;
- std::vector<IndexType> 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<IndexType>& scvfIndicesOfScv(IndexType scvIdx) const
- { return scvfIndicesOfScv_[scvIdx]; }
-
- const std::vector<IndexType>& 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<GridView> elementMap_;
- std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
- std::vector<std::vector<IndexType>> neighborVolVarIndices_;
+ // TODO: implement without caching
};
} // end namespace
diff --git a/dumux/discretization/staggered/globalfacevariables.hh b/dumux/discretization/staggered/globalfacevariables.hh
index 1e63de4291..74c013decc 100644
--- a/dumux/discretization/staggered/globalfacevariables.hh
+++ b/dumux/discretization/staggered/globalfacevariables.hh
@@ -34,21 +34,31 @@ 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 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)
+ StaggeredGlobalFaceVariables(const Problem& problem) : problemPtr_(&problem) {}
+
+ void update(const FVGridGeometry& fvGridGeometry, const SolutionVector& sol)
{
- problemPtr_ = &problem;
- faceVariables_.resize(problem.model().numFaceDofs());
- assert(faceVariables_.size() == sol.size());
+ const auto& faceSol = sol[faceIdx];
+
+ const auto numFaceDofs = fvGridGeometry.gridView().size(1);
+
+ faceVariables_.resize(numFaceDofs);
+ assert(faceVariables_.size() == faceSol.size());
- for(int i = 0; i < problem.model().numFaceDofs(); ++i)
+ for(int i = 0; i < numFaceDofs; ++i)
{
- faceVariables_[i].update(sol[i]);
+ faceVariables_[i].update(faceSol[i]);
}
}
diff --git a/dumux/discretization/staggered/globalfluxvariablescache.hh b/dumux/discretization/staggered/globalfluxvariablescache.hh
index 5e4dfd38cf..9e9a194ba7 100644
--- a/dumux/discretization/staggered/globalfluxvariablescache.hh
+++ b/dumux/discretization/staggered/globalfluxvariablescache.hh
@@ -47,32 +47,38 @@ class StaggeredGlobalFluxVariablesCache<TypeTag, true>
friend StaggeredElementFluxVariablesCache<TypeTag, true>;
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<FluxVariablesCache> fluxVarsCache_;
diff --git a/dumux/discretization/staggered/globalvolumevariables.hh b/dumux/discretization/staggered/globalvolumevariables.hh
index 04477eb7bd..003d844109 100644
--- a/dumux/discretization/staggered/globalvolumevariables.hh
+++ b/dumux/discretization/staggered/globalvolumevariables.hh
@@ -43,8 +43,7 @@ class StaggeredGlobalVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/true>
{
// The local class needs to access and change volVars
friend StaggeredElementVolumeVariables<TypeTag, true>;
- // 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<TypeTag>;
friend typename GET_PROP_TYPE(TypeTag, PrimaryVariableSwitch);
@@ -52,6 +51,7 @@ class StaggeredGlobalVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/true>
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 +71,17 @@ class StaggeredGlobalVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/true>
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 +89,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 +99,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 +108,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 +118,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 +143,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> volumeVariables_;
diff --git a/dumux/freeflow/staggered/fluxvariables.hh b/dumux/freeflow/staggered/fluxvariables.hh
index 0d941b9df6..4635bb46e1 100644
--- a/dumux/freeflow/staggered/fluxvariables.hh
+++ b/dumux/freeflow/staggered/fluxvariables.hh
@@ -128,7 +128,6 @@ public:
}
void computeCellCenterToCellCenterStencil(Stencil& stencil,
- const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
@@ -141,7 +140,6 @@ public:
}
void computeCellCenterToFaceStencil(Stencil& stencil,
- const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
@@ -150,7 +148,6 @@ public:
}
void computeFaceToCellCenterStencil(Stencil& stencil,
- const Problem& problem,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
{
@@ -167,7 +164,6 @@ public:
}
void computeFaceToFaceStencil(Stencil& stencil,
- const Problem& problem,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
{
diff --git a/dumux/freeflow/staggered/localresidual.hh b/dumux/freeflow/staggered/localresidual.hh
index 4d9ab7b95c..1a4a112db5 100644
--- a/dumux/freeflow/staggered/localresidual.hh
+++ b/dumux/freeflow/staggered/localresidual.hh
@@ -115,10 +115,8 @@ class StaggeredNavierStokesResidualImpl<TypeTag, false> : public Dumux::Staggere
static constexpr bool normalizePressure = GET_PROP_VALUE(TypeTag, NormalizePressure);
public:
- // copying the local residual class is not a good idea
- StaggeredNavierStokesResidualImpl(const StaggeredNavierStokesResidualImpl &) = delete;
- StaggeredNavierStokesResidualImpl() = default;
+ using ParentType::ParentType;
CellCenterPrimaryVariables computeFluxForCellCenter(const Element &element,
diff --git a/dumux/freeflow/staggered/problem.hh b/dumux/freeflow/staggered/problem.hh
index d9d5b54486..081bebb8ad 100644
--- a/dumux/freeflow/staggered/problem.hh
+++ b/dumux/freeflow/staggered/problem.hh
@@ -23,9 +23,9 @@
#ifndef DUMUX_NAVIERSTOKES_PROBLEM_HH
#define DUMUX_NAVIERSTOKES_PROBLEM_HH
-#include <dumux/implicit/problem.hh>
-
+#include <dumux/implicit/properties.hh>
#include "properties.hh"
+#include <dumux/common/staggeredfvproblem.hh>
namespace Dumux
{
@@ -37,21 +37,22 @@ namespace Dumux
* This implements gravity (if desired) and a function returning the temperature.
*/
template<class TypeTag>
-class NavierStokesProblem : public ImplicitProblem<TypeTag>
+class NavierStokesProblem : public StaggeredFVProblem<TypeTag>
{
- typedef ImplicitProblem<TypeTag> ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Implementation;
+ using ParentType = StaggeredFVProblem<TypeTag>;
+ using Implementation = typename GET_PROP_TYPE(TypeTag, Problem);
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
- typedef typename GridView::Grid Grid;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager);
+ using Grid = typename GridView::Grid;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::Intersection Intersection;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Intersection = typename GridView::Intersection;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
@@ -75,11 +76,11 @@ class NavierStokesProblem : public ImplicitProblem<TypeTag>
typename DofTypeIndices::FaceIdx faceIdx;
public:
- NavierStokesProblem(TimeManager &timeManager, const GridView &gridView)
- : ParentType(timeManager, gridView),
+ NavierStokesProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry),
gravity_(0)
{
- if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
+ if (getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity"))
gravity_[dim-1] = -9.81;
}
diff --git a/dumux/freeflow/staggered/propertydefaults.hh b/dumux/freeflow/staggered/propertydefaults.hh
index 4206e1fa9c..0aad12938c 100644
--- a/dumux/freeflow/staggered/propertydefaults.hh
+++ b/dumux/freeflow/staggered/propertydefaults.hh
@@ -38,6 +38,7 @@
#include "fluxvariablescache.hh"
#include "velocityoutput.hh"
#include "vtkoutputmodule.hh"
+#include "vtkoutputfields.hh"
#include "boundarytypes.hh"
#include <dumux/implicit/staggered/localresidual.hh>
@@ -158,6 +159,8 @@ SET_PROP(NavierStokes, FluidState){
// disable velocity output by default
SET_BOOL_PROP(NavierStokes, VtkAddVelocity, true);
+SET_TYPE_PROP(NavierStokes, VtkOutputFields, NavierStokesVtkOutputFields<TypeTag>);
+
// enable gravity by default
SET_BOOL_PROP(NavierStokes, ProblemEnableGravity, true);
diff --git a/dumux/freeflow/staggered/velocityoutput.hh b/dumux/freeflow/staggered/velocityoutput.hh
index 2fa4415f30..36951d7ef3 100644
--- a/dumux/freeflow/staggered/velocityoutput.hh
+++ b/dumux/freeflow/staggered/velocityoutput.hh
@@ -57,6 +57,10 @@ class StaggeredFreeFlowVelocityOutput
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+
static const bool isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox);
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
@@ -75,11 +79,17 @@ public:
*
* \param problem The problem to be solved
*/
- StaggeredFreeFlowVelocityOutput(const Problem& problem)
+ StaggeredFreeFlowVelocityOutput(const Problem& problem,
+ const FVGridGeometry& fvGridGeometry,
+ const GridVariables& gridVariables,
+ const SolutionVector& sol)
: problem_(problem)
+ , fvGridGeometry_(fvGridGeometry)
+ , gridVariables_(gridVariables)
+ , sol_(sol)
{
// check, if velocity output can be used (works only for cubes so far)
- velocityOutput_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddVelocity);
+ velocityOutput_ = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Vtk.AddVelocity");
}
bool enableOutput()
@@ -103,7 +113,7 @@ public:
for (auto&& scvf : scvfs(fvGeometry))
{
- auto& origFaceVars = problem_.model().curGlobalFaceVars().faceVars(scvf.dofIndex());
+ auto& origFaceVars = gridVariables_.curGridFaceVars().faceVars(scvf.dofIndex());
auto dirIdx = scvf.directionIndex();
velocity[dofIdxGlobal][dirIdx] += 0.5*origFaceVars.velocity();
}
@@ -113,6 +123,9 @@ public:
private:
const Problem& problem_;
+ const FVGridGeometry& fvGridGeometry_;
+ const GridVariables& gridVariables_;
+ const SolutionVector& sol_;
bool velocityOutput_;
};
diff --git a/dumux/freeflow/staggered/vtkoutputfields.hh b/dumux/freeflow/staggered/vtkoutputfields.hh
new file mode 100644
index 0000000000..6708e6a8bd
--- /dev/null
+++ b/dumux/freeflow/staggered/vtkoutputfields.hh
@@ -0,0 +1,61 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Adds vtk output fields specific to the twop model
+ */
+#ifndef DUMUX_NAVIER_STOKES_VTK_OUTPUT_FIELDS_HH
+#define DUMUX_NAVIER_STOKES_VTK_OUTPUT_FIELDS_HH
+
+#include <dumux/implicit/properties.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoP, InputOutput
+ * \brief Adds vtk output fields specific to the twop model
+ */
+template<class TypeTag>
+class NavierStokesVtkOutputFields
+{
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+public:
+ template <class VtkOutputModule>
+ static void init(VtkOutputModule& vtk)
+ {
+ // vtk.addSecondaryVariable("Sw", [](const VolumeVariables& v){ return v.saturation(Indices::wPhaseIdx); });
+ // vtk.addSecondaryVariable("Sn", [](const VolumeVariables& v){ return v.saturation(Indices::nPhaseIdx); });
+ // vtk.addSecondaryVariable("pw", [](const VolumeVariables& v){ return v.pressure(Indices::wPhaseIdx); });
+ // vtk.addSecondaryVariable("pn", [](const VolumeVariables& v){ return v.pressure(Indices::nPhaseIdx); });
+ vtk.addSecondaryVariable("p", [](const VolumeVariables& v){ return v.pressure(); });
+ // vtk.addSecondaryVariable("pc", [](const VolumeVariables& v){ return v.capillaryPressure(); });
+ // vtk.addSecondaryVariable("rhoW", [](const VolumeVariables& v){ return v.density(Indices::wPhaseIdx); });
+ // vtk.addSecondaryVariable("rhoN", [](const VolumeVariables& v){ return v.density(Indices::nPhaseIdx); });
+ // vtk.addSecondaryVariable("mobW", [](const VolumeVariables& v){ return v.mobility(Indices::wPhaseIdx); });
+ // vtk.addSecondaryVariable("mobN", [](const VolumeVariables& v){ return v.mobility(Indices::nPhaseIdx); });
+ // vtk.addSecondaryVariable("temperature", [](const VolumeVariables& v){ return v.temperature(); });
+ // vtk.addSecondaryVariable("porosity", [](const VolumeVariables& v){ return v.porosity(); });
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/implicit/staggered/assembler.hh b/dumux/implicit/staggered/assembler.hh
deleted file mode 100644
index 30a84437f4..0000000000
--- a/dumux/implicit/staggered/assembler.hh
+++ /dev/null
@@ -1,174 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \brief An assembler for the global Jacobian matrix for fully implicit models.
- */
-#ifndef DUMUX_STAGGERED_ASSEMBLER_HH
-#define DUMUX_STAGGERED_ASSEMBLER_HH
-
-#include <dumux/implicit/properties.hh>
-#include <dumux/implicit/assembler.hh>
-#include <dune/istl/matrixindexset.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup ImplicitModel
- * \brief An assembler for the global Jacobian matrix for fully implicit models.
- */
-template<class TypeTag>
-class StaggeredAssembler : public ImplicitAssembler<TypeTag>
-{
- typedef ImplicitAssembler<TypeTag> ParentType;
- friend class ImplicitAssembler<TypeTag>;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
-
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::IndexSet::IndexType IndexType;
-
- typedef typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockCCToCC CCToCCMatrixBlock;
- typedef typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockCCToFace CCToFaceMatrixBlock;
-
- typedef typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockFaceToFace FaceToFaceMatrixBlock;
- typedef typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockFaceToCC FaceToCCMatrixBlock;
-
- using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
- typename DofTypeIndices::CellCenterIdx cellCenterIdx;
- typename DofTypeIndices::FaceIdx faceIdx;
-
-public:
-
- /*!
- * \brief Initialize the jacobian assembler.
- *
- * At this point we can assume that all objects in the problem and
- * the model have been allocated. We can not assume that they are
- * fully initialized, though.
- *
- * \param problem The problem object
- */
- void init(Problem& problem)
- {
- std::cout << "init(Problem& problem)" << std::endl;
- this->problemPtr_ = &problem;
-
- // initialize the BCRS matrix
- createMatrix_();
-
- // initialize the jacobian matrix with zeros
- *this->matrix_ = 0;
-
- // allocate the residual vector
- this->residual_[cellCenterIdx].resize(problem.model().numCellCenterDofs());
- this->residual_[faceIdx].resize(problem.model().numFaceDofs());
-// printmatrix(std::cout, matrix(), "", "");
- }
-
-
-private:
-
- // Construct the multitype matrix for the global jacobian
- void createMatrix_()
- {
- // create the multitype matrix
- this->matrix_ = std::make_shared<JacobianMatrix>();
-
- // get sub matrix sizes
- const auto cellCenterSize = this->problem_().model().numCellCenterDofs();
- const auto faceSize = this->problem_().model().numFaceDofs();
-
- // allocate the sub matrices (BCRS matrices)
- auto A11 = CCToCCMatrixBlock(cellCenterSize, cellCenterSize, CCToCCMatrixBlock::random);
- auto A12 = CCToFaceMatrixBlock(cellCenterSize, faceSize, CCToFaceMatrixBlock::random);
-
- auto A21 = FaceToCCMatrixBlock(faceSize, cellCenterSize, FaceToCCMatrixBlock::random);
- auto A22 = FaceToFaceMatrixBlock(faceSize, faceSize, FaceToFaceMatrixBlock::random);
-
- setupMatrices_(A11, A12, A21, A22);
-
- (*this->matrix_)[cellCenterIdx][cellCenterIdx] = A11;
- (*this->matrix_)[cellCenterIdx][faceIdx] = A12;
- (*this->matrix_)[faceIdx][cellCenterIdx] = A21;
- (*this->matrix_)[faceIdx][faceIdx] = A22;
-
-// printmatrix(std::cout, A11, "A11", "");
-// printmatrix(std::cout, A12, "A12", "");
-// printmatrix(std::cout, A21, "A21", "");
-// printmatrix(std::cout, A22, "A22", "");
- }
-
- void setupMatrices_(CCToCCMatrixBlock &A11, CCToFaceMatrixBlock &A12,
- FaceToCCMatrixBlock &A21, FaceToFaceMatrixBlock &A22)
- {
- // get sub matrix sizes
- const auto numDofsCC = this->problem_().model().numCellCenterDofs();
- const auto numDofsFace = this->problem_().model().numFaceDofs();
-
- // get occupation pattern of the matrix
- 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& assemblyMap = this->problem_().model().localJacobian().assemblyMap();
-
- for (const auto& element : elements(this->gridView_()))
- {
- // the global index of the element at hand
- const auto ccGlobalI = this->elementMapper_().index(element);
-
- for (auto&& ccGlobalJ : assemblyMap(cellCenterIdx, cellCenterIdx, ccGlobalI))
- occupationPatternA11.add(ccGlobalI, ccGlobalJ);
- for (auto&& faceGlobalJ : assemblyMap(cellCenterIdx, faceIdx, ccGlobalI))
- occupationPatternA12.add(ccGlobalI, faceGlobalJ);
-
- auto fvGeometry = localView(this->problem_().model().fvGridGeometry());
- fvGeometry.bindElement(element);
-
- // loop over sub control faces
- for (auto&& scvf : scvfs(fvGeometry))
- {
- const auto faceGlobalI = scvf.dofIndex();
- for (auto&& ccGlobalJ : assemblyMap(faceIdx, cellCenterIdx, scvf.index()))
- occupationPatternA21.add(faceGlobalI, ccGlobalJ);
- for (auto&& faceGlobalJ : assemblyMap(faceIdx, faceIdx, scvf.index()))
- occupationPatternA22.add(faceGlobalI, faceGlobalJ);
- }
- }
- // export patterns to matrices
- occupationPatternA11.exportIdx(A11);
- occupationPatternA12.exportIdx(A12);
- occupationPatternA21.exportIdx(A21);
- occupationPatternA22.exportIdx(A22);
- }
-
-
-};
-
-} // namespace Dumux
-
-#endif
diff --git a/dumux/implicit/staggered/gridvariables.hh b/dumux/implicit/staggered/gridvariables.hh
new file mode 100644
index 0000000000..1140c791f1
--- /dev/null
+++ b/dumux/implicit/staggered/gridvariables.hh
@@ -0,0 +1,115 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Class storing scv and scvf variables
+ */
+#ifndef DUMUX_STAGGERED_GRID_VARIABLES_HH
+#define DUMUX_STAGGERED_GRID_VARIABLES_HH
+
+#include <dumux/implicit/gridvariables.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief Class storing scv and scvf variables
+ */
+template<class TypeTag>
+class StaggeredGridVariables : public GridVariables<TypeTag>
+{
+ using ParentType = GridVariables<TypeTag>;
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GlobalVolumeVariables);
+ using GridFaceVariables = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
+ using GridFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, GlobalFluxVariablesCache);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+
+public:
+ //! Constructor
+ StaggeredGridVariables(std::shared_ptr<Problem> problem, std::shared_ptr<FVGridGeometry> fvGridGeometry)
+ : ParentType(problem, fvGridGeometry)
+ , fvGridGeometry_(fvGridGeometry)
+ , curGridFaceVariables_(*problem)
+ , prevGridFaceVariables_(*problem)
+ {}
+
+ //! update all variables
+ void update(const SolutionVector& curSol)
+ {
+ ParentType::update(curSol);
+ curGridFaceVariables_.update(*fvGridGeometry_, curSol);
+ }
+
+ //! initialize all variables (stationary case)
+ void init(const SolutionVector& curSol)
+ {
+ ParentType::init(curSol);
+ curGridFaceVariables_.update(*fvGridGeometry_, curSol);
+ }
+
+ //! initialize all variables (instationary case)
+ void init(const SolutionVector& curSol, const SolutionVector& initSol)
+ {
+ ParentType::init(curSol, initSol);
+ curGridFaceVariables_.update(*fvGridGeometry_, curSol);
+ prevGridFaceVariables_.update(*fvGridGeometry_, initSol);
+ }
+
+ //! Sets the current state as the previous for next time step
+ //! this has to be called at the end of each time step
+ void advanceTimeStep()
+ {
+ ParentType::advanceTimeStep();
+ prevGridFaceVariables_ = curGridFaceVariables_;
+ }
+
+ //! resets state to the one before time integration
+ void resetTimeStep(const SolutionVector& solution)
+ {
+ ParentType::resetTimeStep(solution);
+ curGridFaceVariables_ = prevGridFaceVariables_;
+ }
+
+ const GridFaceVariables& curGridFaceVars() const
+ { return curGridFaceVariables_; }
+
+ const GridFaceVariables& prevGridFaceVars() const
+ { return prevGridFaceVariables_; }
+
+ GridFaceVariables& curGridFaceVars()
+ { return curGridFaceVariables_; }
+
+ GridFaceVariables& prevGridFaceVars()
+ { return prevGridFaceVariables_; }
+
+private:
+
+ std::shared_ptr<FVGridGeometry> fvGridGeometry_;
+
+ GridFaceVariables curGridFaceVariables_;
+ GridFaceVariables prevGridFaceVariables_;
+
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/implicit/staggered/localjacobian.hh b/dumux/implicit/staggered/localjacobian.hh
deleted file mode 100644
index 571956299e..0000000000
--- a/dumux/implicit/staggered/localjacobian.hh
+++ /dev/null
@@ -1,594 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \brief Caculates the Jacobian of the local residual for fully-implicit models
- */
-#ifndef DUMUX_STAGGERED_LOCAL_JACOBIAN_HH
-#define DUMUX_STAGGERED_LOCAL_JACOBIAN_HH
-
-#include <dune/istl/io.hh>
-#include <dune/istl/matrix.hh>
-
-#include <dumux/common/math.hh>
-#include <dumux/common/valgrind.hh>
-
-#include <dumux/implicit/properties.hh>
-#include "primaryvariables.hh"
-#include "assemblymap.hh"
-
-namespace Dumux
-{
-/*!
- * \ingroup ImplicitLocalJacobian
- * \brief Calculates the Jacobian of the local residual for fully-implicit models
- *
- * The default behavior is to use numeric differentiation, i.e.
- * forward or backward differences (2nd order), or central
- * differences (3rd order). The method used is determined by the
- * "NumericDifferenceMethod" property:
- *
- * - if the value of this property is smaller than 0, backward
- * differences are used, i.e.:
- * \f[
- \frac{\partial f(x)}{\partial x} \approx \frac{f(x) - f(x - \epsilon)}{\epsilon}
- * \f]
- *
- * - if the value of this property is 0, central
- * differences are used, i.e.:
- * \f[
- \frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x - \epsilon)}{2 \epsilon}
- * \f]
- *
- * - if the value of this property is larger than 0, forward
- * differences are used, i.e.:
- * \f[
- \frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x)}{\epsilon}
- * \f]
- *
- * Here, \f$ f \f$ is the residual function for all equations, \f$x\f$
- * is the value of a sub-control volume's primary variable at the
- * evaluation point and \f$\epsilon\f$ is a small value larger than 0.
- */
-template<class TypeTag>
-class StaggeredLocalJacobian
-{
- using Implementation = typename GET_PROP_TYPE(TypeTag, LocalJacobian);
- using JacobianAssembler = typename GET_PROP_TYPE(TypeTag, JacobianAssembler);
- using LocalResidual = typename GET_PROP_TYPE(TypeTag, LocalResidual);
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using Model = typename GET_PROP_TYPE(TypeTag, Model);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
-
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
- using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
- using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
- using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
- using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using Element = typename GridView::template Codim<0>::Entity;
- using IndexType = typename GridView::IndexSet::IndexType;
-
- enum {
- numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter),
- numEqFace = GET_PROP_VALUE(TypeTag, NumEqFace),
- };
-
- using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
- typename DofTypeIndices::CellCenterIdx cellCenterIdx;
- typename DofTypeIndices::FaceIdx faceIdx;
-
- using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
- using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
- using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
-
- using FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
-
- using PriVarIndices = typename Dumux::PriVarIndices<TypeTag>;
-
- using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
-
- using AssemblyMap = Dumux::StaggeredAssemblyMap<TypeTag>;
-
-
-public:
-
- StaggeredLocalJacobian(const StaggeredLocalJacobian &) = delete;
-
- StaggeredLocalJacobian()
- {
- numericDifferenceMethod_ = GET_PARAM_FROM_GROUP(TypeTag, int, Implicit, NumericDifferenceMethod);
- }
-
- /*!
- * \brief Initialize the local Jacobian object.
- *
- * At this point we can assume that everything has been allocated,
- * although some objects may not yet be completely initialized.
- *
- * \param problem The problem which we want to simulate.
- */
- void init(Problem &problem)
- {
- problemPtr_ = &problem;
- localResidual_.init(problem);
- assemblyMap_.init(problem);
- }
-
- /*!
- * \brief Assemble an element's local Jacobian matrix of the
- * defect.
- *
- * \param element The DUNE Codim<0> entity which we look at.
- */
- void assemble(const Element& element, JacobianMatrix& matrix, SolutionVector& residual)
- {
- const bool isGhost = (element.partitionType() == Dune::GhostEntity);
- if(isGhost)
- DUNE_THROW(Dune::NotImplemented, "Support for ghost cells not implemented");
-
- // prepare the volvars/fvGeometries in case caching is disabled
- auto fvGeometry = localView(this->model_().fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(this->model_().curGlobalVolVars());
- curElemVolVars.bind(element, fvGeometry, this->model_().curSol());
-
- auto prevElemVolVars = localView(this->model_().prevGlobalVolVars());
- prevElemVolVars.bindElement(element, fvGeometry, this->model_().prevSol());
-
- auto elemFluxVarsCache = localView(this->model_().globalFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
-
- // set the actual cell center dof index
- ccGlobalI_ = this->problem_().elementMapper().index(element);
-
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(this->problem_(), element, fvGeometry);
-
- auto& curGlobalFaceVars = getCurrNonConstGlobalFaceVars_();
- auto& prevGlobalFaceVars = getPrevGlobalFaceVars_();
-
- // calculate the local residual for all dofs of this element
- this->localResidual().eval(element, fvGeometry,
- prevElemVolVars, curElemVolVars,
- prevGlobalFaceVars, curGlobalFaceVars,
- elemBcTypes, elemFluxVarsCache);
- // store the cell center residual in global container
- residual[cellCenterIdx][ccGlobalI_] = this->localResidual().ccResidual();
-
- // 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))
- {
- residual[faceIdx][scvf.dofIndex()] += this->localResidual().faceResidual(scvf.localFaceIdx());
- faceResidualCache[scvf.localFaceIdx()] = this->localResidual().faceResidual(scvf.localFaceIdx());
- }
-
- this->model_().updatePVWeights(fvGeometry);
-
- // calculate derivatives of all dofs in stencil with respect to the dofs in the element
- evalPartialDerivatives_(element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- prevGlobalFaceVars,
- curGlobalFaceVars,
- elemFluxVarsCache,
- elemBcTypes,
- matrix,
- residual[cellCenterIdx][ccGlobalI_],
- faceResidualCache);
- }
-
- /*!
- * \brief Returns a reference to the object which calculates the
- * local residual.
- */
- const LocalResidual &localResidual() const
- { return localResidual_; }
-
- /*!
- * \brief Returns a reference to the object which calculates the
- * local residual.
- */
- LocalResidual &localResidual()
- { return localResidual_; }
-
- const AssemblyMap& assemblyMap() const
- { return assemblyMap_; }
-
-protected:
- void evalPartialDerivatives_(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& prevElemVolVars,
- ElementVolumeVariables& curElemVolVars,
- const GlobalFaceVars& prevGlobalFaceVars,
- GlobalFaceVars& curGlobalFaceVars,
- ElementFluxVariablesCache& elemFluxVarsCache,
- const ElementBoundaryTypes& elemBcTypes,
- JacobianMatrix& matrix,
- const CellCenterPrimaryVariables& ccResidual,
- const FaceSolutionVector& faceResidualCache)
- {
- // compute the derivatives of the cell center residuals with respect to cell center dofs
- dCCdCC_(element, fvGeometry, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, ccResidual);
-
- // compute the derivatives of the cell center residuals with respect to face dofs
- dCCdFace_(element, fvGeometry, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, ccResidual);
-
- // compute the derivatives of the face residuals with respect to cell center dofs
- dFacedCC_(element, fvGeometry, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, faceResidualCache);
-
- // compute the derivatives of the face residuals with respect to face dofs
- dFacedFace_(element, fvGeometry, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, faceResidualCache);
- }
-
- /*!
- * \brief Computes the derivatives of the cell center residuals with respect to cell center dofs
- */
- void dCCdCC_(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& prevElemVolVars,
- ElementVolumeVariables& curElemVolVars,
- const GlobalFaceVars& prevGlobalFaceVars,
- const GlobalFaceVars& curGlobalFaceVars,
- ElementFluxVariablesCache& elemFluxVarsCache,
- const ElementBoundaryTypes& elemBcTypes,
- JacobianMatrix& matrix,
- const CellCenterPrimaryVariables& ccResidual)
- {
- // build derivatives with for cell center dofs w.r.t. cell center dofs
- auto&& scvI = fvGeometry.scv(ccGlobalI_);
-
- for(const auto& globalJ : assemblyMap_(cellCenterIdx, cellCenterIdx, ccGlobalI_))
- {
- // 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 = getCurVolVars(curElemVolVars, scvJ);
- VolumeVariables origVolVars(curVolVars);
-
- for(auto pvIdx : PriVarIndices(cellCenterIdx))
- {
- PrimaryVariables priVars(CellCenterPrimaryVariables(this->model_().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, this->problem_(), elementJ, scvJ);
-
- this->localResidual().evalCellCenter(element, fvGeometry, scvI,
- prevElemVolVars, curElemVolVars,
- prevGlobalFaceVars, curGlobalFaceVars,
- elemBcTypes, elemFluxVarsCache);
-
- auto partialDeriv = (this->localResidual().ccResidual() - ccResidual);
- partialDeriv /= eps;
-
- // update the global jacobian matrix with the current partial derivatives
- this->updateGlobalJacobian_(matrix[cellCenterIdx][cellCenterIdx], ccGlobalI_, globalJ, pvIdx, partialDeriv);
-
- // restore the original volVars
- curVolVars = origVolVars;
- }
- }
- }
-
- /*!
- * \brief Computes the derivatives of the cell center residuals with respect to face dofs
- */
- void dCCdFace_(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& prevElemVolVars,
- const ElementVolumeVariables& curElemVolVars,
- const GlobalFaceVars& prevGlobalFaceVars,
- GlobalFaceVars& curGlobalFaceVars,
- ElementFluxVariablesCache& elemFluxVarsCache,
- const ElementBoundaryTypes& elemBcTypes,
- JacobianMatrix& matrix,
- const CellCenterPrimaryVariables& ccResidual)
- {
- // build derivatives with for cell center dofs w.r.t. face dofs
- auto&& scvI = fvGeometry.scv(ccGlobalI_);
-
- for(const auto& globalJ : assemblyMap_(cellCenterIdx, faceIdx, ccGlobalI_))
- {
- // get the faceVars of the face with respect to which we are going to build the derivative
- auto origFaceVars = curGlobalFaceVars.faceVars(globalJ);
- auto& curFaceVars = curGlobalFaceVars.faceVars(globalJ);
-
- for(auto pvIdx : PriVarIndices(faceIdx))
- {
- PrimaryVariables priVars(CellCenterPrimaryVariables(0.0), FacePrimaryVariables(this->model_().curSol()[faceIdx][globalJ]));
-
- const Scalar eps = numericEpsilon(priVars[pvIdx], cellCenterIdx, faceIdx);
- priVars[pvIdx] += eps;
- curFaceVars.update(priVars[faceIdx]);
-
- this->localResidual().evalCellCenter(element, fvGeometry, scvI,
- prevElemVolVars, curElemVolVars,
- prevGlobalFaceVars, curGlobalFaceVars,
- elemBcTypes, elemFluxVarsCache);
-
- auto partialDeriv = (this->localResidual().ccResidual() - ccResidual);
- partialDeriv /= eps;
-
- // update the global jacobian matrix with the current partial derivatives
- this->updateGlobalJacobian_(matrix[cellCenterIdx][faceIdx], ccGlobalI_, globalJ, pvIdx - Indices::faceOffset, partialDeriv);
-
- // restore the original faceVars
- curFaceVars = origFaceVars;
- }
- }
- }
-
- /*!
- * \brief Computes the derivatives of the face residuals with respect to cell center dofs
- */
- void dFacedCC_(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& prevElemVolVars,
- ElementVolumeVariables& curElemVolVars,
- const GlobalFaceVars& prevGlobalFaceVars,
- const GlobalFaceVars& curGlobalFaceVars,
- ElementFluxVariablesCache& elemFluxVarsCache,
- const ElementBoundaryTypes& elemBcTypes,
- JacobianMatrix& matrix,
- const FaceSolutionVector& cachedResidual)
- {
- 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 : assemblyMap_(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 = getCurVolVars(curElemVolVars, scvJ);
- VolumeVariables origVolVars(curVolVars);
-
- for(auto pvIdx : PriVarIndices(cellCenterIdx))
- {
- PrimaryVariables priVars(CellCenterPrimaryVariables(this->model_().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, this->problem_(), elementJ, scvJ);
-
- this->localResidual().evalFace(element, fvGeometry, scvf,
- prevElemVolVars, curElemVolVars,
- prevGlobalFaceVars, curGlobalFaceVars,
- elemBcTypes, elemFluxVarsCache);
-
- auto partialDeriv = (this->localResidual().faceResidual(scvf.localFaceIdx()) - cachedResidual[scvf.localFaceIdx()]);
- partialDeriv /= eps;
- // update the global jacobian matrix with the current partial derivatives
- this->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
- */
- void dFacedFace_(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& prevElemVolVars,
- const ElementVolumeVariables& curElemVolVars,
- const GlobalFaceVars& prevGlobalFaceVars,
- GlobalFaceVars& curGlobalFaceVars,
- ElementFluxVariablesCache& elemFluxVarsCache,
- const ElementBoundaryTypes& elemBcTypes,
- JacobianMatrix& matrix,
- const FaceSolutionVector& cachedResidual)
- {
- 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 : assemblyMap_(faceIdx, faceIdx, scvf.index()))
- {
- // get the faceVars of the face with respect to which we are going to build the derivative
- auto origFaceVars = curGlobalFaceVars.faceVars(globalJ);
- auto& curFaceVars = curGlobalFaceVars.faceVars(globalJ);
-
- for(auto pvIdx : PriVarIndices(faceIdx))
- {
- PrimaryVariables priVars(CellCenterPrimaryVariables(0.0), FacePrimaryVariables(this->model_().curSol()[faceIdx][globalJ]));
-
- const Scalar eps = numericEpsilon(priVars[pvIdx], faceIdx, faceIdx);
- priVars[pvIdx] += eps;
- curFaceVars.update(priVars[faceIdx]);
-
- this->localResidual().evalFace(element, fvGeometry, scvf,
- prevElemVolVars, curElemVolVars,
- prevGlobalFaceVars, curGlobalFaceVars,
- elemBcTypes, elemFluxVarsCache);
-
- auto partialDeriv = (this->localResidual().faceResidual(scvf.localFaceIdx()) - cachedResidual[scvf.localFaceIdx()]);
- partialDeriv /= eps;
-
- // update the global jacobian matrix with the current partial derivatives
- this->updateGlobalJacobian_(matrix[faceIdx][faceIdx], faceGlobalI, globalJ, pvIdx - Indices::faceOffset, partialDeriv);
-
- // restore the original faceVars
- curFaceVars = origFaceVars;
- }
- }
- }
- }
-
- /*!
- * \brief Returns a reference to the problem.
- */
- const Problem &problem_() const
- {
- Valgrind::CheckDefined(problemPtr_);
- return *problemPtr_;
- }
-
- /*!
- * \brief Returns a reference to the problem.
- */
- Problem &problem_()
- {
- Valgrind::CheckDefined(problemPtr_);
- return *problemPtr_;
- }
-
- /*!
- * \brief Returns a reference to the grid view.
- */
- const GridView &gridView_() const
- { return problem_().gridView(); }
-
- /*!
- * \brief Returns a reference to the model.
- */
- const Model &model_() const
- { return problem_().model(); }
-
- /*!
- * \brief Returns a reference to the model.
- */
- Model &model_()
- { return problem_().model(); }
-
- /*!
- * \brief Returns a reference to the jacobian assembler.
- */
- const JacobianAssembler &jacAsm_() const
- { return model_().jacobianAssembler(); }
-
- /*!
- * \brief Return the epsilon used to calculate the numeric derivative of a localResidual w.r.t a certain priVar
- *
- * \param priVar The the value of primary varible w.r.t which to derivative of the localResidual is calculated
- * \param idx1 Indicates whether the the derivative is build for a cellCenter or face localResidual
- * \param idx2 Indicates whether the the derivative is build w.r.t a priVar living on a cellCenter or face
- */
- Scalar numericEpsilon(const Scalar priVar, const int idx1, const int idx2) const
- {
- // 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<Scalar>(std::numeric_limits<Scalar>::epsilon(), 1.0);
- */
-
- static const Scalar baseEps = baseEps_[idx1][idx2];
- assert(std::numeric_limits<Scalar>::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<class SubMatrix, class CCOrFacePrimaryVariables>
- 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]);
- }
- }
-
- //! If the global vol vars caching is enabled we have to modify the global volvar object
- template<class T = TypeTag>
- typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), VolumeVariables>::type&
- getCurVolVars(ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
- { return this->model_().nonConstCurGlobalVolVars().volVars(scv); }
-
- //! When global volume variables caching is disabled, return the local volvar object
- template<class T = TypeTag>
- typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), VolumeVariables>::type&
- getCurVolVars(ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
- { return elemVolVars[scv]; }
-
-
- //! Convenience function to get the current non-const global face vars. This is necessary for classes that inherit from this class.
- auto& getCurrNonConstGlobalFaceVars_()
- {
- return this->model_().nonConstCurFaceVars();
- }
-
- //! Convenience function to get the previous global face vars
- auto& getPrevGlobalFaceVars_()
- {
- return this->model_().prevGlobalFaceVars();
- }
-
- IndexType ccGlobalI_;
- int numericDifferenceMethod_;
-
- // The problem we would like to solve
- Problem *problemPtr_;
-
- LocalResidual localResidual_;
-
- AssemblyMap assemblyMap_;
-
- using BaseEpsilon = typename GET_PROP(TypeTag, BaseEpsilon);
- const std::array<std::array<Scalar, 2>, 2> baseEps_ = BaseEpsilon::getEps();
-};
-
-}
-
-#endif
diff --git a/dumux/implicit/staggered/localresidual.hh b/dumux/implicit/staggered/localresidual.hh
index fe799358cf..4b8d1be3e9 100644
--- a/dumux/implicit/staggered/localresidual.hh
+++ b/dumux/implicit/staggered/localresidual.hh
@@ -26,6 +26,8 @@
#include <dune/istl/matrix.hh>
#include <dumux/common/valgrind.hh>
+#include <dumux/common/capabilities.hh>
+#include <dumux/common/timeloop.hh>
#include "properties.hh"
@@ -42,7 +44,6 @@ namespace Dumux
template<class TypeTag>
class StaggeredLocalResidual
{
- using ParentType = ImplicitLocalResidual<TypeTag>;
friend class ImplicitLocalResidual<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
@@ -65,6 +66,7 @@ class StaggeredLocalResidual
using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
@@ -77,25 +79,17 @@ class StaggeredLocalResidual
dimWorld = GridView::dimensionworld
};
-public:
- // copying the local residual class is not a good idea
- StaggeredLocalResidual(const StaggeredLocalResidual &) = delete;
-
- StaggeredLocalResidual() = default;
+ using TimeLoop = TimeLoopBase<Scalar>;
- /*!
- * \brief Initialize the local residual.
- *
- * This assumes that all objects of the simulation have been fully
- * allocated but not necessarily initialized completely.
- *
- * \param problem The representation of the physical problem to be
- * solved.
- */
- void init(Problem &problem)
- { problemPtr_ = &problem; }
+public:
+ //! the constructor for stationary problems
+ StaggeredLocalResidual() : prevSol_(nullptr) {}
+ StaggeredLocalResidual(std::shared_ptr<TimeLoop> timeLoop)
+ : timeLoop_(timeLoop)
+ , prevSol_(nullptr)
+ {}
/*!
* \name User interface
@@ -111,33 +105,33 @@ public:
* \param element The DUNE Codim<0> entity for which the residual
* ought to be calculated
*/
- void eval(const Element &element)
- {
- // make sure FVElementGeometry and volume variables are bound to the element
- auto fvGeometry = localView(this->problem().model().fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(problem().model().curGlobalVolVars());
- curElemVolVars.bind(element, fvGeometry, problem().model().curSol());
-
- auto prevElemVolVars = localView(problem().model().prevGlobalVolVars());
- prevElemVolVars.bindElement(element, fvGeometry, problem().model().prevSol());
-
- auto elemFluxVarsCache = localView(problem().model().globalFluxVarsCache());
- elemFluxVarsCache.bindElement(element, fvGeometry, curElemVolVars);
-
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem(), element, fvGeometry);
-
- auto& curGlobalFaceVars = problem().model().curGlobalFaceVars();
- auto& prevGlobalFaceVars = problem().model().prevGlobalFaceVars();
-
-
- asImp_().eval(element, fvGeometry,
- prevElemVolVars, curElemVolVars,
- prevGlobalFaceVars, curGlobalFaceVars,
- bcTypes, elemFluxVarsCache);
- }
+ // void eval(const Element &element)
+ // {
+ // // make sure FVElementGeometry and volume variables are bound to the element
+ // auto fvGeometry = localView(this->problem().model().fvGridGeometry());
+ // fvGeometry.bind(element);
+ //
+ // auto curElemVolVars = localView(problem().model().curGlobalVolVars());
+ // curElemVolVars.bind(element, fvGeometry, problem().model().curSol());
+ //
+ // auto prevElemVolVars = localView(problem().model().prevGlobalVolVars());
+ // prevElemVolVars.bindElement(element, fvGeometry, problem().model().prevSol());
+ //
+ // auto elemFluxVarsCache = localView(problem().model().globalFluxVarsCache());
+ // elemFluxVarsCache.bindElement(element, fvGeometry, curElemVolVars);
+ //
+ // ElementBoundaryTypes bcTypes;
+ // bcTypes.update(problem(), element, fvGeometry);
+ //
+ // auto& curGlobalFaceVars = problem().model().curGlobalFaceVars();
+ // auto& prevGlobalFaceVars = problem().model().prevGlobalFaceVars();
+ //
+ //
+ // asImp_().eval(element, fvGeometry,
+ // prevElemVolVars, curElemVolVars,
+ // prevGlobalFaceVars, curGlobalFaceVars,
+ // bcTypes, elemFluxVarsCache);
+ // }
/*!
* \brief Compute the local residual, i.e. the deviation of the
@@ -155,7 +149,8 @@ public:
* \param bcTypes The types of the boundary conditions for all
* vertices of the element
*/
- void eval(const Element &element,
+ void eval(const Problem& problem,
+ const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& prevElemVolVars,
const ElementVolumeVariables& curElemVolVars,
@@ -176,12 +171,12 @@ public:
for (auto&& scv : scvs(fvGeometry))
{
// treat the cell center dof
- evalCellCenter(element, fvGeometry, scv, prevElemVolVars, curElemVolVars, prevFaceVars, curFaceVars, bcTypes, elemFluxVarsCache);
+ evalCellCenter(problem, element, fvGeometry, scv, prevElemVolVars, curElemVolVars, prevFaceVars, curFaceVars, bcTypes, elemFluxVarsCache);
// now, treat the dofs on the facets:
for(auto&& scvf : scvfs(fvGeometry))
{
- evalFace(element, fvGeometry, scvf, prevElemVolVars, curElemVolVars, prevFaceVars, curFaceVars, bcTypes, elemFluxVarsCache);
+ evalFace(problem, element, fvGeometry, scvf, prevElemVolVars, curElemVolVars, prevFaceVars, curFaceVars, bcTypes, elemFluxVarsCache);
}
}
}
@@ -202,7 +197,8 @@ public:
* \param bcTypes The types of the boundary conditions for all
* vertices of the element
*/
- void evalCellCenter(const Element &element,
+ void evalCellCenter(const Problem& problem,
+ const Element &element,
const FVElementGeometry& fvGeometry,
const SubControlVolume& scv,
const ElementVolumeVariables& prevElemVolVars,
@@ -216,9 +212,9 @@ public:
ccResidual_ = 0.0;
ccStorageTerm_ = 0.0;
- asImp_().evalVolumeTermForCellCenter_(element, fvGeometry, scv, prevElemVolVars, curElemVolVars, prevFaceVars, curFaceVars, bcTypes);
- asImp_().evalFluxesForCellCenter_(element, fvGeometry, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
- asImp_().evalBoundaryForCellCenter_(element, fvGeometry, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
+ asImp_().evalVolumeTermForCellCenter_(problem, element, fvGeometry, scv, prevElemVolVars, curElemVolVars, prevFaceVars, curFaceVars, bcTypes);
+ asImp_().evalFluxesForCellCenter_(problem, element, fvGeometry, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
+ asImp_().evalBoundaryForCellCenter_(problem, element, fvGeometry, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
}
/*!
@@ -237,7 +233,8 @@ public:
* \param bcTypes The types of the boundary conditions for all
* vertices of the element
*/
- void evalFace(const Element &element,
+ void evalFace(const Problem& problem,
+ const Element &element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const ElementVolumeVariables& prevElemVolVars,
@@ -258,24 +255,12 @@ public:
faceResiduals_[scvf.localFaceIdx()] = 0.0;
faceStorageTerms_[scvf.localFaceIdx()] = 0.0;
- asImp_().evalVolumeTermForFace_(element, fvGeometry, scvf, prevElemVolVars, curElemVolVars, prevFaceVars, curFaceVars, bcTypes);
- asImp_().evalFluxesForFace_(element, fvGeometry, scvf, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
- asImp_().evalBoundaryForFace_(element, fvGeometry, scvf, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
+ asImp_().evalVolumeTermForFace_(problem, element, fvGeometry, scvf, prevElemVolVars, curElemVolVars, prevFaceVars, curFaceVars, bcTypes);
+ asImp_().evalFluxesForFace_(problem, element, fvGeometry, scvf, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
+ asImp_().evalBoundaryForFace_(problem, element, fvGeometry, scvf, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
}
- /*!
- * \brief Return the problem we are solving. Only call this after init()!
- */
- const Problem& problem() const
- { return *problemPtr_; }
-
- /*!
- * \brief Return the problem we are solving. Only call this after init()!
- */
- Problem& problem()
- { return *problemPtr_; }
-
const auto& ccResidual() const
{ return ccResidual_; }
@@ -291,7 +276,8 @@ protected:
/*!
* \brief Evaluate the flux terms for cell center dofs
*/
- void evalFluxesForCellCenter_(const Element& element,
+ void evalFluxesForCellCenter_(const Problem& problem,
+ const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const GlobalFaceVars& faceVars,
@@ -301,14 +287,15 @@ protected:
for (auto&& scvf : scvfs(fvGeometry))
{
if(!scvf.boundary())
- ccResidual_ += asImp_().computeFluxForCellCenter(element, fvGeometry, elemVolVars, faceVars, scvf, elemFluxVarsCache);
+ ccResidual_ += asImp_().computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars, faceVars, scvf, elemFluxVarsCache);
}
}
/*!
* \brief Evaluate the flux terms for face dofs
*/
- void evalFluxesForFace_(const Element& element,
+ void evalFluxesForFace_(const Problem& problem,
+ const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const ElementVolumeVariables& elemVolVars,
@@ -317,13 +304,14 @@ protected:
const ElementFluxVariablesCache& elemFluxVarsCache)
{
if(!scvf.boundary())
- faceResiduals_[scvf.localFaceIdx()] += asImp_().computeFluxForFace(element, scvf, fvGeometry, elemVolVars, globalFaceVars, elemFluxVarsCache);
+ faceResiduals_[scvf.localFaceIdx()] += asImp_().computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, globalFaceVars, elemFluxVarsCache);
}
/*!
* \brief Evaluate boundary conditions
*/
- void evalBoundary_(const Element& element,
+ void evalBoundary_(const Problem& problem,
+ const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const GlobalFaceVars& faceVars,
@@ -340,19 +328,20 @@ protected:
*/
template<class P = Problem>
typename std::enable_if<Dumux::Capabilities::isStationary<P>::value, void>::type
- evalVolumeTermForCellCenter_(const Element &element,
- const FVElementGeometry& fvGeometry,
- const SubControlVolume& scv,
- const ElementVolumeVariables& prevElemVolVars,
- const ElementVolumeVariables& curElemVolVars,
- const GlobalFaceVars& prevFaceVars,
- const GlobalFaceVars& curFaceVars,
- const ElementBoundaryTypes &bcTypes)
+ evalVolumeTermForCellCenter_(const Problem& problem,
+ const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolume& scv,
+ const ElementVolumeVariables& prevElemVolVars,
+ const ElementVolumeVariables& curElemVolVars,
+ const GlobalFaceVars& prevFaceVars,
+ const GlobalFaceVars& curFaceVars,
+ const ElementBoundaryTypes &bcTypes)
{
const auto curExtrusionFactor = curElemVolVars[scv].extrusionFactor();
// subtract the source term from the local rate
- CellCenterPrimaryVariables source = asImp_().computeSourceForCellCenter(element, fvGeometry, curElemVolVars, curFaceVars, scv);
+ CellCenterPrimaryVariables source = asImp_().computeSourceForCellCenter(problem, element, fvGeometry, curElemVolVars, curFaceVars, scv);
source *= scv.volume()*curExtrusionFactor;
ccResidual_ -= source;
@@ -363,17 +352,18 @@ protected:
*/
template<class P = Problem>
typename std::enable_if<Dumux::Capabilities::isStationary<P>::value, void>::type
- evalVolumeTermForFace_(const Element &element,
- const FVElementGeometry& fvGeometry,
- const SubControlVolumeFace& scvf,
- const ElementVolumeVariables& prevElemVolVars,
- const ElementVolumeVariables& curElemVolVars,
- const GlobalFaceVars& prevFaceVars,
- const GlobalFaceVars& curFaceVars,
- const ElementBoundaryTypes &bcTypes)
+ evalVolumeTermForFace_(const Problem& problem,
+ const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf,
+ const ElementVolumeVariables& prevElemVolVars,
+ const ElementVolumeVariables& curElemVolVars,
+ const GlobalFaceVars& prevFaceVars,
+ const GlobalFaceVars& curFaceVars,
+ const ElementBoundaryTypes &bcTypes)
{
// the source term:
- auto faceSource = asImp_().computeSourceForFace(scvf, curElemVolVars, curFaceVars);
+ auto faceSource = asImp_().computeSourceForFace(problem, scvf, curElemVolVars, curFaceVars);
const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
const auto curExtrusionFactor = curElemVolVars[scv].extrusionFactor();
faceSource *= 0.5*scv.volume()*curExtrusionFactor;
@@ -385,7 +375,8 @@ protected:
*/
template<class P = Problem>
typename std::enable_if<!Dumux::Capabilities::isStationary<P>::value, void>::type
- evalVolumeTermForCellCenter_(const Element &element,
+ evalVolumeTermForCellCenter_(const Problem& problem,
+ const Element &element,
const FVElementGeometry& fvGeometry,
const SubControlVolume& scv,
const ElementVolumeVariables& prevElemVolVars,
@@ -403,8 +394,8 @@ protected:
//
// We might need a more explicit way for
// doing the time discretization...
- auto prevCCStorage = asImp_().computeStorageForCellCenter(scv, prevVolVars);
- auto curCCStorage = asImp_().computeStorageForCellCenter(scv, curVolVars);
+ auto prevCCStorage = asImp_().computeStorageForCellCenter(problem, scv, prevVolVars);
+ auto curCCStorage = asImp_().computeStorageForCellCenter(problem, scv, curVolVars);
prevCCStorage *= prevVolVars.extrusionFactor();
curCCStorage *= curVolVars.extrusionFactor();
@@ -412,13 +403,13 @@ protected:
ccStorageTerm_ = std::move(curCCStorage);
ccStorageTerm_ -= std::move(prevCCStorage);
ccStorageTerm_ *= scv.volume();
- ccStorageTerm_ /= problem().timeManager().timeStepSize();
+ ccStorageTerm_ /= timeLoop_->timeStepSize();
// add the storage term to the residual
ccResidual_ += ccStorageTerm_;
// subtract the source term from the local rate
- CellCenterPrimaryVariables source = asImp_().computeSourceForCellCenter(element, fvGeometry, curElemVolVars, curFaceVars, scv);
+ CellCenterPrimaryVariables source = asImp_().computeSourceForCellCenter(problem, element, fvGeometry, curElemVolVars, curFaceVars, scv);
source *= scv.volume()*curVolVars.extrusionFactor();
ccResidual_ -= source;
@@ -429,7 +420,8 @@ protected:
*/
template<class P = Problem>
typename std::enable_if<!Dumux::Capabilities::isStationary<P>::value, void>::type
- evalVolumeTermForFace_(const Element &element,
+ evalVolumeTermForFace_(const Problem& problem,
+ const Element &element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const ElementVolumeVariables& prevElemVolVars,
@@ -441,18 +433,18 @@ protected:
const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
const auto& curVolVars = curElemVolVars[scv];
const auto& prevVolVars = prevElemVolVars[scv];
- auto prevFaceStorage = asImp_().computeStorageForFace(scvf, prevVolVars, prevFaceVars);
- auto curFaceStorage = asImp_().computeStorageForFace(scvf, curVolVars, curFaceVars);
+ auto prevFaceStorage = asImp_().computeStorageForFace(problem, scvf, prevVolVars, prevFaceVars);
+ auto curFaceStorage = asImp_().computeStorageForFace(problem, scvf, curVolVars, curFaceVars);
// the storage term
faceStorageTerms_[scvf.localFaceIdx()] = std::move(curFaceStorage);
faceStorageTerms_[scvf.localFaceIdx()] -= std::move(prevFaceStorage);
faceStorageTerms_[scvf.localFaceIdx()] *= (scv.volume()/2.0);
- faceStorageTerms_[scvf.localFaceIdx()] /= problem().timeManager().timeStepSize();
+ faceStorageTerms_[scvf.localFaceIdx()] /= timeLoop_->timeStepSize();
faceResiduals_[scvf.localFaceIdx()] += faceStorageTerms_[scvf.localFaceIdx()];
// the source term:
- auto faceSource = asImp_().computeSourceForFace(scvf, curElemVolVars, curFaceVars);
+ auto faceSource = asImp_().computeSourceForFace(problem, scvf, curElemVolVars, curFaceVars);
faceSource *= 0.5*scv.volume()*curVolVars.extrusionFactor();
faceResiduals_[scvf.localFaceIdx()] -= faceSource;
}
@@ -468,8 +460,46 @@ protected:
FaceSolutionVector faceResiduals_;
FaceSolutionVector faceStorageTerms_;
+
+
+ /*!
+ * \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)
+ { prevSol_ = &u; }
+
+ /*!
+ * \brief Return the solution that has been set as the previous one.
+ */
+ const SolutionVector& prevSol() const
+ {
+ assert(prevSol_ && "no solution set for storage term evaluation");
+ return *prevSol_;
+ }
+
+ /*!
+ * \brief If no solution has been set, we treat the problem as stationary.
+ */
+ bool isStationary() const
+ { return !prevSol_; }
+
+protected:
+ TimeLoop& timeLoop()
+ { return *timeLoop_; }
+
+ const TimeLoop& timeLoop() const
+ { return *timeLoop_; }
+
+ Implementation &asImp()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp() const
+ { return *static_cast<const Implementation*>(this); }
+
private:
- Problem* problemPtr_;
+ std::shared_ptr<TimeLoop> timeLoop_;
+ const SolutionVector* prevSol_;
};
diff --git a/dumux/implicit/staggered/newtoncontroller.hh b/dumux/implicit/staggered/newtoncontroller.hh
index 1e4cf1dd4a..d142371dce 100644
--- a/dumux/implicit/staggered/newtoncontroller.hh
+++ b/dumux/implicit/staggered/newtoncontroller.hh
@@ -32,7 +32,7 @@
#include <dumux/nonlinear/newtoncontroller.hh>
#include <dumux/linear/linearsolveracceptsmultitypematrix.hh>
#include <dumux/linear/matrixconverter.hh>
-#include "newtonconvergencewriter.hh"
+// #include "newtonconvergencewriter.hh"
namespace Dumux {
@@ -59,7 +59,7 @@ template <class TypeTag>
class StaggeredNewtonController : public NewtonController<TypeTag>
{
typedef NewtonController<TypeTag> ParentType;
- typedef NewtonConvergenceWriter<TypeTag> StaggeredNewtonConvergenceWriter;
+ // typedef NewtonConvergenceWriter<TypeTag> StaggeredNewtonConvergenceWriter;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
diff --git a/dumux/implicit/staggered/propertydefaults.hh b/dumux/implicit/staggered/propertydefaults.hh
index b1dc5044a9..b3d7c27bd0 100644
--- a/dumux/implicit/staggered/propertydefaults.hh
+++ b/dumux/implicit/staggered/propertydefaults.hh
@@ -30,6 +30,7 @@
#include <dumux/implicit/propertydefaults.hh>
#include <dumux/discretization/staggered/fvgridgeometry.hh>
#include <dumux/discretization/staggered/fvelementgeometry.hh>
+// #include <dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh>
#include <dumux/implicit/staggered/properties.hh>
#include <dumux/discretization/methods.hh>
@@ -52,14 +53,13 @@
#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dumux/common/intersectionmapper.hh>
-#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 "gridvariables.hh"
namespace Dumux {
@@ -82,6 +82,9 @@ SET_TYPE_PROP(StaggeredModel, FVGridGeometry, StaggeredFVGridGeometry<TypeTag, G
//! Set the default for the local finite volume geometry
SET_TYPE_PROP(StaggeredModel, FVElementGeometry, StaggeredFVElementGeometry<TypeTag, GET_PROP_VALUE(TypeTag, EnableFVGridGeometryCache)>);
+// SET_TYPE_PROP(StaggeredModel, FVElementGeometry, Dumux::CCTpfaFVElementGeometry<TypeTag, GET_PROP_VALUE(TypeTag, EnableFVGridGeometryCache)>);
+//TODO: try to use CC
+SET_TYPE_PROP(StaggeredModel, GridVariables, StaggeredGridVariables<TypeTag>);
//! The sub control volume
SET_PROP(StaggeredModel, SubControlVolume)
@@ -108,12 +111,6 @@ SET_TYPE_PROP(StaggeredModel, GlobalVolumeVariables, Dumux::StaggeredGlobalVolum
//! The global flux variables cache vector class
SET_TYPE_PROP(StaggeredModel, GlobalFluxVariablesCache, Dumux::StaggeredGlobalFluxVariablesCache<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFluxVariablesCache)>);
-//! The local jacobian operator
-SET_TYPE_PROP(StaggeredModel, LocalJacobian, Dumux::StaggeredLocalJacobian<TypeTag>);
-
-//! Assembler for the global jacobian matrix
-SET_TYPE_PROP(StaggeredModel, JacobianAssembler, Dumux::StaggeredAssembler<TypeTag>);
-
//! The local flux variables cache vector class
SET_TYPE_PROP(StaggeredModel, ElementFluxVariablesCache, Dumux::StaggeredElementFluxVariablesCache<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFluxVariablesCache)>);
@@ -244,7 +241,7 @@ public:
};
//! use the plain newton convergence writer by default
-SET_TYPE_PROP(StaggeredModel, NewtonConvergenceWriter, StaggeredNewtonConvergenceWriter<TypeTag>);
+// SET_TYPE_PROP(StaggeredModel, NewtonConvergenceWriter, StaggeredNewtonConvergenceWriter<TypeTag>);
//! Write separate vtp files for face variables by default
SET_BOOL_PROP(StaggeredModel, VtkWriteFaceData, true);
diff --git a/dumux/io/pointcloudvtkwriter.hh b/dumux/io/pointcloudvtkwriter.hh
index c008658e0a..8b9e546f51 100644
--- a/dumux/io/pointcloudvtkwriter.hh
+++ b/dumux/io/pointcloudvtkwriter.hh
@@ -25,7 +25,7 @@
#include <dune/common/fvector.hh>
-#include <dumux/io/vtkoutputmodulebase.hh>
+#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dune/grid/io/file/vtk/common.hh>
#include <dune/common/path.hh>
diff --git a/dumux/io/staggeredvtkoutputmodule.hh b/dumux/io/staggeredvtkoutputmodule.hh
index 4d8fcec2d2..d0b40646c8 100644
--- a/dumux/io/staggeredvtkoutputmodule.hh
+++ b/dumux/io/staggeredvtkoutputmodule.hh
@@ -25,7 +25,7 @@
#include <dune/common/fvector.hh>
-#include <dumux/io/vtkoutputmodulebase.hh>
+#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/pointcloudvtkwriter.hh>
#include <dumux/io/vtksequencewriter.hh>
@@ -44,10 +44,10 @@ namespace Dumux
* Specialization for staggered grids with dofs on faces.
*/
template<typename TypeTag>
-class StaggeredVtkOutputModule : public VtkOutputModuleBase<TypeTag>
+class StaggeredVtkOutputModule : public VtkOutputModule<TypeTag>
{
- friend class VtkOutputModuleBase<TypeTag>;
- using ParentType = VtkOutputModuleBase<TypeTag>;
+ friend class VtkOutputModule<TypeTag>;
+ using ParentType = VtkOutputModule<TypeTag>;
using Implementation = typename GET_PROP_TYPE(TypeTag, VtkOutputModule);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
@@ -68,215 +68,222 @@ class StaggeredVtkOutputModule : public VtkOutputModuleBase<TypeTag>
using Positions = std::vector<Scalar>;
using Data = std::vector<std::vector<Scalar>>;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+
public:
StaggeredVtkOutputModule(const Problem& problem,
- Dune::VTK::DataMode dm = Dune::VTK::conforming) : ParentType(problem, dm),
- faceWriter_(std::make_shared<PointCloudVtkWriter<Scalar, dim>>(coordinates_)),
- sequenceWriter_(faceWriter_, problem.name() + "-face", "","",
- problem.gridView().comm().rank(),
- problem.gridView().comm().size() )
-
- {
- writeFaceVars_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, WriteFaceData);
- coordinatesInitialized_ = false;
- }
-
- //////////////////////////////////////////////////////////////////////////////////////////////
- //! Methods to conveniently add primary and secondary variables upon problem initialization
- //! Do not call these methods after initialization
- //////////////////////////////////////////////////////////////////////////////////////////////
-
-
- //! Output a scalar field
- //! \param name The name of the vtk field
- //! \returns A reference to the resized scalar field to be filled with the actual data
- std::vector<Scalar>& createFaceScalarField(const std::string& name)
- {
- faceScalarFields_.emplace_back(std::make_pair(std::vector<Scalar>(this->problem().model().numFaceDofs()), name));
- return faceScalarFields_.back().first;
- }
-
- //! Output a vector field
- //! \param name The name of the vtk field
- //! \returns A reference to the resized vector field to be filled with the actual data
- std::vector<GlobalPosition>& createFaceVectorField(const std::string& name)
- {
- faceVectorFields_.emplace_back(std::make_pair(std::vector<GlobalPosition>(this->problem().model().numFaceDofs()), name));
- return faceVectorFields_.back().first;
- }
-
- //! Output a scalar primary variable
- //! \param name The name of the vtk field
- //! \param pvIdx The index in the primary variables vector
- void addFacePrimaryVariable(const std::string& name, unsigned int pvIdx)
- {
- priVarScalarDataInfo_.push_back(PriVarScalarDataInfo{pvIdx, name});
- }
-
- //! Output a vector primary variable
- //! \param name The name of the vtk field
- //! \param pvIndices A vector of indices in the primary variables vector to group for vector visualization
- void addFacePrimaryVariable(const std::string& name, std::vector<unsigned int> pvIndices)
- {
- assert(pvIndices.size() < 4 && "Vtk doesn't support vector dimensions greater than 3!");
- priVarVectorDataInfo_.push_back(PriVarVectorDataInfo{pvIndices, name});
- }
-
- void write(double time, Dune::VTK::OutputType type = Dune::VTK::ascii)
- {
- ParentType::write(time, type);
- if(writeFaceVars_)
- getFaceDataAndWrite_();
- }
-
-protected:
-
- /*!
- * \brief Returns the number of cell center dofs
- */
- unsigned int numDofs_() const
- {
- return this->problem().model().numCellCenterDofs();
- }
-
- /*!
- * \brief Returns priVar data from dofs not on the face.
- *
- * \param dofIdxGlobal The global dof index
- * \param pvIdx The primary variable index
- */
- auto getPriVarData_(const std::size_t dofIdxGlobal, const std::size_t pvIdx)
- {
- return this->problem().model().curSol()[cellCenterIdx][dofIdxGlobal][pvIdx];
- }
-
- std::vector<PriVarScalarDataInfo> priVarScalarDataInfo_;
- std::vector<PriVarVectorDataInfo> priVarVectorDataInfo_;
- std::vector<PriVarScalarDataInfo> secondVarScalarDataInfo_;
- std::vector<PriVarVectorDataInfo> secondVarVectorDataInfo_;
-
-private:
-
- void updateCoordinates_()
- {
- std::cout << "updating coordinates" << std::endl;
- coordinates_.resize(this->problem().model().numFaceDofs());
- for(auto&& facet : facets(this->problem().gridView()))
- {
- const int dofIdxGlobal = this->problem().gridView().indexSet().index(facet);
- coordinates_[dofIdxGlobal] = facet.geometry().center();
- }
- coordinatesInitialized_ = true;
- }
-
- /*!
- * \brief Gathers all face-related data and invokes the face vtk-writer using these data.
- */
- void getFaceDataAndWrite_()
- {
- const int numPoints = this->problem().model().numFaceDofs();
-
- // make sure not to iterate over the same dofs twice
- std::vector<bool> dofVisited(numPoints, false);
-
- // get fields for all primary coordinates and variables
- if(!coordinatesInitialized_)
- updateCoordinates_();
-
- Data priVarScalarData(priVarScalarDataInfo_.size(), std::vector<Scalar>(numPoints));
-
- Data priVarVectorData(priVarVectorDataInfo_.size());
- for (std::size_t i = 0; i < priVarVectorDataInfo_.size(); ++i)
- priVarVectorData[i].resize(numPoints*priVarVectorDataInfo_[i].pvIdx.size());
-
- for(auto&& element : elements(this->problem().gridView()))
- {
- auto fvGeometry = localView(this->problem().model().fvGridGeometry());
- fvGeometry.bindElement(element);
- for(auto && scvf : scvfs(fvGeometry))
- {
- if(dofVisited[scvf.dofIndex()])
- continue;
-
- asImp_().getPrivarScalarData_(priVarScalarData, scvf);
- asImp_().getPrivarVectorData_(priVarVectorData, scvf);
-
- dofVisited[scvf.dofIndex()] = true;
- }
- }
-
- // transfer priVar scalar data to writer
- for(int i = 0; i < priVarScalarDataInfo_.size(); ++i)
- faceWriter_->addPointData(priVarScalarData[i], priVarScalarDataInfo_[i].name);
-
- // transfer priVar vector data to writer
- for(int i = 0; i < priVarVectorDataInfo_.size(); ++i)
- faceWriter_->addPointData(priVarVectorData[i], priVarVectorDataInfo_[i].name, priVarVectorDataInfo_[i].pvIdx.size());
-
- // transfer custom scalar data to writer
- for(auto&& scalarField : faceScalarFields_)
- faceWriter_->addPointData(scalarField.first, scalarField.second);
-
- // transfer custom vector data to writer
- for(auto&& vectorField : faceVectorFields_)
- faceWriter_->addPointData(vectorField.first, vectorField.second, 3);
-
- sequenceWriter_.write(this->problem().timeManager().time() + 1.0);
- faceScalarFields_.clear();
- faceVectorFields_.clear();
- }
-
-
- /*!
- * \brief Retrives scalar-valued data from the face.
- *
- * \param priVarScalarData Container to store the data
- * \param face The face
- */
- template<class Face>
- void getPrivarScalarData_(Data& priVarScalarData, const Face& face)
- {
- const int dofIdxGlobal = face.dofIndex();
- for(int pvIdx = 0; pvIdx < priVarScalarDataInfo_.size(); ++pvIdx)
- priVarScalarData[pvIdx][dofIdxGlobal] = this->problem().model().curSol()[faceIdx][dofIdxGlobal][pvIdx];
- }
-
- /*!
- * \brief Retrives vector-valued data from the face.
- *
- * \param priVarVectorData Container to store the data
- * \param face The face
- */
- template<class Face>
- void getPrivarVectorData_(Data& priVarVectorData, const Face& face)
- {
- const int dofIdxGlobal = face.dofIndex();
- for (int i = 0; i < priVarVectorDataInfo_.size(); ++i)
- for (int j = 0; j < priVarVectorDataInfo_[i].pvIdx.size(); ++j)
- priVarVectorData[i][dofIdxGlobal*priVarVectorDataInfo_[i].pvIdx.size() + j]
- = this->problem().model().curSol()[faceIdx][dofIdxGlobal][priVarVectorDataInfo_[i].pvIdx[j]];
- }
-
- std::shared_ptr<PointCloudVtkWriter<Scalar, dim>> faceWriter_;
-
- VTKSequenceWriter<PointCloudVtkWriter<Scalar, dim>> sequenceWriter_;
-
- bool writeFaceVars_;
-
- std::vector<GlobalPosition> coordinates_;
- bool coordinatesInitialized_;
-
- std::list<std::pair<std::vector<Scalar>, std::string>> faceScalarFields_;
- std::list<std::pair<std::vector<GlobalPosition>, std::string>> faceVectorFields_;
-
- //! Returns the implementation of the problem (i.e. static polymorphism)
- Implementation &asImp_()
- { return *static_cast<Implementation *>(this); }
-
- //! \copydoc asImp_()
- const Implementation &asImp_() const
- { return *static_cast<const Implementation *>(this); }
+ const FVGridGeometry& fvGridGeometry,
+ const GridVariables& gridVariables,
+ const SolutionVector& sol,
+ const std::string& name,
+ bool verbose = true,
+ Dune::VTK::DataMode dm = Dune::VTK::conforming)
+ : ParentType(problem, fvGridGeometry, gridVariables, sol, name, verbose, dm)
+
+{}
+// {
+// writeFaceVars_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, WriteFaceData);
+// coordinatesInitialized_ = false;
+// }
+//
+// //////////////////////////////////////////////////////////////////////////////////////////////
+// //! Methods to conveniently add primary and secondary variables upon problem initialization
+// //! Do not call these methods after initialization
+// //////////////////////////////////////////////////////////////////////////////////////////////
+//
+//
+// //! Output a scalar field
+// //! \param name The name of the vtk field
+// //! \returns A reference to the resized scalar field to be filled with the actual data
+// std::vector<Scalar>& createFaceScalarField(const std::string& name)
+// {
+// faceScalarFields_.emplace_back(std::make_pair(std::vector<Scalar>(this->problem().model().numFaceDofs()), name));
+// return faceScalarFields_.back().first;
+// }
+//
+// //! Output a vector field
+// //! \param name The name of the vtk field
+// //! \returns A reference to the resized vector field to be filled with the actual data
+// std::vector<GlobalPosition>& createFaceVectorField(const std::string& name)
+// {
+// faceVectorFields_.emplace_back(std::make_pair(std::vector<GlobalPosition>(this->problem().model().numFaceDofs()), name));
+// return faceVectorFields_.back().first;
+// }
+//
+// //! Output a scalar primary variable
+// //! \param name The name of the vtk field
+// //! \param pvIdx The index in the primary variables vector
+// void addFacePrimaryVariable(const std::string& name, unsigned int pvIdx)
+// {
+// priVarScalarDataInfo_.push_back(PriVarScalarDataInfo{pvIdx, name});
+// }
+//
+// //! Output a vector primary variable
+// //! \param name The name of the vtk field
+// //! \param pvIndices A vector of indices in the primary variables vector to group for vector visualization
+// void addFacePrimaryVariable(const std::string& name, std::vector<unsigned int> pvIndices)
+// {
+// assert(pvIndices.size() < 4 && "Vtk doesn't support vector dimensions greater than 3!");
+// priVarVectorDataInfo_.push_back(PriVarVectorDataInfo{pvIndices, name});
+// }
+//
+// void write(double time, Dune::VTK::OutputType type = Dune::VTK::ascii)
+// {
+// ParentType::write(time, type);
+// if(writeFaceVars_)
+// getFaceDataAndWrite_();
+// }
+//
+// protected:
+//
+// /*!
+// * \brief Returns the number of cell center dofs
+// */
+// unsigned int numDofs_() const
+// {
+// return fvGridGeometry
+// }
+//
+// /*!
+// * \brief Returns priVar data from dofs not on the face.
+// *
+// * \param dofIdxGlobal The global dof index
+// * \param pvIdx The primary variable index
+// */
+// auto getPriVarData_(const std::size_t dofIdxGlobal, const std::size_t pvIdx)
+// {
+// return this->problem().model().curSol()[cellCenterIdx][dofIdxGlobal][pvIdx];
+// }
+//
+// std::vector<PriVarScalarDataInfo> priVarScalarDataInfo_;
+// std::vector<PriVarVectorDataInfo> priVarVectorDataInfo_;
+// std::vector<PriVarScalarDataInfo> secondVarScalarDataInfo_;
+// std::vector<PriVarVectorDataInfo> secondVarVectorDataInfo_;
+//
+// private:
+//
+// void updateCoordinates_()
+// {
+// std::cout << "updating coordinates" << std::endl;
+// coordinates_.resize(this->problem().model().numFaceDofs());
+// for(auto&& facet : facets(this->problem().gridView()))
+// {
+// const int dofIdxGlobal = this->problem().gridView().indexSet().index(facet);
+// coordinates_[dofIdxGlobal] = facet.geometry().center();
+// }
+// coordinatesInitialized_ = true;
+// }
+//
+// /*!
+// * \brief Gathers all face-related data and invokes the face vtk-writer using these data.
+// */
+// void getFaceDataAndWrite_()
+// {
+// const int numPoints = this->problem().model().numFaceDofs();
+//
+// // make sure not to iterate over the same dofs twice
+// std::vector<bool> dofVisited(numPoints, false);
+//
+// // get fields for all primary coordinates and variables
+// if(!coordinatesInitialized_)
+// updateCoordinates_();
+//
+// Data priVarScalarData(priVarScalarDataInfo_.size(), std::vector<Scalar>(numPoints));
+//
+// Data priVarVectorData(priVarVectorDataInfo_.size());
+// for (std::size_t i = 0; i < priVarVectorDataInfo_.size(); ++i)
+// priVarVectorData[i].resize(numPoints*priVarVectorDataInfo_[i].pvIdx.size());
+//
+// for(auto&& element : elements(this->problem().gridView()))
+// {
+// auto fvGeometry = localView(this->problem().model().fvGridGeometry());
+// fvGeometry.bindElement(element);
+// for(auto && scvf : scvfs(fvGeometry))
+// {
+// if(dofVisited[scvf.dofIndex()])
+// continue;
+//
+// asImp_().getPrivarScalarData_(priVarScalarData, scvf);
+// asImp_().getPrivarVectorData_(priVarVectorData, scvf);
+//
+// dofVisited[scvf.dofIndex()] = true;
+// }
+// }
+//
+// // transfer priVar scalar data to writer
+// for(int i = 0; i < priVarScalarDataInfo_.size(); ++i)
+// faceWriter_->addPointData(priVarScalarData[i], priVarScalarDataInfo_[i].name);
+//
+// // transfer priVar vector data to writer
+// for(int i = 0; i < priVarVectorDataInfo_.size(); ++i)
+// faceWriter_->addPointData(priVarVectorData[i], priVarVectorDataInfo_[i].name, priVarVectorDataInfo_[i].pvIdx.size());
+//
+// // transfer custom scalar data to writer
+// for(auto&& scalarField : faceScalarFields_)
+// faceWriter_->addPointData(scalarField.first, scalarField.second);
+//
+// // transfer custom vector data to writer
+// for(auto&& vectorField : faceVectorFields_)
+// faceWriter_->addPointData(vectorField.first, vectorField.second, 3);
+//
+// sequenceWriter_.write(this->problem().timeManager().time() + 1.0);
+// faceScalarFields_.clear();
+// faceVectorFields_.clear();
+// }
+//
+//
+// /*!
+// * \brief Retrives scalar-valued data from the face.
+// *
+// * \param priVarScalarData Container to store the data
+// * \param face The face
+// */
+// template<class Face>
+// void getPrivarScalarData_(Data& priVarScalarData, const Face& face)
+// {
+// const int dofIdxGlobal = face.dofIndex();
+// for(int pvIdx = 0; pvIdx < priVarScalarDataInfo_.size(); ++pvIdx)
+// priVarScalarData[pvIdx][dofIdxGlobal] = this->problem().model().curSol()[faceIdx][dofIdxGlobal][pvIdx];
+// }
+//
+// /*!
+// * \brief Retrives vector-valued data from the face.
+// *
+// * \param priVarVectorData Container to store the data
+// * \param face The face
+// */
+// template<class Face>
+// void getPrivarVectorData_(Data& priVarVectorData, const Face& face)
+// {
+// const int dofIdxGlobal = face.dofIndex();
+// for (int i = 0; i < priVarVectorDataInfo_.size(); ++i)
+// for (int j = 0; j < priVarVectorDataInfo_[i].pvIdx.size(); ++j)
+// priVarVectorData[i][dofIdxGlobal*priVarVectorDataInfo_[i].pvIdx.size() + j]
+// = this->problem().model().curSol()[faceIdx][dofIdxGlobal][priVarVectorDataInfo_[i].pvIdx[j]];
+// }
+//
+// std::shared_ptr<PointCloudVtkWriter<Scalar, dim>> faceWriter_;
+//
+// VTKSequenceWriter<PointCloudVtkWriter<Scalar, dim>> sequenceWriter_;
+//
+// bool writeFaceVars_;
+//
+// std::vector<GlobalPosition> coordinates_;
+// bool coordinatesInitialized_;
+//
+// std::list<std::pair<std::vector<Scalar>, std::string>> faceScalarFields_;
+// std::list<std::pair<std::vector<GlobalPosition>, std::string>> faceVectorFields_;
+//
+// //! Returns the implementation of the problem (i.e. static polymorphism)
+// Implementation &asImp_()
+// { return *static_cast<Implementation *>(this); }
+//
+// //! \copydoc asImp_()
+// const Implementation &asImp_() const
+// { return *static_cast<const Implementation *>(this); }
};
} // end namespace Dumux
diff --git a/test/freeflow/staggered/doneatestproblem.hh b/test/freeflow/staggered/doneatestproblem.hh
index ba22ea2b50..5aaff6433a 100644
--- a/test/freeflow/staggered/doneatestproblem.hh
+++ b/test/freeflow/staggered/doneatestproblem.hh
@@ -26,13 +26,14 @@
#ifndef DUMUX_DONEA_TEST_PROBLEM_HH
#define DUMUX_DONEA_TEST_PROBLEM_HH
+#include <dumux/freeflow/staggered/problem.hh>
#include <dumux/implicit/staggered/properties.hh>
-#include <dumux/freeflow/staggered/model.hh>
-#include <dumux/implicit/problem.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/fluidsystems/liquidphase.hh>
#include <dumux/material/components/constant.hh>
+#include <dumux/freeflow/staggered/propertydefaults.hh>
+
namespace Dumux
{
@@ -87,13 +88,13 @@ SET_BOOL_PROP(DoneaTestProblem, EnableInertiaTerms, false);
template <class TypeTag>
class DoneaTestProblem : public NavierStokesProblem<TypeTag>
{
- typedef NavierStokesProblem<TypeTag> ParentType;
+ using ParentType = NavierStokesProblem<TypeTag>;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
// copy some indices for convenience
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum {
// Grid and world dimension
dim = GridView::dimension,
@@ -109,16 +110,17 @@ class DoneaTestProblem : public NavierStokesProblem<TypeTag>
velocityYIdx = Indices::velocityYIdx
};
- typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager);
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::Intersection Intersection;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Intersection = typename GridView::Intersection;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
@@ -132,18 +134,12 @@ class DoneaTestProblem : public NavierStokesProblem<TypeTag>
typename DofTypeIndices::FaceIdx faceIdx;
public:
- DoneaTestProblem(TimeManager &timeManager, const GridView &gridView)
- : ParentType(timeManager, gridView), eps_(1e-6)
+ DoneaTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry), eps_(1e-6)
{
- name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag,
- std::string,
- Problem,
- Name);
-
- printL2Error_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag,
- bool,
- Problem,
- PrintL2Error);
+ name_ = getParam<std::string>("Problem.Name");
+
+ printL2Error_ = getParam<bool>("Problem.PrintL2Error");
}
/*!
@@ -278,8 +274,8 @@ public:
InitialValues initialAtPos(const GlobalPosition& globalPos) const
{
InitialValues values;
- values[pressureIdx] = 0.0;
- values[velocityXIdx] = 0.0;
+ values[pressureIdx] = 123.0;
+ values[velocityXIdx] = 44.0;
values[velocityYIdx] = 0.0;
return values;
diff --git a/test/freeflow/staggered/test_donea.cc b/test/freeflow/staggered/test_donea.cc
index c471bfe00f..a0a1d70449 100644
--- a/test/freeflow/staggered/test_donea.cc
+++ b/test/freeflow/staggered/test_donea.cc
@@ -21,9 +21,37 @@
*
* \brief Test for the staggered grid Navier-Stokes model (Kovasznay 1947)
*/
-#include <config.h>
+
+ #include <config.h>
+
+ #include <ctime>
+ #include <iostream>
+
+ #include <dune/common/parallel/mpihelper.hh>
+ #include <dune/common/timer.hh>
+ #include <dune/grid/io/file/dgfparser/dgfexception.hh>
+ #include <dune/grid/io/file/vtk.hh>
+ #include <dune/istl/io.hh>
+
+
#include "doneatestproblem.hh"
-#include <dumux/common/start.hh>
+
+#include <dumux/common/propertysystem.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/nonlinear/newtonmethod.hh>
+#include <dumux/nonlinear/newtoncontroller.hh>
+
+#include <dumux/assembly/staggeredfvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/methods.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
/*!
* \brief Provides an interface for customizing error messages associated with
@@ -59,6 +87,85 @@ void usage(const char *progName, const std::string &errorMsg)
int main(int argc, char** argv)
{
- typedef TTAG(DoneaTestProblem) ProblemTypeTag;
- return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(DoneaTestProblem);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv, usage);
+
+ // try to create a grid (from the given grid file or the input file)
+ using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+ GridCreator::makeGrid(Parameters::getTree());
+ GridCreator::loadBalance();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = GridCreator::grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // the solution vector
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+ const auto numDofsCellCenter = leafGridView.size(0);
+ const auto numDofsFace = leafGridView.size(1);
+ SolutionVector x;
+ x[cellCenterIdx].resize(numDofsCellCenter);
+ x[faceIdx].resize(numDofsFace);
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // intialize the vtk output module
+ using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+ VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+ VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+ vtkWriter.write(0.0);
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+ assembler->setLinearSystem();
}
diff --git a/test/freeflow/staggered/test_donea.input b/test/freeflow/staggered/test_donea.input
index 59941d9810..3ade63c047 100644
--- a/test/freeflow/staggered/test_donea.input
+++ b/test/freeflow/staggered/test_donea.input
@@ -1,4 +1,4 @@
-[TimeManager]
+[TimeLoop]
DtInitial = 1 # [s]
TEnd = 2 # [s]
@@ -11,7 +11,14 @@ Name = test_donea # name passed to the output routines
LiquidDensity = 1
EnableGravity = false
PrintL2Error = false
+LiquidKinematicViscosity = 1
+
+[Component]
+Name = fu
[ Newton ]
MaxSteps = 10
MaxRelativeShift = 1e-5
+
+[Vtk]
+AddVelocity = 1
--
GitLab