diff --git a/dumux/discretization/staggered/globalvolumevariables.hh b/dumux/discretization/staggered/globalvolumevariables.hh
index cea72bcffc194021be12fe54072cf643eccc8ea2..fd5d2492ce1eae6fc0a3fc0588feb16b77d93e9f 100644
--- a/dumux/discretization/staggered/globalvolumevariables.hh
+++ b/dumux/discretization/staggered/globalvolumevariables.hh
@@ -57,6 +57,10 @@ class StaggeredGlobalVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/true>
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using IndexType = typename GridView::IndexSet::IndexType;
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
static const int dim = GridView::dimension;
using Element = typename GridView::template Codim<0>::Entity;
@@ -75,7 +79,7 @@ public:
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
- volumeVariables_[scv.index()].update(sol[scv.dofIndex()], problem, element, scv);
+ volumeVariables_[scv.index()].update(sol[cellCenterIdx][scv.dofIndex()], problem, element, scv);
// handle the boundary volume variables
for (auto&& scvf : scvfs(fvGeometry))
diff --git a/dumux/freeflow/staggered/model.hh b/dumux/freeflow/staggered/model.hh
index bae02ce4593f2cb83a42950724429d1270353dd3..1df03e8fb914fd25aab32bafcce1620b83c0be28 100644
--- a/dumux/freeflow/staggered/model.hh
+++ b/dumux/freeflow/staggered/model.hh
@@ -57,8 +57,11 @@ template<class TypeTag >
class NavierStokesModel : public GET_PROP_TYPE(TypeTag, BaseModel)
{
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, GlobalFVGeometry) GlobalFVGeometry;
// typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
enum { dim = GridView::dimension };
@@ -69,7 +72,14 @@ class NavierStokesModel : public GET_PROP_TYPE(TypeTag, BaseModel)
using StencilsVector = typename GET_PROP_TYPE(TypeTag, StencilsVector);
using Element = typename GridView::template Codim<0>::Entity;
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
public:
+
+
+
/*!
* \brief \copybrief Dumux::ImplicitModel::addOutputVtkFields
*
@@ -132,61 +142,6 @@ public:
writer.attachCellData(*rank, "process rank");
}
- /*!
- * \brief Returns the number of global degrees of freedoms (DOFs)
- */
- size_t numDofs() const
- {
- return this->gridView_().size(0) + this->gridView_().size(1);
- }
-
- /*!
- * \brief Returns the number of cell center degrees of freedoms (DOFs)
- */
- size_t numCellCenterDofs() const
- {
- return this->gridView_().size(0);
- }
-
- /*!
- * \brief Returns the number of cell center degrees of freedoms (DOFs)
- */
- size_t numFaceDofs() const
- {
- return this->gridView_().size(1);
- }
-
- /*!
- * \brief Returns the size of a complete face dof stencil
- */
- size_t fullFaceToCellCenterStencilSize(const int idx) const
- {
- return this->stencilsVector_.fullFaceToCellCenterStencilSize(idx);
- }
-
- /*!
- * \brief Returns the size of a complete face dof stencil
- */
- size_t fullfaceToFaceStencilSize(const int idx) const
- {
- return this->stencilsVector_.fullfaceToFaceStencilSize(idx);
- }
-
- /*!
- * \brief Returns a unique pointer to the complete face dof stencils which is used once for setting up the global matrix and deleted afterwards
- */
- auto getFullFaceToCellCenterStencilsPtr()
- {
- return this->stencilsVector_.getFullFaceToCellCenterStencilsPtr();
- }
-
- /*!
- * \brief Returns a unique pointer to the complete face dof stencils which is used once for setting up the global matrix and deleted afterwards
- */
- auto getFullfaceToFaceStencilsPtr()
- {
- return this->stencilsVector_.getFullfaceToFaceStencilsPtr();
- }
};
}
diff --git a/dumux/implicit/staggered/assembler.hh b/dumux/implicit/staggered/assembler.hh
index e1f528ed885eac6cb7e0c48a1447be76cf2011da..e7a334f6e95bca9e3be940425e77dd105f268bbf 100644
--- a/dumux/implicit/staggered/assembler.hh
+++ b/dumux/implicit/staggered/assembler.hh
@@ -55,6 +55,37 @@ class StaggeredAssembler : public ImplicitAssembler<TypeTag>
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_()
{
diff --git a/dumux/implicit/staggered/localjacobian.hh b/dumux/implicit/staggered/localjacobian.hh
index 0059b955c2fa60c7694d319bf26c58cb806334fe..939de60938f69bce58fbdee27d152b3342c05965 100644
--- a/dumux/implicit/staggered/localjacobian.hh
+++ b/dumux/implicit/staggered/localjacobian.hh
@@ -90,6 +90,10 @@ class StaggeredLocalJacobian : public ImplicitLocalJacobian<TypeTag>
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
using AssemblyMap = std::vector<std::vector<std::vector<IndexType>>>;
public:
@@ -145,14 +149,14 @@ public:
if (isGhost)
{
this->residual_ = 0.0;
- residual[globalI_] = 0.0;
+ residual[cellCenterIdx][globalI_] = 0.0;
}
else
{
this->localResidual().eval(element, fvGeometry, prevElemVolVars, curElemVolVars, elemBcTypes, elemFluxVarsCache);
this->residual_ = this->localResidual().residual();
// store residual in global container as well
- residual[globalI_] = this->localResidual().residual(0);
+ residual[cellCenterIdx][globalI_] = this->localResidual().residual(0);
}
this->model_().updatePVWeights(fvGeometry);
diff --git a/dumux/implicit/staggered/model.hh b/dumux/implicit/staggered/model.hh
new file mode 100644
index 0000000000000000000000000000000000000000..afdc5d04667efaae1c62d41df85d6e8dd8ccf6ff
--- /dev/null
+++ b/dumux/implicit/staggered/model.hh
@@ -0,0 +1,496 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Base class for all models which use the one-phase,
+ * fully implicit model.
+ * Adaption of the fully implicit scheme to the one-phase flow model.
+ */
+
+#ifndef DUMUX_STAGGERED_BASEMODEL_HH
+#define DUMUX_STAGGERED_BASEMODEL_HH
+
+// #include <dumux/porousmediumflow/implicit/velocityoutput.hh>
+#include <dumux/implicit/model.hh>
+#include "properties.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup NavierStokesModel
+ * \brief A single-phase, isothermal flow model using the fully implicit scheme.
+ *
+ * Single-phase, isothermal flow model, which uses a standard Darcy approach as the
+ * equation for the conservation of momentum:
+ * \f[
+ v = - \frac{\textbf K}{\mu}
+ \left(\textbf{grad}\, p - \varrho {\textbf g} \right)
+ * \f]
+ *
+ * and solves the mass continuity equation:
+ * \f[
+ \phi \frac{\partial \varrho}{\partial t} + \text{div} \left\lbrace
+ - \varrho \frac{\textbf K}{\mu} \left( \textbf{grad}\, p -\varrho {\textbf g} \right) \right\rbrace = q,
+ * \f]
+ * All equations are discretized using a vertex-centered finite volume (box)
+ * or cell-centered finite volume scheme as spatial
+ * and the implicit Euler method as time discretization.
+ * The model supports compressible as well as incompressible fluids.
+ */
+template<class TypeTag >
+class StaggeredBaseModel : public ImplicitModel<TypeTag>
+{
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, GlobalFVGeometry) GlobalFVGeometry;
+// typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+ typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ enum { dim = GridView::dimension };
+ enum { dimWorld = GridView::dimensionworld };
+
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
+ using StencilsVector = typename GET_PROP_TYPE(TypeTag, StencilsVector);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Implementation = typename GET_PROP_TYPE(TypeTag, Model);
+ using NewtonMethod = typename GET_PROP_TYPE(TypeTag, NewtonMethod);
+ using NewtonController = typename GET_PROP_TYPE(TypeTag, NewtonController);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+ static_assert(!isBox, "must not be box!");
+
+public:
+
+ /*!
+ * \brief Apply the initial conditions to the model.
+ *
+ * \param problem The object representing the problem which needs to
+ * be simulated.
+ */
+ void init(Problem &problem)
+ {
+ this->problemPtr_ = &problem;
+
+ this->updateBoundaryIndices_();
+
+ this->globalFvGeometryPtr_ = std::make_shared<GlobalFVGeometry>(problem.gridView());
+ this->globalFvGeometryPtr_->update(problem);
+
+ this->uCur_[cellCenterIdx].resize(asImp_().numCellCenterDofs());
+ this->uCur_[faceIdx].resize(asImp_().numFaceDofs());
+
+ // apply initial solution
+ // for compositional models initial the phase presence herein
+ asImp_().applyInitialSolution_();
+
+ // resize and update the volVars with the initial solution
+ this->curGlobalVolVars_.update(problem, this->curSol());
+
+ // update stencils
+ this->stencilsVector_.update(problem);
+
+ // update the flux variables caches
+ this->globalfluxVarsCache_.update(problem);
+
+ // initialize assembler and create matrix
+ this->localJacobian_.init(problem);
+ this->jacAsm_ = std::make_shared<JacobianAssembler>();
+ this->jacAsm_->init(problem);
+
+ // also set the solution of the "previous" time step to the
+ // initial solution.
+ this->uPrev_ = this->uCur_;
+ this->prevGlobalVolVars_ = this->curGlobalVolVars_;
+ }
+
+ /*!
+ * \brief Try to progress the model to the next timestep.
+ *
+ * \param solver The non-linear solver
+ * \param controller The controller which specifies the behaviour
+ * of the non-linear solver
+ */
+ bool update(NewtonMethod &solver,
+ NewtonController &controller)
+ {
+#if HAVE_VALGRIND
+ for (size_t i = 0; i < this->curSol()[cellCenterIdx].size(); ++i)
+ Valgrind::CheckDefined(this->curSol()[cellCenterIdx][i]);
+ for (size_t i = 0; i < this->curSol()[faceIdx].size(); ++i)
+ Valgrind::CheckDefined(this->curSol()[faceIdx][i]);
+#endif // HAVE_VALGRIND
+
+ asImp_().updateBegin();
+
+ int converged = solver.execute(controller);
+
+ if (this->gridView_().comm().size() > 1)
+ {
+ converged = this->gridView_().comm().min(converged);
+ }
+ if (converged) {
+ asImp_().updateSuccessful();
+ }
+ else
+ asImp_().updateFailed();
+
+#if HAVE_VALGRIND
+ for (size_t i = 0; i < this->curSol()[cellCenterIdx].size(); ++i)
+ Valgrind::CheckDefined(this->curSol()[cellCenterIdx][i]);
+ for (size_t i = 0; i < this->curSol()[faceIdx].size(); ++i)
+ Valgrind::CheckDefined(this->curSol()[faceIdx][i]);
+#endif // HAVE_VALGRIND
+
+ return converged;
+ }
+
+ /*!
+ * \brief Applies the initial solution for all vertices of the grid.
+ *
+ * \todo the initial condition needs to be unique for
+ * each vertex. we should think about the API...
+ */
+ void applyInitialSolution_()
+ {
+ // first set the whole domain to zero
+ this->uCur_ = Scalar(0.0);
+
+ // iterate through leaf grid and evaluate initial
+ // condition at the center of each sub control volume
+ for (const auto& element : elements(this->gridView_()))
+ {
+ // deal with the current element, bind FVGeometry to it
+ auto fvGeometry = localView(this->globalFvGeometry());
+ 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 = this->problem_().initial(scv);
+
+ auto dofIdxGlobal = scv.dofIndex();
+ this->uCur_[cellCenterIdx][dofIdxGlobal] += initPriVars;
+ }
+ }
+ }
+
+ /*!
+ * \brief Write the current solution for a vertex to a restart
+ * file.
+ *
+ * \param outstream The stream into which the vertex data should
+ * be serialized to
+ * \param entity The entity which's data should be
+ * serialized, i.e. a vertex for the box method
+ * and an element for the cell-centered method
+ */
+ template <class Entity>
+ void serializeEntity(std::ostream &outstream,
+ const Entity &entity)
+ {
+ DUNE_THROW(Dune::NotImplemented, "deserializeEntity() not implemented yet");
+// int dofIdxGlobal = dofMapper().index(entity);
+// int dofIdxGlobal = dofMapper().index(entity);
+//
+// // write phase state
+// if (!outstream.good()) {
+// DUNE_THROW(Dune::IOError,
+// "Could not serialize vertex "
+// << dofIdxGlobal);
+// }
+//
+// for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+// outstream << curSol()[dofIdxGlobal][eqIdx] << " ";
+// }
+ }
+
+ /*!
+ * \brief Reads the current solution variables for a vertex from a
+ * restart file.
+ *
+ * \param instream The stream from which the vertex data should
+ * be deserialized from
+ * \param entity The entity which's data should be
+ * serialized, i.e. a vertex for the box method
+ * and an element for the cell-centered method
+ */
+ template <class Entity>
+ void deserializeEntity(std::istream &instream,
+ const Entity &entity)
+ {
+ DUNE_THROW(Dune::NotImplemented, "deserializeEntity() not implemented yet");
+// int dofIdxGlobal = dofMapper().index(entity);
+//
+// for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+// if (!instream.good())
+// DUNE_THROW(Dune::IOError,
+// "Could not deserialize vertex "
+// << dofIdxGlobal);
+// instream >> curSol()[dofIdxGlobal][eqIdx];
+// }
+ }
+
+ /*!
+ * \brief \copybrief Dumux::ImplicitModel::addOutputVtkFields
+ *
+ * Specialization for the NavierStokesModel, adding the pressure and
+ * the process rank to the VTK writer.
+ */
+ template<class MultiWriter>
+ void addOutputVtkFields(const SolutionVector &sol,
+ MultiWriter &writer)
+ {
+ // typedef Dune::BlockVector<Dune::FieldVector<double, dimWorld> > VectorField;
+
+ // create the required scalar fields
+ unsigned numDofs = this->numDofs();
+ auto *p = writer.allocateManagedBuffer(numDofs);
+ // VectorField *velocity = writer.template allocateManagedBuffer<double, dimWorld>(numDofs);
+ // ImplicitVelocityOutput<TypeTag> velocityOutput(this->problem_());
+
+ // if (velocityOutput.enableOutput())
+ // {
+ // // initialize velocity field
+ // for (unsigned int i = 0; i < numDofs; ++i)
+ // {
+ // (*velocity)[i] = double(0);
+ // }
+ // }
+
+ unsigned numElements = this->gridView_().size(0);
+ auto *rank = writer.allocateManagedBuffer(numElements);
+
+ for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
+ {
+ auto eIdx = this->elementMapper().index(element);
+ (*rank)[eIdx] = this->gridView_().comm().rank();
+
+ // get the local fv geometry
+ auto fvGeometry = localView(this->globalFvGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(this->curGlobalVolVars());
+ elemVolVars.bindElement(element, fvGeometry, this->curSol());
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+ auto dofIdxGlobal = scv.dofIndex();
+
+ (*p)[dofIdxGlobal] = volVars.pressure();
+ }
+
+ // velocity output
+ //velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, element, /*phaseIdx=*/0);
+ }
+
+ writer.attachDofData(*p, "p", isBox);
+ // if (velocityOutput.enableOutput())
+ // {
+ // writer.attachDofData(*velocity, "velocity", isBox, dim);
+ // }
+ writer.attachCellData(*rank, "process rank");
+ }
+
+ /*!
+ * \brief Add the vector fields for analysing the convergence of
+ * the newton method to the a VTK multi writer.
+ *
+ * \tparam MultiWriter The type of the VTK multi writer
+ *
+ * \param writer The VTK multi writer object on which the fields should be added.
+ * \param u The solution function
+ * \param deltaU The delta of the solution function before and after the Newton update
+ */
+ template <class MultiWriter>
+ void addConvergenceVtkFields(MultiWriter &writer,
+ const SolutionVector &u,
+ const SolutionVector &deltaU)
+ {
+// typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
+//
+// SolutionVector residual(u);
+// asImp_().globalResidual(residual, u);
+//
+// // create the required scalar fields
+// unsigned numDofs = asImp_().numDofs();
+//
+// // global defect of the two auxiliary equations
+// ScalarField* def[numEq];
+// ScalarField* delta[numEq];
+// ScalarField* x[numEq];
+// for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+// x[eqIdx] = writer.allocateManagedBuffer(numDofs);
+// delta[eqIdx] = writer.allocateManagedBuffer(numDofs);
+// def[eqIdx] = writer.allocateManagedBuffer(numDofs);
+// }
+//
+// for (unsigned int dofIdxGlobal = 0; dofIdxGlobal < u.size(); dofIdxGlobal++)
+// {
+// for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+// {
+// (*x[eqIdx])[dofIdxGlobal] = u[dofIdxGlobal][eqIdx];
+// (*delta[eqIdx])[dofIdxGlobal] = - deltaU[dofIdxGlobal][eqIdx];
+// (*def[eqIdx])[dofIdxGlobal] = residual[dofIdxGlobal][eqIdx];
+// }
+// }
+//
+// for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+// std::ostringstream oss;
+// oss.str(""); oss << "x_" << eqIdx;
+// if (isBox)
+// writer.attachVertexData(*x[eqIdx], oss.str());
+// else
+// writer.attachCellData(*x[eqIdx], oss.str());
+// oss.str(""); oss << "delta_" << eqIdx;
+// if (isBox)
+// writer.attachVertexData(*delta[eqIdx], oss.str());
+// else
+// writer.attachCellData(*delta[eqIdx], oss.str());
+// oss.str(""); oss << "defect_" << eqIdx;
+// if (isBox)
+// writer.attachVertexData(*def[eqIdx], oss.str());
+// else
+// writer.attachCellData(*def[eqIdx], oss.str());
+// }
+//
+// asImp_().addOutputVtkFields(u, writer);
+ }
+
+ /*!
+ * \brief Compute the global residual for an arbitrary solution
+ * vector.
+ *
+ * \param residual Stores the result
+ * \param u The solution for which the residual ought to be calculated
+ */
+ Scalar globalResidual(SolutionVector &residual,
+ const SolutionVector &u)
+ {
+ SolutionVector tmp(this->curSol());
+ this->curSol() = u;
+ Scalar res = globalResidual(residual);
+ this->curSol() = tmp;
+ return res;
+ }
+
+ /*!
+ * \brief Compute the global residual for the current solution
+ * vector.
+ *
+ * \param residual Stores the result
+ */
+ Scalar globalResidual(SolutionVector &residual)
+ {
+ residual = 0;
+
+ for (const auto& element : elements(this->gridView_())) {
+ this->localResidual().eval(element);
+
+
+ int globalI = this->elementMapper().index(element);
+ residual[cellCenterIdx][globalI] = this->localResidual().residual(0);
+ }
+
+ // calculate the square norm of the residual
+ Scalar result2 = residual.two_norm2();
+ if (this->gridView_().comm().size() > 1)
+ result2 = this->gridView_().comm().sum(result2);
+
+ return std::sqrt(result2);
+ }
+
+ /*!
+ * \brief Returns the number of global degrees of freedoms (DOFs)
+ */
+ size_t numDofs() const
+ {
+ return this->gridView_().size(0) + this->gridView_().size(1);
+ }
+
+ /*!
+ * \brief Returns the number of cell center degrees of freedoms (DOFs)
+ */
+ size_t numCellCenterDofs() const
+ {
+ return this->gridView_().size(0);
+ }
+
+ /*!
+ * \brief Returns the number of cell center degrees of freedoms (DOFs)
+ */
+ size_t numFaceDofs() const
+ {
+ return this->gridView_().size(1);
+ }
+
+ /*!
+ * \brief Returns the size of a complete face dof stencil
+ */
+ size_t fullFaceToCellCenterStencilSize(const int idx) const
+ {
+ return this->stencilsVector_.fullFaceToCellCenterStencilSize(idx);
+ }
+
+ /*!
+ * \brief Returns the size of a complete face dof stencil
+ */
+ size_t fullfaceToFaceStencilSize(const int idx) const
+ {
+ return this->stencilsVector_.fullfaceToFaceStencilSize(idx);
+ }
+
+ /*!
+ * \brief Returns a unique pointer to the complete face dof stencils which is used once for setting up the global matrix and deleted afterwards
+ */
+ auto getFullFaceToCellCenterStencilsPtr()
+ {
+ return this->stencilsVector_.getFullFaceToCellCenterStencilsPtr();
+ }
+
+ /*!
+ * \brief Returns a unique pointer to the complete face dof stencils which is used once for setting up the global matrix and deleted afterwards
+ */
+ auto getFullfaceToFaceStencilsPtr()
+ {
+ return this->stencilsVector_.getFullfaceToFaceStencilsPtr();
+ }
+private:
+
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+
+};
+}
+
+#include "propertydefaults.hh"
+
+#endif
diff --git a/dumux/implicit/staggered/newtoncontroller.hh b/dumux/implicit/staggered/newtoncontroller.hh
index 75d750399c7a41e27059cb089431ae37399582f2..791d8779e59428998480b40c0557ff892ee95b17 100644
--- a/dumux/implicit/staggered/newtoncontroller.hh
+++ b/dumux/implicit/staggered/newtoncontroller.hh
@@ -49,6 +49,10 @@ class StaggeredNewtonController : public NewtonController<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
public:
StaggeredNewtonController(const Problem &problem)
: ParentType(problem)
@@ -119,6 +123,88 @@ public:
}
}
+ /*!
+ * \brief Update the current solution with a delta vector.
+ *
+ * The error estimates required for the newtonConverged() and
+ * newtonProceed() methods should be updated inside this method.
+ *
+ * Different update strategies, such as line search and chopped
+ * updates can be implemented. The default behavior is just to
+ * subtract deltaU from uLastIter, i.e.
+ * \f[ u^{k+1} = u^k - \Delta u^k \f]
+ *
+ * \param uCurrentIter The solution vector after the current iteration
+ * \param uLastIter The solution vector after the last iteration
+ * \param deltaU The delta as calculated from solving the linear
+ * system of equations. This parameter also stores
+ * the updated solution.
+ */
+ void newtonUpdate(SolutionVector &uCurrentIter,
+ const SolutionVector &uLastIter,
+ const SolutionVector &deltaU)
+ {
+ if (this->enableShiftCriterion_)
+ this->newtonUpdateShift(uLastIter, deltaU);
+
+ this->writeConvergence_(uLastIter, deltaU);
+
+ if (this->useLineSearch_)
+ {
+ this->lineSearchUpdate_(uCurrentIter, uLastIter, deltaU);
+ }
+ else {
+ for (unsigned int i = 0; i < uLastIter[cellCenterIdx].size(); ++i) {
+ uCurrentIter[cellCenterIdx][i] = uLastIter[cellCenterIdx][i];
+ uCurrentIter[cellCenterIdx][i] -= deltaU[cellCenterIdx][i];
+ }
+ for (unsigned int i = 0; i < uLastIter[faceIdx].size(); ++i) {
+ uCurrentIter[faceIdx][i] = uLastIter[faceIdx][i];
+ uCurrentIter[faceIdx][i] -= deltaU[faceIdx][i];
+ }
+
+ if (this->enableResidualCriterion_)
+ {
+ SolutionVector tmp(uLastIter);
+ this->reduction_ = this->method().model().globalResidual(tmp, uCurrentIter);
+ this->reduction_ /= this->initialResidual_;
+ }
+ }
+ }
+
+ /*!
+ * \brief Update the maximum relative shift of the solution compared to
+ * the previous iteration.
+ *
+ * \param uLastIter The current iterative solution
+ * \param deltaU The difference between the current and the next solution
+ */
+ void newtonUpdateShift(const SolutionVector &uLastIter,
+ const SolutionVector &deltaU)
+ {
+ this->shift_ = 0;
+
+ for (int i = 0; i < int(uLastIter[cellCenterIdx].size()); ++i) {
+ auto uNewI = uLastIter[cellCenterIdx][i];
+ uNewI -= deltaU[cellCenterIdx][i];
+
+ Scalar shiftAtDof = this->model_().relativeShiftAtDof(uLastIter[cellCenterIdx][i],
+ uNewI);
+ this->shift_ = std::max(this->shift_, shiftAtDof);
+ }
+ for (int i = 0; i < int(uLastIter[faceIdx].size()); ++i) {
+ auto uNewI = uLastIter[faceIdx][i];
+ uNewI -= deltaU[faceIdx][i];
+
+ Scalar shiftAtDof = this->model_().relativeShiftAtDof(uLastIter[faceIdx][i],
+ uNewI);
+ this->shift_ = std::max(this->shift_, shiftAtDof);
+ }
+
+ if (this->gridView_().comm().size() > 1)
+ this->shift_ = this->gridView_().comm().max(this->shift_);
+ }
+
};
}
diff --git a/dumux/implicit/staggered/properties.hh b/dumux/implicit/staggered/properties.hh
index bcf6c90cbd0352beebb197b14ae87243db05f3ca..7710767b869fda0653a36fe0eebc48b25a198a26 100644
--- a/dumux/implicit/staggered/properties.hh
+++ b/dumux/implicit/staggered/properties.hh
@@ -45,6 +45,9 @@ NEW_TYPE_TAG(StaggeredModel, INHERITS_FROM(ImplicitBase));
NEW_PROP_TAG(CellCenterPrimaryVariables); //!< A vector of primary variables for cell center dofs
NEW_PROP_TAG(FacePrimaryVariables); //!< A vector of primary variables for face dofs
+NEW_PROP_TAG(CellCenterSolutionVector); //!< Vector containing all cell centered primary variables
+NEW_PROP_TAG(FaceSolutionVector); //!< Vector containing all face primary variables
+
NEW_PROP_TAG(NumEqCellCenter); //!< Number of equations per cell center dof
NEW_PROP_TAG(NumEqFace); //!< Number of equations per face dof
diff --git a/dumux/implicit/staggered/propertydefaults.hh b/dumux/implicit/staggered/propertydefaults.hh
index 9e10b99ec2e99ac201795de15774053a73a8eafb..315d3a305c48478a184965b5b5f4b01e99d62424 100644
--- a/dumux/implicit/staggered/propertydefaults.hh
+++ b/dumux/implicit/staggered/propertydefaults.hh
@@ -61,6 +61,7 @@
#include "localjacobian.hh"
#include "properties.hh"
#include "newtoncontroller.hh"
+#include "model.hh"
namespace Dumux {
@@ -74,6 +75,10 @@ SET_PROP(StaggeredModel, DiscretizationMethod)
static const DiscretizationMethods value = DiscretizationMethods::Staggered;
};
+
+SET_TYPE_PROP(StaggeredModel, BaseModel, StaggeredBaseModel<TypeTag>);
+
+
//! Set the default for the global finite volume geometry
SET_TYPE_PROP(StaggeredModel, GlobalFVGeometry, StaggeredGlobalFVGeometry<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFVGeometryCache)>);
@@ -146,8 +151,39 @@ SET_TYPE_PROP(StaggeredModel, NewtonController, StaggeredNewtonController<TypeTa
SET_INT_PROP(StaggeredModel, NumEqCellCenter, 1);
SET_INT_PROP(StaggeredModel, NumEqFace, 1);
+//! A vector of primary variables
+SET_TYPE_PROP(StaggeredModel,
+ CellCenterPrimaryVariables,
+ Dune::FieldVector<typename GET_PROP_TYPE(TypeTag, Scalar),
+ GET_PROP_VALUE(TypeTag, NumEqCellCenter)>);
+//! A vector of primary variables
+SET_TYPE_PROP(StaggeredModel,
+ FacePrimaryVariables,
+ Dune::FieldVector<typename GET_PROP_TYPE(TypeTag, Scalar),
+ GET_PROP_VALUE(TypeTag, NumEqFace)>);
+
+//! The type of a solution for the whole grid at a fixed time
+SET_TYPE_PROP(StaggeredModel,
+ CellCenterSolutionVector,
+ Dune::BlockVector<typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables)>);
+
+//! The type of a solution for the whole grid at a fixed time
+SET_TYPE_PROP(StaggeredModel,
+ FaceSolutionVector,
+ Dune::BlockVector<typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables)>);
+
+//! default property value for the solution vector only used for monolithic solver
+SET_PROP(StaggeredModel, SolutionVector)
+{
+private:
+ using CellCenterSolutionVector = typename GET_PROP_TYPE(TypeTag, CellCenterSolutionVector);
+ using FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
+public:
+ typedef typename Dune::MultiTypeBlockVector<CellCenterSolutionVector, FaceSolutionVector> type;
+};
+
//! Set the type of a global jacobian matrix from the solution types
-SET_PROP(StaggeredModel, JacobianMatrix/*New*/)
+SET_PROP(StaggeredModel, JacobianMatrix)
{
private:
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);