diff --git a/dumux/freeflow/staggered/model.hh b/dumux/freeflow/staggered/model.hh
index 08ef270ff2b222e8f21f18f2655fbdeac3503700..bae02ce4593f2cb83a42950724429d1270353dd3 100644
--- a/dumux/freeflow/staggered/model.hh
+++ b/dumux/freeflow/staggered/model.hh
@@ -137,7 +137,23 @@ public:
*/
size_t numDofs() const
{
- return this->gridView_().size(0) + this->gridView_().size(1);;
+ 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);
}
/*!
diff --git a/dumux/implicit/staggered/assembler.hh b/dumux/implicit/staggered/assembler.hh
index 3d9689cbcc76022a1038d15a9a0f5a8ad836086a..e1f528ed885eac6cb7e0c48a1447be76cf2011da 100644
--- a/dumux/implicit/staggered/assembler.hh
+++ b/dumux/implicit/staggered/assembler.hh
@@ -40,11 +40,55 @@ class StaggeredAssembler : public 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;
- void setRowSizes_()
+ 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;
+
+ // 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);
+
+
+ setRowSizes_(A11, A12, A21, A22);
+ addIndices_(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 setRowSizes_(CCToCCMatrixBlock &A11, CCToFaceMatrixBlock &A12,
+ FaceToCCMatrixBlock &A21, FaceToFaceMatrixBlock &A22)
{
for (const auto& element : elements(this->gridView_()))
{
@@ -52,22 +96,27 @@ class StaggeredAssembler : public ImplicitAssembler<TypeTag>
const auto globalI = this->elementMapper_().index(element);
const auto& cellCenterToCellCenterStencil = this->model_().stencils(element).cellCenterToCellCenterStencil();
const auto& cellCenterToFaceStencil = this->model_().stencils(element).cellCenterToFaceStencil();
- const auto size = cellCenterToCellCenterStencil.size() + cellCenterToFaceStencil.size();
- this->matrix().setrowsize(globalI, size);
+ A11.setrowsize(globalI, cellCenterToCellCenterStencil.size());
+ A12.setrowsize(globalI, cellCenterToFaceStencil.size());
}
+ A11.endrowsizes();
+ A12.endrowsizes();
- for(int facetIdx = 0; facetIdx < this->gridView_().size(1); ++facetIdx)
+ for(int faceDofxIdx = 0; faceDofxIdx < this->gridView_().size(1); ++faceDofxIdx)
{
- const int faceDofxIdx = facetIdx + this->gridView_().size(0);
- auto size = this->model_().fullFaceToCellCenterStencilSize(facetIdx);
- size += this->model_().fullfaceToFaceStencilSize(facetIdx);
- this->matrix().setrowsize(faceDofxIdx, size);
+ const auto faceToCellSize = this->model_().fullFaceToCellCenterStencilSize(faceDofxIdx);
+ const auto faceToFaceSize = this->model_().fullfaceToFaceStencilSize(faceDofxIdx);
+
+ A21.setrowsize(faceDofxIdx, faceToCellSize);
+ A22.setrowsize(faceDofxIdx, faceToFaceSize);
}
- this->matrix().endrowsizes();
+ A21.endrowsizes();
+ A22.endrowsizes();
}
- void addIndices_()
+ void addIndices_(CCToCCMatrixBlock &A11, CCToFaceMatrixBlock &A12,
+ FaceToCCMatrixBlock &A21, FaceToFaceMatrixBlock &A22)
{
for (const auto& element : elements(this->gridView_()))
{
@@ -78,36 +127,39 @@ class StaggeredAssembler : public ImplicitAssembler<TypeTag>
for (auto&& globalJ : cellCenterToCellCenterStencil)
- this->matrix().addindex(globalI, globalJ);
+ A11.addindex(globalI, globalJ);
for (auto&& globalJ : cellCenterToFaceStencil)
- this->matrix().addindex(globalI, globalJ);
+ A12.addindex(globalI, globalJ - this->gridView_().size(0));
}
+ A11.endindices();
+ A12.endindices();
auto ptr1 = this->model_().getFullFaceToCellCenterStencilsPtr();
auto ptr2 = this->model_().getFullfaceToFaceStencilsPtr();
- if(ptr1 && ptr2)
- std::cout << "success!!!" << std::endl;
+ assert((ptr1 && ptr2) && "pointers to full face stencils are empty!");
const auto& fullFaceToCellCenterStencils = ptr1.get()[0];
const auto& fullfaceToFaceStencils = ptr2.get()[0];
- int globalI = this->gridView_().size(0);
+ int globalI = 0;
for(const auto& stencil : fullFaceToCellCenterStencils)
{
for(auto&& globalJ : stencil)
- this->matrix().addindex(globalI, globalJ);
+ A21.addindex(globalI, globalJ);
++globalI;
}
- globalI = this->gridView_().size(0);
+
+ globalI = 0;
for(const auto& stencil : fullfaceToFaceStencils)
{
for(auto&& globalJ : stencil)
- this->matrix().addindex(globalI, globalJ);
+ A22.addindex(globalI, globalJ - this->gridView_().size(0));
++globalI;
}
- this->matrix().endindices();
+ A21.endindices();
+ A22.endindices();
}
};
diff --git a/dumux/implicit/staggered/newtoncontroller.hh b/dumux/implicit/staggered/newtoncontroller.hh
new file mode 100644
index 0000000000000000000000000000000000000000..75d750399c7a41e27059cb089431ae37399582f2
--- /dev/null
+++ b/dumux/implicit/staggered/newtoncontroller.hh
@@ -0,0 +1,125 @@
+// -*- 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 A 2p1cni specific controller for the newton solver.
+ *
+ * This controller 'knows' what a 'physically meaningful' solution is
+ * which allows the newton method to abort quicker if the solution is
+ * way out of bounds.
+ */
+#ifndef DUMUX_STAGGERED_NEWTON_CONTROLLER_HH
+#define DUMUX_STAGGERED_NEWTON_CONTROLLER_HH
+
+#include "properties.hh"
+
+#include <dumux/nonlinear/newtoncontroller.hh>
+
+namespace Dumux {
+/*!
+ * \ingroup PNMModel
+ * \brief A PNM specific controller for the newton solver.
+ *
+ * This controller 'knows' what a 'physically meaningful' solution is
+ * which allows the newton method to abort quicker if the solution is
+ * way out of bounds.
+ */
+template <class TypeTag>
+class StaggeredNewtonController : public NewtonController<TypeTag>
+{
+ typedef NewtonController<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+public:
+ StaggeredNewtonController(const Problem &problem)
+ : ParentType(problem)
+ {}
+
+ /*!
+ * \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$.
+ *
+ * Throws Dumux::NumericalProblem if the linear solver didn't
+ * converge.
+ *
+ * \param A The matrix of the linear system of equations
+ * \param x The vector which solves the linear system
+ * \param b The right hand side of the linear system
+ */
+ void newtonSolveLinear(JacobianMatrix &A,
+ SolutionVector &x,
+ SolutionVector &b)
+ {
+ try {
+ if (this->numSteps_ == 0)
+ {
+ Scalar norm2 = b.two_norm2();
+ if (this->gridView_().comm().size() > 1)
+ norm2 = this->gridView_().comm().sum(norm2);
+
+ this->initialResidual_ = std::sqrt(norm2);
+ }
+
+ int converged = 1;//this->linearSolver_.solve(A, x, b); //TODO: fix
+
+ // make sure all processes converged
+ int convergedRemote = converged;
+ if (this->gridView_().comm().size() > 1)
+ convergedRemote = this->gridView_().comm().min(converged);
+
+ if (!converged) {
+ DUNE_THROW(NumericalProblem,
+ "Linear solver did not converge");
+ }
+ else if (!convergedRemote) {
+ DUNE_THROW(NumericalProblem,
+ "Linear solver did not converge on a remote process");
+ }
+ }
+ catch (Dune::MatrixBlockError e) {
+ // make sure all processes converged
+ int converged = 0;
+ if (this->gridView_().comm().size() > 1)
+ converged = this->gridView_().comm().min(converged);
+
+ Dumux::NumericalProblem p;
+ std::string msg;
+ std::ostringstream ms(msg);
+// ms << e.what() << "M=" << A[e.r][e.c];
+// p.message(ms.str());
+// throw p;
+ }
+ catch (const Dune::Exception &e) {
+ // make sure all processes converged
+ int converged = 0;
+ if (this->gridView_().comm().size() > 1)
+ converged = this->gridView_().comm().min(converged);
+
+ Dumux::NumericalProblem p;
+ p.message(e.what());
+ throw p;
+ }
+ }
+
+};
+}
+
+#endif
diff --git a/dumux/implicit/staggered/properties.hh b/dumux/implicit/staggered/properties.hh
index 73510d65ab1ba9252a7b232c957866f335509b60..bcf6c90cbd0352beebb197b14ae87243db05f3ca 100644
--- a/dumux/implicit/staggered/properties.hh
+++ b/dumux/implicit/staggered/properties.hh
@@ -42,6 +42,15 @@ namespace Properties
//! The type tag for models based on staggered scheme
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(NumEqCellCenter); //!< Number of equations per cell center dof
+NEW_PROP_TAG(NumEqFace); //!< Number of equations per face dof
+NEW_PROP_TAG(DofTypeIndices); //!< Indices to choose between cell center and face dofs
+
+
}
}
diff --git a/dumux/implicit/staggered/propertydefaults.hh b/dumux/implicit/staggered/propertydefaults.hh
index ae7e1ab4e2ad95e4970aed84f6610535d6b8f6cd..9e10b99ec2e99ac201795de15774053a73a8eafb 100644
--- a/dumux/implicit/staggered/propertydefaults.hh
+++ b/dumux/implicit/staggered/propertydefaults.hh
@@ -46,12 +46,21 @@
#include <dumux/freeflow/staggered/fluxvariables.hh>
#include <dumux/freeflow/staggered/fluxvariablescache.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/indices.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/multitypeblockvector.hh>
+#include <dune/istl/multitypeblockmatrix.hh>
+
#include "assembler.hh"
#include "localresidual.hh"
#include "localjacobian.hh"
+#include "properties.hh"
+#include "newtoncontroller.hh"
namespace Dumux {
@@ -132,6 +141,53 @@ SET_TYPE_PROP(StaggeredModel, FluxVariables, FreeFlowFluxVariables<TypeTag>);
//! The flux variables cache class, by default the one for porous media
SET_TYPE_PROP(StaggeredModel, FluxVariablesCache, FreeFlowFluxVariablesCache<TypeTag>);
+SET_TYPE_PROP(StaggeredModel, NewtonController, StaggeredNewtonController<TypeTag>);
+
+SET_INT_PROP(StaggeredModel, NumEqCellCenter, 1);
+SET_INT_PROP(StaggeredModel, NumEqFace, 1);
+
+//! Set the type of a global jacobian matrix from the solution types
+SET_PROP(StaggeredModel, JacobianMatrix/*New*/)
+{
+private:
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ enum {
+ numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter),
+ numEqFace = GET_PROP_VALUE(TypeTag, NumEqFace)
+ };
+
+public:
+ // the sub-blocks
+ using MatrixLittleBlockCCToCC = typename Dune::FieldMatrix<Scalar, numEqCellCenter, numEqCellCenter>; // 2x2
+ using MatrixLittleBlockCCToFace = typename Dune::FieldMatrix<Scalar, numEqFace, numEqCellCenter>; // 1x2
+
+ using MatrixLittleBlockFaceToFace = typename Dune::FieldMatrix<Scalar, numEqFace, numEqFace>; // 1x1
+ using MatrixLittleBlockFaceToCC = typename Dune::FieldMatrix<Scalar, numEqCellCenter, numEqFace>; // 2x1
+
+ // the BCRS matrices of the subproblems as big blocks
+ using MatrixBlockCCToCC = typename Dune::BCRSMatrix<MatrixLittleBlockCCToCC>;
+ using MatrixBlockCCToFace = typename Dune::BCRSMatrix<MatrixLittleBlockCCToFace>;
+
+
+ using MatrixBlockFaceToFace = typename Dune::BCRSMatrix<MatrixLittleBlockFaceToFace>;
+ using MatrixBlockFaceToCC = typename Dune::BCRSMatrix<MatrixLittleBlockFaceToCC>;
+
+ // the row types
+ using RowCellCenter = typename Dune::MultiTypeBlockVector<MatrixLittleBlockCCToCC, MatrixLittleBlockCCToFace>;
+ using RowFace = typename Dune::MultiTypeBlockVector<MatrixLittleBlockFaceToCC, MatrixLittleBlockFaceToFace>;
+
+ // the jacobian matrix
+ using type = typename Dune::MultiTypeBlockMatrix<RowCellCenter, RowFace>;
+};
+
+//! Definition of the indices for cell center and face dofs in the global solution vector
+SET_PROP(StaggeredModel, DofTypeIndices)
+{
+ using CellCenterIdx = Dune::index_constant<0>;
+ using FaceIdx = Dune::index_constant<1>;
+};
+
+
} // namespace Properties
} // namespace Dumux