diff --git a/dumux/multidomain/common/multidomainassembler.hh b/dumux/multidomain/common/multidomainassembler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..32f7bec7a8091bb665548e3116039b2d7e91b89e
--- /dev/null
+++ b/dumux/multidomain/common/multidomainassembler.hh
@@ -0,0 +1,248 @@
+// -*- 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/>. *
+ *****************************************************************************/
+#ifndef DUMUX_MULTIDOMAIN_ASSEMBLER_HH
+#define DUMUX_MULTIDOMAIN_ASSEMBLER_HH
+
+#include "multidomainproperties.hh"
+#include "multidomainpropertydefaults.hh"
+#include <dune/pdelab/constraints/constraintsparameters.hh>
+#include <dune/pdelab/multidomain/constraints.hh>
+
+namespace Dumux {
+
+//! Prevents the setting of a dirichlet constraint anywhere
+struct NoDirichletConstraints :
+ public Dune::PDELab::DirichletConstraintsParameters
+{
+ template<typename I>
+ bool isDirichlet(const I& intersection, const Dune::FieldVector<typename I::ctype, I::dimension-1> & coord) const
+ {
+ return false;
+ }
+};
+
+template<class TypeTag>
+class MultiDomainAssembler
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, MDGrid) MDGrid;
+
+ typedef typename GET_PROP_TYPE(TypeTag, SubProblem1TypeTag) SubTypeTag1;
+ typedef typename GET_PROP_TYPE(TypeTag, SubProblem2TypeTag) SubTypeTag2;
+
+ typedef typename GET_PROP_TYPE(SubTypeTag1, Problem) SubProblem1;
+ typedef typename GET_PROP_TYPE(SubTypeTag2, Problem) SubProblem2;
+
+ typedef typename GET_PROP_TYPE(SubTypeTag1, LocalFEMSpace) FEM1;
+ typedef typename GET_PROP_TYPE(SubTypeTag2, LocalFEMSpace) FEM2;
+
+ typedef typename GET_PROP_TYPE(SubTypeTag1, ScalarGridFunctionSpace) ScalarGridFunctionSpace1;
+ typedef typename GET_PROP_TYPE(SubTypeTag2, ScalarGridFunctionSpace) ScalarGridFunctionSpace2;
+
+ typedef typename GET_PROP_TYPE(SubTypeTag1, GridFunctionSpace) GridFunctionSpace1;
+ typedef typename GET_PROP_TYPE(SubTypeTag2, GridFunctionSpace) GridFunctionSpace2;
+
+ typedef typename GET_PROP_TYPE(SubTypeTag1, LocalOperator) LocalOperator1;
+ typedef typename GET_PROP_TYPE(SubTypeTag2, LocalOperator) LocalOperator2;
+
+ typedef typename GET_PROP_TYPE(TypeTag, MDGridFunctionSpace) MDGridFunctionSpace;
+ typedef typename GET_PROP_TYPE(TypeTag, MDCondition) MDCondition;
+ typedef typename GET_PROP_TYPE(TypeTag, MDSubProblem1) MDSubProblem1;
+ typedef typename GET_PROP_TYPE(TypeTag, MDSubProblem2) MDSubProblem2;
+ typedef typename GET_PROP_TYPE(TypeTag, MDCouplingLocalOperator) MDCouplingLocalOperator;
+ typedef typename GET_PROP_TYPE(TypeTag, MDCoupling) MDCoupling;
+ typedef typename GET_PROP_TYPE(TypeTag, MDConstraintsTrafo) MDConstraintsTrafo;
+ typedef typename GET_PROP_TYPE(TypeTag, MDGridOperator) MDGridOperator;
+
+ typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
+
+ typedef typename GET_PROP_TYPE(SubTypeTag1, Constraints) Constraints1;
+ typedef typename GET_PROP_TYPE(SubTypeTag2, Constraints) Constraints2;
+
+ // copying the jacobian assembler is not a good idea
+ MultiDomainAssembler(const MultiDomainAssembler &);
+
+public:
+ MultiDomainAssembler()
+ {
+ globalProblem_ = 0;
+ problem1_= 0;
+ problem2_= 0;
+ }
+
+ ~MultiDomainAssembler()
+ { }
+
+ /*!
+ * \brief docme
+ *
+ * \param problem docme
+ *
+ */
+ void init(Problem& problem)
+ {
+ globalProblem_ = &problem;
+ problem1_ = &globalProblem_->subProblem1();
+ problem2_ = &globalProblem_->subProblem2();
+
+ fem1_ = Dune::make_shared<FEM1>();
+ fem2_ = Dune::make_shared<FEM2>();
+
+ constraints1_ = Dune::make_shared<Constraints1>();
+ constraints2_ = Dune::make_shared<Constraints2>();
+
+ scalarGridFunctionSpace1_ = Dune::make_shared<ScalarGridFunctionSpace1>(problem1_->gridView(),
+ *fem1_,
+ *constraints1_);
+ scalarGridFunctionSpace2_ = Dune::make_shared<ScalarGridFunctionSpace2>(problem2_->gridView(),
+ *fem2_,
+ *constraints2_);
+ // constraints store indices of ghost dofs
+ constraints1_->compute_ghosts(*scalarGridFunctionSpace1_);
+ constraints2_->compute_ghosts(*scalarGridFunctionSpace2_);
+ Valgrind::CheckDefined(*constraints1_);
+ Valgrind::CheckDefined(*constraints2_);
+// std::cerr << __FILE__ << ":" << __LINE__ << "\n";
+
+ gridFunctionSpace1_ = Dune::make_shared<GridFunctionSpace1>(*scalarGridFunctionSpace1_);
+ gridFunctionSpace2_ = Dune::make_shared<GridFunctionSpace2>(*scalarGridFunctionSpace2_);
+
+ mdGridFunctionSpace_ = Dune::make_shared<MDGridFunctionSpace>(globalProblem_->mdGrid(),
+ *gridFunctionSpace1_,
+ *gridFunctionSpace2_);
+
+ localOperator1_ = Dune::make_shared<LocalOperator1>(problem1_->model());
+ localOperator2_ = Dune::make_shared<LocalOperator2>(problem2_->model());
+
+ condition1_ = Dune::make_shared<MDCondition>(0);
+ condition2_ = Dune::make_shared<MDCondition>(1);
+
+ mdSubProblem1_ = Dune::make_shared<MDSubProblem1>(*localOperator1_, *condition1_);
+ mdSubProblem2_ = Dune::make_shared<MDSubProblem2>(*localOperator2_, *condition2_);
+
+ couplingLocalOperator_ = Dune::make_shared<MDCouplingLocalOperator>(*globalProblem_);
+ mdCoupling_ = Dune::make_shared<MDCoupling>(*mdSubProblem1_, *mdSubProblem2_, *couplingLocalOperator_);
+
+ // TODO proper constraints stuff
+ constraintsTrafo_ = Dune::make_shared<MDConstraintsTrafo>();
+
+ NoDirichletConstraints dirichletVal;
+ auto constraints = Dune::PDELab::MultiDomain::constraints<Scalar>(*mdGridFunctionSpace_,
+ Dune::PDELab::MultiDomain::constrainSubProblem(*mdSubProblem1_,dirichletVal),
+ Dune::PDELab::MultiDomain::constrainSubProblem(*mdSubProblem2_,dirichletVal));
+ constraints.assemble(*constraintsTrafo_);
+
+ mdGridOperator_ = Dune::make_shared<MDGridOperator>(*mdGridFunctionSpace_, *mdGridFunctionSpace_,
+ *constraintsTrafo_, *constraintsTrafo_,
+ *mdSubProblem1_, *mdSubProblem2_, *mdCoupling_);
+
+ matrix_ = Dune::make_shared<JacobianMatrix>(*mdGridOperator_);
+ *matrix_ = 0;
+
+ residual_.resize(matrix_->N());
+ }
+
+ /*!
+ * \brief docme
+ */
+ void assemble()
+ {
+// std::cerr << __FILE__ << ":" << __LINE__ << "\n";
+
+ // assemble the matrix
+ *matrix_ = 0;
+
+ residual_ = 0;
+ mdGridOperator_->jacobian(globalProblem_->model().curSol(), *matrix_);
+// printmatrix(std::cout, matrix_->base(), "global stiffness matrix", "row", 11, 3);
+
+ // calculate the global residual
+ residual_ = 0;
+ mdGridOperator_->residual(globalProblem_->model().curSol(), residual_);
+// printvector(std::cout, residual_, "residual", "row", 200, 1, 3);
+ }
+
+ //! return const reference to matrix
+ const JacobianMatrix &matrix() const
+ { return *matrix_; }
+ //! return reference to matrix
+ // This is not very nice, but required for the AMG solver
+ JacobianMatrix &matrix()
+ { return *matrix_; }
+
+ //! return const reference to residual
+ const SolutionVector &residual() const
+ { return residual_; }
+ SolutionVector &residual()
+ { return residual_; }
+
+ //! return the multidomain gridfunctionspace
+ MDGridFunctionSpace &gridFunctionSpace() const
+ { return *mdGridFunctionSpace_; }
+ MDGridFunctionSpace &mdGridFunctionSpace() const
+ { return *mdGridFunctionSpace_; }
+
+ //! return the multidomain constraints trafo
+ MDConstraintsTrafo &constraintsTrafo() const
+ { return *constraintsTrafo_; }
+
+private:
+ Problem *globalProblem_;
+ SubProblem1 *problem1_;
+ SubProblem2 *problem2_;
+
+ Dune::shared_ptr<FEM1> fem1_;
+ Dune::shared_ptr<FEM2> fem2_;
+
+ Dune::shared_ptr<Constraints1> constraints1_;
+ Dune::shared_ptr<Constraints2> constraints2_;
+
+ Dune::shared_ptr<ScalarGridFunctionSpace1> scalarGridFunctionSpace1_;
+ Dune::shared_ptr<ScalarGridFunctionSpace2> scalarGridFunctionSpace2_;
+
+ Dune::shared_ptr<GridFunctionSpace1> gridFunctionSpace1_;
+ Dune::shared_ptr<GridFunctionSpace2> gridFunctionSpace2_;
+ Dune::shared_ptr<MDGridFunctionSpace> mdGridFunctionSpace_;
+
+ Dune::shared_ptr<LocalOperator1> localOperator1_;
+ Dune::shared_ptr<LocalOperator2> localOperator2_;
+
+ Dune::shared_ptr<MDCondition> condition1_;
+ Dune::shared_ptr<MDCondition> condition2_;
+
+ Dune::shared_ptr<MDSubProblem1> mdSubProblem1_;
+ Dune::shared_ptr<MDSubProblem2> mdSubProblem2_;
+
+ Dune::shared_ptr<MDCouplingLocalOperator> couplingLocalOperator_;
+ Dune::shared_ptr<MDCoupling> mdCoupling_;
+
+ Dune::shared_ptr<MDConstraintsTrafo> constraintsTrafo_;
+ Dune::shared_ptr<MDGridOperator> mdGridOperator_;
+
+ Dune::shared_ptr<JacobianMatrix> matrix_;
+
+ SolutionVector residual_;
+};
+
+} // namespace Dumux
+
+#endif
+
diff --git a/dumux/multidomain/common/multidomainconvergencewriter.hh b/dumux/multidomain/common/multidomainconvergencewriter.hh
new file mode 100644
index 0000000000000000000000000000000000000000..75cf4673563e786edcc590530bed011b15a79ba9
--- /dev/null
+++ b/dumux/multidomain/common/multidomainconvergencewriter.hh
@@ -0,0 +1,163 @@
+// -*- 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 Reference implementation of a newton convergence writer for coupled problems.
+*/
+#ifndef DUMUX_COUPLED_NEWTON_CONVERGENCEWRITER_HH
+#define DUMUX_COUPLED_NEWTON_CONVERGENCEWRITER_HH
+
+#include <dune/pdelab/backend/istlsolverbackend.hh>
+#include <dumux/io/vtkmultiwriter.hh>
+#include <dune/grid/multidomaingrid.hh>
+
+#include "coupledcommon.hh"
+#include "couplednewtoncontroller.hh"
+
+/*!
+ * \file
+ * \brief Forward declaration of properties required for the coupled Newton convergence writer
+ */
+namespace Dumux
+{
+namespace Properties
+{
+ // forward declaration
+ NEW_PROP_TAG(SubProblem1TypeTag);
+ NEW_PROP_TAG(SubProblem2TypeTag);
+ NEW_PROP_TAG(Grid);
+ NEW_PROP_TAG(GridView);
+}
+
+
+//! \cond INTERNAL
+/*!
+ * \brief Writes the intermediate solutions during
+ * the Newton scheme
+ */
+template <class TypeTag>
+struct CoupledNewtonConvergenceWriter
+{
+ //typedef typename GET_PROP_TYPE(TypeTag, Grid) HostGrid;
+ //typedef Dune::mdgrid::FewSubDomainsTraits<HostGrid::dimension,4> MDGridTraits;
+ //typedef Dune::MultiDomainGrid<HostGrid, MDGridTraits> MDGrid;
+ //typedef typename MDGrid::LeafGridView MDGridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, NewtonController) NewtonController;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, SubProblem1TypeTag) SubProblem1TypeTag;
+ typedef typename GET_PROP_TYPE(TypeTag, SubProblem2TypeTag) SubProblem2TypeTag;
+
+ typedef typename GET_PROP_TYPE(SubProblem1TypeTag, GridView) GridView1;
+ typedef typename GET_PROP_TYPE(SubProblem2TypeTag, GridView) GridView2;
+
+ typedef typename GET_PROP_TYPE(SubProblem1TypeTag, SolutionVector) SolutionVector1;
+ typedef typename GET_PROP_TYPE(SubProblem2TypeTag, SolutionVector) SolutionVector2;
+
+ typedef Dumux::VtkMultiWriter<GridView1> VtkMultiWriter1;
+ typedef Dumux::VtkMultiWriter<GridView2> VtkMultiWriter2;
+
+ CoupledNewtonConvergenceWriter(NewtonController &ctl)
+ : ctl_(ctl)
+ {
+ timeStepIndex_ = 0;
+ iteration_ = 0;
+ vtkMultiWriter1_ = 0;
+ vtkMultiWriter2_ = 0;
+ }
+
+ ~CoupledNewtonConvergenceWriter()
+ {
+ delete vtkMultiWriter1_;
+ delete vtkMultiWriter2_;
+ };
+
+ void beginTimestep()
+ {
+ ++timeStepIndex_;
+ iteration_ = 0;
+ if (!vtkMultiWriter1_)
+ vtkMultiWriter1_ = new VtkMultiWriter1(problem_().subProblem1().gridView(), "convergence1");
+
+ if (!vtkMultiWriter2_)
+ vtkMultiWriter2_ = new VtkMultiWriter2(problem_().subProblem2().gridView(), "convergence2");
+ };
+
+ void beginIteration(const GridView1 &gridView1,
+ const GridView2 &gridView2)
+ {
+ ++ iteration_;
+ vtkMultiWriter1_->beginWrite(timeStepIndex_ + iteration_ / 100.0);
+ vtkMultiWriter2_->beginWrite(timeStepIndex_ + iteration_ / 100.0);
+ };
+
+ void writeFields(const SolutionVector &uLastIter,
+ const SolutionVector &deltaU)
+ {
+ SolutionVector1 uLastIter1;
+ SolutionVector2 uLastIter2;
+ SolutionVector1 deltaU1;
+ SolutionVector2 deltaU2;
+
+ uLastIter1.resize(ctl_.method().model().subModel1().numDofs());
+ uLastIter2.resize(ctl_.method().model().subModel2().numDofs());
+ deltaU1.resize(ctl_.method().model().subModel1().numDofs());
+ deltaU2.resize(ctl_.method().model().subModel2().numDofs());
+
+ typedef Dumux::CoupledCommon<TypeTag> Common;
+
+ Common::splitSolVector(uLastIter, uLastIter1, uLastIter2);
+ Common::splitSolVector(deltaU, deltaU1, deltaU2);
+
+
+ std::cout << "\n writing convergence file of current Newton iteration \n";
+ ctl_.method().model().subModel1().addConvergenceVtkFields(*vtkMultiWriter1_, uLastIter1, deltaU1);
+ ctl_.method().model().subModel2().addConvergenceVtkFields(*vtkMultiWriter2_, uLastIter2, deltaU2);
+ };
+
+ void endIteration()
+ {
+ vtkMultiWriter1_->endWrite();
+ vtkMultiWriter2_->endWrite();
+ };
+
+ void endTimestep()
+ {
+ ++timeStepIndex_;
+ iteration_ = 0;
+ };
+
+private:
+ const Problem &problem_() const
+ { return ctl_.method().problem(); }
+
+ int timeStepIndex_;
+ int iteration_;
+ VtkMultiWriter1 *vtkMultiWriter1_;
+ VtkMultiWriter2 *vtkMultiWriter2_;
+ NewtonController &ctl_;
+};
+
+} // namespace Dumux
+
+
+#endif
diff --git a/dumux/multidomain/common/multidomainlocaloperator.hh b/dumux/multidomain/common/multidomainlocaloperator.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3662e08082d62ba5b361b25f4133c3ea31b163fb
--- /dev/null
+++ b/dumux/multidomain/common/multidomainlocaloperator.hh
@@ -0,0 +1,119 @@
+// -*- 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/>. *
+ *****************************************************************************/
+#ifndef DUMUX_MULTIDOMAIN_BOX_LOCAL_OPERATOR_HH
+#define DUMUX_MULTIDOMAIN_BOX_LOCAL_OPERATOR_HH
+
+#include<dune/pdelab/localoperator/pattern.hh>
+#include<dune/pdelab/localoperator/flags.hh>
+
+#include <dumux/implicit/box/boxproperties.hh>
+
+namespace Dumux {
+
+namespace PDELab {
+
+template<class TypeTag>
+class MultiDomainBoxLocalOperator
+:
+public Dune::PDELab::FullVolumePattern,
+public Dune::PDELab::LocalOperatorDefaultFlags
+{
+ // copying the local operator for PDELab is not a good idea
+ MultiDomainBoxLocalOperator(const MultiDomainBoxLocalOperator &);
+
+ typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+ typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
+ typedef typename Grid::Traits::template Codim<0>::EntityPointer EntityPointer;
+
+ enum{numEq = GET_PROP_VALUE(TypeTag, NumEq)};
+
+public:
+ // pattern assembly flags
+ enum { doPatternVolume = true };
+
+ // residual assembly flags
+ enum { doAlphaVolume = true };
+
+ MultiDomainBoxLocalOperator(Model &model)
+ : model_(model)
+ {}
+
+ /*!
+ * \brief Volume integral depending on test and ansatz functions
+ *
+ * \param eg docme
+ * \param lfsu docme
+ * \param x docme
+ * \param lfsv docme
+ * \param r docme
+ *
+ */
+ template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
+ void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x,
+ const LFSV& lfsv, R& r) const
+ {
+ typedef typename LFSU::Traits::SizeType size_type;
+
+ const EntityPointer elementPtr = model_.gridView().grid().subDomainEntityPointer(eg.entity());
+ model_.localResidual().eval(*elementPtr);
+
+ int numVertices = x.size()/numEq;
+ for (size_type comp = 0; comp < r.size(); comp++)
+ r.accumulate(lfsv, comp, model_.localResidual().residual(comp%numVertices)[comp/numVertices]);
+ }
+
+ /*!
+ * \brief Jacobian of volume term
+ *
+ * \param eg docme
+ * \param lfsu docme, is basis
+ * \param x docme
+ * \param lfsv docme, is test
+ * \param mat docme
+ *
+ */
+ template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
+ void jacobian_volume (const EG& eg,
+ const LFSU& lfsu,
+ const X& x,
+ const LFSV& lfsv,
+ M& mat) const
+ {
+ typedef typename LFSU::Traits::SizeType size_typeU;
+ typedef typename LFSV::Traits::SizeType size_typeV;
+
+ const EntityPointer elementPtr = model_.gridView().grid().subDomainEntityPointer(eg.entity());
+ model_.localJacobian().assemble(*elementPtr);
+
+ int numVertices = x.size()/numEq;
+ for (size_typeV j=0; j<lfsv.size(); j++) {
+ for (size_typeU i=0; i<lfsu.size(); i++) {
+ mat.accumulate(lfsv, i, lfsu, j, (model_.localJacobian().mat(i%numVertices,j%numVertices))[i/numVertices][j/numVertices]);
+ }
+ }
+ }
+
+private:
+ Model& model_;
+};
+
+} // namespace PDELab
+} // namespace Dumux
+
+#endif
diff --git a/dumux/multidomain/common/multidomainmodel.hh b/dumux/multidomain/common/multidomainmodel.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ad1124dee9e731fe800d70938f2eba660b800f49
--- /dev/null
+++ b/dumux/multidomain/common/multidomainmodel.hh
@@ -0,0 +1,393 @@
+// -*- 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/>. *
+ *****************************************************************************/
+#ifndef DUMUX_COUPLED_MODEL_HH
+#define DUMUX_COUPLED_MODEL_HH
+
+#include "coupledproperties.hh"
+#include "coupledpropertydefaults.hh"
+#include "coupledproblem.hh"
+#include "couplednewtoncontroller.hh"
+//#include "coupledjacobianassembler.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \defgroup ModelCoupling Coupled implicit models
+ */
+
+
+/*!
+ * \ingroup ModelCoupling
+ *
+ * \brief The base class of models which consist of two arbitrary
+ * sub-models which are coupled
+ */
+template<class TypeTag>
+class CoupledModel
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Model) Implementation;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, NewtonMethod) NewtonMethod;
+ typedef typename GET_PROP_TYPE(TypeTag, NewtonController) NewtonController;
+ typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler;
+ typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, SubProblem1TypeTag) SubTypeTag1;
+ typedef typename GET_PROP_TYPE(TypeTag, SubProblem2TypeTag) SubTypeTag2;
+
+ typedef typename GET_PROP_TYPE(SubTypeTag1, Problem) SubProblem1;
+ typedef typename GET_PROP_TYPE(SubTypeTag2, Problem) SubProblem2;
+
+ typedef typename GET_PROP_TYPE(SubTypeTag1, Model) SubModel1;
+ typedef typename GET_PROP_TYPE(SubTypeTag2, Model) SubModel2;
+
+ enum {
+ numEq1 = GET_PROP_VALUE(TypeTag, NumEq1),
+ numEq2 = GET_PROP_VALUE(TypeTag, NumEq2)
+ };
+
+public:
+ /*!
+ * \brief Apply the initial conditions to the model.
+ *
+ * \param problem docme
+ *
+ */
+ void init(Problem &problem)
+ {
+ problem_ = &problem;
+
+ // the two sub models have already been initialized by the
+ // sub-problems!
+ jacAsm_ = new JacobianAssembler();
+ jacAsm_->init(problem);
+
+// uCur_.resize(asImp_().numDofs());
+// uPrev_.resize(asImp_().numDofs());
+
+ uCur_.resize(jacAsm_->residual().size());
+ uPrev_.resize(jacAsm_->residual().size());
+
+ uCur_= 0;
+ uPrev_= 0;
+
+ typedef Dumux::CoupledCommon<TypeTag> Common;
+ Common::spliceSolVectors(subModel1().curSol(),
+ subModel2().curSol(),
+ uCur_);
+ Common::spliceSolVectors(subModel1().prevSol(),
+ subModel2().prevSol(),
+ uPrev_);
+ }
+
+ Scalar globalResidual(const SolutionVector &u, SolutionVector &tmp)
+ {
+ DUNE_THROW(Dune::NotImplemented, "");
+#if 0
+ SolutionVector tmpU(asImp_(), 0.0);
+ tmpU = uCur_;
+ uCur_ = u;
+ localJacobian_.evalGlobalResidual(tmp);
+
+ Scalar result = tmp.two_norm();
+ /*
+ Scalar result = 0;
+ for (int i = 0; i < (*tmp).size(); ++i) {
+ for (int j = 0; j < numEq; ++j)
+ result += std::abs((*tmp)[i][j]);
+ }
+ */
+ uCur_ = tmpU;
+ return result;
+#endif
+ };
+
+ /*!
+ * \brief Reference to the current solution as a block vector.
+ */
+ const SolutionVector &curSol() const
+ { return uCur_; }
+
+ /*!
+ * \brief Reference to the current solution as a block vector.
+ */
+ SolutionVector &curSol()
+ { return uCur_; }
+
+ /*!
+ * \brief Reference to the previous solution as a block vector.
+ */
+ const SolutionVector &prevSol() const
+ { return uPrev_; }
+
+ /*!
+ * \brief Reference to the previous solution as a block vector.
+ */
+ SolutionVector &prevSol()
+ { return uPrev_; }
+
+ /*!
+ * \brief Returns the operator assembler for the global jacobian of
+ * the problem.
+ */
+ JacobianAssembler &jacobianAssembler()
+ { return *jacAsm_; }
+ const JacobianAssembler &jacobianAssembler() const
+ { return *jacAsm_; }
+
+ /*!
+ * \brief A reference to the problem on which the model is applied.
+ */
+ Problem &problem()
+ { return *problem_; }
+ /*!
+ * \copydoc problem()
+ */
+ const Problem &problem() const
+ { return *problem_; }
+
+ /*!
+ * \brief A reference to the problem on which the model is applied.
+ */
+ SubProblem1 &subProblem1()
+ { return problem().subProblem1(); }
+ /*!
+ * \copydoc subProblem1()
+ */
+ const SubProblem1 &subProblem1() const
+ { return problem().subProblem1(); }
+
+ /*!
+ * \brief A reference to the problem on which the model is applied.
+ */
+ SubProblem2 &subProblem2()
+ { return problem().subProblem2(); }
+ /*!
+ * \copydoc subProblem2()
+ */
+ const SubProblem2 &subProblem2() const
+ { return problem().subProblem2(); }
+
+ /*!
+ * \brief A reference to the first sub-problem's model.
+ */
+ SubModel1 &subModel1()
+ { return subProblem1().model(); }
+ /*!
+ * \copydoc subModel1()
+ */
+ const SubModel1 &subModel1() const
+ { return subProblem1().model(); }
+
+ /*!
+ * \brief A reference to the second sub-problem's model.
+ */
+ SubModel2 &subModel2()
+ { return subProblem2().model(); }
+ /*!
+ * \copydoc subModel2()
+ */
+ const SubModel2 &subModel2() const
+ { return subProblem2().model(); }
+
+ /*!
+ * \brief Try to progress the model to the next timestep.
+ *
+ * \param solver docme
+ * \param controller docme
+ *
+ */
+ bool update(NewtonMethod &solver,
+ NewtonController &controller)
+ {
+#if HAVE_VALGRIND
+ for (size_t i = 0; i < curSol().size(); ++i)
+ Valgrind::CheckDefined(curSol()[i]);
+#endif // HAVE_VALGRIND
+
+ asImp_().updateBegin();
+
+ bool converged = solver.execute(controller);
+ if (!converged)
+ asImp_().updateFailed();
+ else
+ asImp_().updateSuccessful();
+
+#if HAVE_VALGRIND
+ for (size_t i = 0; i < curSol().size(); ++i) {
+ Valgrind::CheckDefined(curSol()[i]);
+ }
+#endif // HAVE_VALGRIND
+
+ return converged;
+ }
+
+
+ /*!
+ * \brief Called by the update() method before it tries to
+ * apply the newton method. This is primary a hook
+ * which the actual model can overload.
+ */
+ void updateBegin()
+ {
+ subModel1().updateBegin();
+ subModel2().updateBegin();
+
+ typedef Dumux::CoupledCommon<TypeTag> Common;
+ Common::spliceSolVectors(subModel1().curSol(), subModel2().curSol(), uCur_);
+ }
+
+
+ /*!
+ * \brief Called by the update() method if it was
+ * successful. This is primary a hook which the actual
+ * model can overload.
+ */
+ void updateSuccessful()
+ {
+ subModel1().updateSuccessful();
+ subModel2().updateSuccessful();
+ }
+
+ /*!
+ * \brief Called by the problem if a timeintegration was
+ * successful, post processing of the solution is done and the
+ * result has been written to disk.
+ *
+ * This should perpare the model for the next time integration.
+ * Note, that the advanceTimeLevel() methods of the sub-models
+ * have already been called by the individual sub problems...
+ */
+ void advanceTimeLevel()
+ {
+ typedef Dumux::CoupledCommon<TypeTag> Common;
+ // splice the two sub-vectors together
+ Common::spliceSolVectors(subModel1().curSol(), subModel2().curSol(), uCur_);
+ Common::spliceSolVectors(subModel1().prevSol(), subModel2().prevSol(), uPrev_);
+ };
+
+ /*!
+ * \brief Called by the update() method if a try was ultimately
+ * unsuccessful. This is primary a hook which the
+ * actual model can overload.
+ */
+ void updateFailed()
+ {
+ subModel1().updateFailed();
+ subModel2().updateFailed();
+
+ typedef Dumux::CoupledCommon<TypeTag> Common;
+ // splice the two sub-vectors together
+ Common::spliceSolVectors(subModel1().curSol(), subModel2().curSol(), uCur_);
+ };
+
+ /*!
+ * \brief Called by the update() method if a try was
+ * unsuccessful. This is primary a hook which the
+ * actual model can overload.
+ */
+ void updateFailedTry()
+ {
+ subModel1().updateFailedTry();
+ subModel2().updateFailedTry();
+
+ typedef Dumux::CoupledCommon<TypeTag> Common;
+ // splice the two sub-vectors together
+ Common::spliceSolVectors(subModel1().curSol(), subModel2().curSol(), uCur_);
+ };
+
+ /*!
+ * \brief Calculate the global residual.
+ *
+ * \param globResidual docme
+ *
+ * The global deflection of the mass balance from zero.
+ */
+ void evalGlobalResidual(SolutionVector &globResidual)
+ {
+ DUNE_THROW(Dune::NotImplemented, "");
+ }
+
+ /*!
+ * \brief Serializes the current state of the model.
+ *
+ * \param res docme
+ *
+ */
+ template <class Restarter>
+ void serialize(Restarter &res)
+ {
+ subProblem1().serialize(res);
+ subProblem2().serialize(res);
+ }
+
+ /*!
+ * \brief Deserializes the state of the model.
+ *
+ * \param res docme
+ *
+ */
+ template <class Restarter>
+ void deserialize(Restarter &res)
+ {
+ subProblem1().deserialize(res);
+ subProblem2().deserialize(res);
+ wasRestarted_ = true;
+ }
+
+ /*!
+ * \brief Returns the number of global degrees of freedoms (DOFs)
+ */
+ size_t numDofs() const
+ {
+ return subModel1().numDofs()*numEq1 + subModel2().numDofs()*numEq2;
+ }
+
+ void resetJacobianAssembler()
+ {
+ delete jacAsm_;
+ jacAsm_ = new JacobianAssembler(asImp_(), problem());
+ }
+
+
+protected:
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+
+ // the problem we want to solve. defines the constitutive
+ // relations, material laws, etc.
+ Problem *problem_;
+
+ // the jacobian assembler
+ JacobianAssembler *jacAsm_;
+
+ // cur is the current solution, prev the solution of the previous
+ // time step
+ SolutionVector uCur_;
+ SolutionVector uPrev_;
+
+ bool wasRestarted_;
+};
+}
+
+#endif
diff --git a/dumux/multidomain/common/multidomainnewtoncontroller.hh b/dumux/multidomain/common/multidomainnewtoncontroller.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d9639d630577a0c185e4167bd237da16f31b1d44
--- /dev/null
+++ b/dumux/multidomain/common/multidomainnewtoncontroller.hh
@@ -0,0 +1,603 @@
+// -*- 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 Reference implementation of a newton controller for coupled problems.
+*/
+#ifndef DUMUX_COUPLED_NEWTON_CONTROLLER_HH
+#define DUMUX_COUPLED_NEWTON_CONTROLLER_HH
+
+#include <dumux/common/exceptions.hh>
+
+#include <dumux/linear/linearsolverproperties.hh>
+#include <dumux/linear/boxlinearsolver.hh>
+
+#include <dune/istl/overlappingschwarz.hh>
+#include <dune/istl/schwarz.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/solvers.hh>
+#include "dune/istl/owneroverlapcopy.hh"
+
+#include <dune/istl/io.hh>
+#include <dune/common/mpihelper.hh>
+#include <iostream>
+#include <boost/format.hpp>
+
+#include <dune/pdelab/backend/istlsolverbackend.hh>
+
+#include "couplednewtonconvergencewriter.hh"
+
+/*!
+ * \file
+ * \brief Additional properties required for the coupled Newton controller
+ */
+namespace Dumux
+{
+template <class TypeTag>
+class CoupledNewtonController;
+
+namespace Properties
+{
+
+//! Specifies the implementation of the Newton controller
+NEW_PROP_TAG(NewtonController);
+
+//! Specifies the type of the actual Newton method
+NEW_PROP_TAG(NewtonMethod);
+
+//! specifies whether the convergence rate and the global residual
+//! gets written out to disk for every newton iteration (default is false)
+NEW_PROP_TAG(NewtonWriteConvergence);
+
+/*!
+ * \brief Specifies whether the update should be done using the line search
+ * method instead of the plain Newton method.
+ *
+ * Whether this property has any effect depends on wether the line
+ * search method is implemented for the actual model's Newton
+ * controller's update() method. By default line search is not used.
+ */
+NEW_PROP_TAG(NewtonUseLineSearch);
+
+//! the value for the relative error below which convergence is declared
+NEW_PROP_TAG(NewtonRelTolerance);
+
+/*!
+ * \brief The number of iterations at which the Newton method
+ * should aim at.
+ *
+ * This is used to control the time step size. The heuristic used
+ * is to scale the last time step size by the deviation of the
+ * number of iterations used from the target steps.
+ */
+NEW_PROP_TAG(NewtonTargetSteps);
+
+//! Number of maximum iterations for the Newton method.
+NEW_PROP_TAG(NewtonMaxSteps);
+
+// set default values for Newton
+// they can be overwritten in the parameter file
+SET_INT_PROP(CoupledModel, NewtonTargetSteps, 8);
+SET_INT_PROP(CoupledModel, NewtonMaxSteps, 15);
+SET_SCALAR_PROP(CoupledModel, NewtonRelTolerance, 1e-5);
+SET_BOOL_PROP(CoupledModel, NewtonWriteConvergence, false);
+}
+
+
+/*!
+ * \brief Reference implementation of a newton controller for coupled problems.
+ *
+ * If you want to specialize only some methods but are happy with
+ * the defaults of the reference controller, derive your
+ * controller from this class and simply overload the required
+ * methods.
+ */
+template <class TypeTag>
+class CoupledNewtonController
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, NewtonController) Implementation;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+ typedef typename GET_PROP_TYPE(TypeTag, NewtonMethod) NewtonMethod;
+ typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, SubProblem1TypeTag) SubTypeTag1;
+ typedef typename GET_PROP_TYPE(TypeTag, SubProblem2TypeTag) SubTypeTag2;
+
+ typedef typename GET_PROP_TYPE(SubTypeTag1, GridView) GridView1;
+ typedef typename GET_PROP_TYPE(SubTypeTag2, GridView) GridView2;
+
+ typedef CoupledNewtonConvergenceWriter<TypeTag> ConvergenceWriter;
+ typedef typename GET_PROP_TYPE(TypeTag, LinearSolver) LinearSolver;
+
+public:
+ CoupledNewtonController(const Problem &problem)
+ : endIterMsgStream_(std::ostringstream::out)
+ , linearSolver_(problem)
+ , convergenceWriter_(asImp_())
+ {
+ verbose_ = true;
+ numSteps_ = 0;
+
+ // maximum tolerated deflection between two iterations
+ setRelTolerance(GET_PARAM_FROM_GROUP(TypeTag, Scalar, Newton, RelTolerance));
+ setTargetSteps(GET_PARAM_FROM_GROUP(TypeTag, Scalar, Newton, TargetSteps));
+ setMaxSteps(GET_PARAM_FROM_GROUP(TypeTag, Scalar, Newton, MaxSteps));
+
+// Writes out, where the relative tolerance is defined
+// std::cout << "ParameterNewtonRelTol= " << PROP_DIAGNOSTIC(TypeTag, NewtonRelTolerance) << std::endl;
+ };
+
+ /*!
+ * \brief Destructor
+ */
+ ~CoupledNewtonController()
+ {
+ };
+
+ /*!
+ * \brief Specifies if the newton method ought to be chatty.
+ *
+ * \param val docme
+ *
+ */
+ void setVerbose(bool val)
+ { verbose_ = val; }
+
+ /*!
+ * \brief Returns true if the newton method ought to be chatty.
+ */
+ bool verbose() const
+ { return verbose_ && gridView_().comm().rank() == 0; }
+
+ /*!
+ * \brief Set the maximum acceptable difference for convergence of
+ * any primary variable between two iterations.
+ *
+ * \param tolerance The maximum relative error between two Newton
+ * iterations at which the scheme is considered
+ * finished
+ */
+ void setRelTolerance(Scalar tolerance)
+ { tolerance_ = tolerance; }
+
+ /*!
+ * \brief Set the maximum acceptable residual norm reduction.
+ *
+ * \param tolerance The maximum reduction of the residual norm
+ * at which the scheme is considered finished
+ */
+// void setAbsTolerance(Scalar tolerance)
+// { absoluteTolerance_ = tolerance; }
+
+ /*!
+ * \brief Set the number of iterations at which the Newton method
+ * should aim at.
+ *
+ * This is used to control the time step size. The heuristic used
+ * is to scale the last time step size by the deviation of the
+ * number of iterations used from the target steps.
+ *
+ * \param targetSteps Number of iterations which are considered "optimal"
+ */
+ void setTargetSteps(int targetSteps)
+ { targetSteps_ = targetSteps; }
+
+ /*!
+ * \brief Set the number of iterations after which the Newton
+ * method gives up.
+ */
+ void setMaxSteps(int maxSteps)
+ { maxSteps_ = maxSteps; }
+
+ /*!
+ * \brief Returns true if another iteration should be done.
+ *
+ * \param uCurrentIter The solution of the current newton iteration
+ */
+ bool newtonProceed(SolutionVector &uCurrentIter)
+ {
+ if (numSteps_ < 2)
+ return true; // we always do at least two iterations
+ else
+ if (numSteps_ >= maxSteps_)
+ {
+ // we have exceeded the allowed number of steps. if the
+ // relative error was reduced by a factor of at least 5,
+ // we proceed anyway
+ return error_*4.0 < lastError_ && !asImp_().newtonConverged();
+ }
+ else
+ if (asImp_().newtonConverged())
+ return false; // we are below the desired tolerance
+
+ return true; // do another round
+ }
+
+ /*!
+ * \brief Returns true if the error of the solution is below the
+ * tolerance.
+ */
+ bool newtonConverged() const
+ { return error_ <= tolerance_; }
+
+ /*!
+ * \brief Called before the newton method is applied to an
+ * non-linear system of equations.
+ *
+ * \param method docme
+ * \param uCurrentIter docme
+ *
+ */
+ void newtonBegin(NewtonMethod &method, const SolutionVector &uCurrentIter)
+ {
+ method_ = &method;
+ numSteps_ = 0;
+
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, WriteConvergence))
+ convergenceWriter_.beginTimestep();
+ }
+
+ /*!
+ * \brief Indidicates the beginning of a newton iteration.
+ */
+ void newtonBeginStep()
+ { lastError_ = error_; }
+
+ /*!
+ * \brief Returns the number of steps done since newtonBegin() was
+ * called.
+ */
+ int newtonNumSteps()
+ { return numSteps_; }
+
+
+ /*!
+ * \brief Update the error of the solution compared to the
+ * previous iteration.
+ *
+ * \param uLastIter docme
+ * \param deltaU docme
+ *
+ */
+ void newtonUpdateRelError(const SolutionVector &uLastIter,
+ const SolutionVector &deltaU)
+ {
+ // calculate the relative error as the maximum relative
+ // deflection in any degree of freedom.
+ error_ = 0;
+
+ SolutionVector uNewI = uLastIter;
+ uNewI -= deltaU;
+
+ for (unsigned int i = 0; i < uLastIter.size(); ++i) {
+ for (unsigned int j = 0; j < uLastIter[i].size(); ++j) {
+ Scalar vertexError = std::abs(deltaU[i][j]);
+ vertexError /= std::max<Scalar>(1.0, std::abs(uLastIter[i][j] + uNewI[i][j])/2);
+
+ error_ = std::max(error_, vertexError);
+ }
+ }
+ }
+
+
+ /*!
+ * \brief Solve the linear system of equations \f$ \mathbf{A}x - b
+ * = 0\f$.
+ *
+ * \param A docme
+ * \param x docme
+ * \param b docme
+ *
+ * Throws Dumux::NumericalProblem if the linear solver didn't
+ * converge.
+ */
+ template <class Matrix, class Vector>
+ void newtonSolveLinear(Matrix &A,
+ Vector &x,
+ Vector &b)
+ {
+ // if the deflection of the newton method is large, we do not
+ // need to solve the linear approximation accurately. Assuming
+ // that the initial value for the delta vector u is quite
+ // close to the final value, a reduction of 6 orders of
+ // magnitude in the defect should be sufficient...
+ try {
+ int converged = linearSolver_.solve(A, x, b);
+
+ // make sure all processes converged
+#if HAVE_MPI
+ int convergedSend = 1;
+ MPI_Allreduce(/*sendBuf=*/&convergedSend,
+ /*recvBuf=*/&converged,
+ /*count=*/1,
+ MPI_INT,
+ MPI_MIN,
+ MPI_COMM_WORLD);
+#endif
+ if (!converged) {
+ DUNE_THROW(NumericalProblem,
+ "Linear solver did not converge");
+ }
+ }
+ catch (const Dune::MatrixBlockError &e) {
+ // make sure all processes converged
+#if HAVE_MPI
+ int convergedSend = 0;
+ int converged;
+
+ MPI_Allreduce(/*sendBuf=*/&convergedSend,
+ /*recvBuf=*/&converged,
+ /*count=*/1,
+ MPI_INT,
+ MPI_MIN,
+ MPI_COMM_WORLD);
+#endif
+
+ 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
+#if HAVE_MPI
+ int convergedSend = 0;
+ int converged;
+
+ MPI_Allreduce(/*sendBuf=*/&convergedSend,
+ /*recvBuf=*/&converged,
+ /*count=*/1,
+ MPI_INT,
+ MPI_MIN,
+ MPI_COMM_WORLD);
+#endif
+
+ Dumux::NumericalProblem p;
+ p.message(e.what());
+ throw p;
+ }
+ }
+
+ /*!
+ * \brief Update the current solution function with a delta vector.
+ *
+ * The error estimates required for the newtonConverged() and
+ * newtonProceed() methods should be updated here.
+ *
+ * Different update strategies, such as line search and chopped
+ * updates can be implemented. The default behaviour is just to
+ * subtract deltaU from uLastIter.
+ *
+ * \param uCurrentIter docme
+ * \param uLastIter The solution of 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 (GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, WriteConvergence)) {
+ writeConvergence_(uLastIter, deltaU);
+ }
+
+ newtonUpdateRelError(uLastIter, deltaU);
+
+ uCurrentIter = uLastIter;
+ uCurrentIter -= deltaU;
+ }
+
+ /*!
+ * \brief Indicates that one newton iteration was finished.
+ *
+ * \param uCurrentIter docme
+ * \param uLastIter The solution of the last iteration
+ *
+ */
+ void newtonEndStep(SolutionVector &uCurrentIter, SolutionVector &uLastIter)
+ {
+ typedef Dumux::CoupledCommon<TypeTag> Common;
+
+ Common::splitSolVector(this->model().curSol(),
+ this->model().subModel1().curSol(),
+ this->model().subModel2().curSol());
+
+ ++numSteps_;
+
+ if (verbose())
+ std::cout << "\rNewton iteration " << numSteps_ << " done: "
+ << "error=" << error_ << endIterMsg().str() << "\n";
+ endIterMsgStream_.str("");
+ }
+
+ /*!
+ * \brief Indicates that we're done solving the non-linear system of equations.
+ */
+ void newtonEnd()
+ {
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, WriteConvergence))
+ convergenceWriter_.endTimestep();
+ }
+
+ /*!
+ * \brief Called if the newton method broke down.
+ *
+ * This method is called _after_ newtonEnd()
+ */
+ void newtonFail()
+ {
+ numSteps_ = std::max(maxSteps_, targetSteps_*2);
+ }
+
+ /*!
+ * \brief Called when the newton method was sucessful.
+ *
+ * This method is called _after_ newtonEnd()
+ */
+ void newtonSucceed()
+ {
+ }
+
+ /*!
+ * \brief Suggest a new time stepsize based on the old time step size.
+ *
+ * The default behaviour is to suggest the old time step size
+ * scaled by the ratio between the target iterations and the
+ * iterations required to actually solve the last time step.
+ *
+ * \param oldTimeStep docme
+ *
+ */
+ Scalar suggestTimeStepSize(Scalar oldTimeStep) const
+ {
+ // be agressive reducing the timestep size but
+ // conservative when increasing it. the rationale is
+ // that we want to avoid failing in the next newton
+ // iteration which would require another linearization
+ // of the problem.
+ if (numSteps_ > targetSteps_) {
+ Scalar percent = ((Scalar) numSteps_ - targetSteps_)/targetSteps_;
+ return oldTimeStep/(1.0 + percent);
+ }
+ else {
+ /*Scalar percent = (Scalar(1))/targetSteps_;
+ return oldTimeStep*(1 + percent);
+ */
+ Scalar percent = ((Scalar) targetSteps_ - numSteps_)/targetSteps_;
+ return oldTimeStep*(1.0 + percent/1.2);
+ }
+ }
+
+ /*!
+ * \brief Returns a reference to the current newton method
+ * which is controlled by this controller.
+ */
+ NewtonMethod &method()
+ { return *method_; }
+
+ /*!
+ * \brief Returns a reference to the current newton method
+ * which is controlled by this controller.
+ */
+ const NewtonMethod &method() const
+ { return *method_; }
+
+ /*!
+ * \brief Returns a reference to the current numeric model.
+ */
+ Model &model()
+ { return method_->model(); }
+
+ /*!
+ * \brief Returns a const reference to the current numeric model.
+ */
+ const Model &model() const
+ { return method_->model(); }
+
+ std::ostringstream &endIterMsg()
+ { return endIterMsgStream_; }
+
+ const GridView1 &gridView1() const
+ { return problem_().gridView1(); }
+ const GridView2 &gridView2() const
+ { return problem_().gridView2(); }
+
+ /*!
+ * \brief the convergence writer produces the output
+ *
+ * \param uLastIter The solution of the last iteration
+ * \param deltaU The delta as calculated from solving the linear
+ * system of equations. This parameter also stores
+ * the updated solution.
+ *
+ */
+ void writeConvergence_(const SolutionVector &uLastIter,
+ const SolutionVector &deltaU)
+ {
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, WriteConvergence)) {
+ convergenceWriter_.beginIteration(gridView1_(), gridView2_());
+ convergenceWriter_.writeFields(uLastIter, deltaU);
+ convergenceWriter_.endIteration();
+ };
+ };
+
+ /*!
+ * \brief Returns a copy of the the grid view.
+ */
+ const GridView gridView_() const
+ { return problem_().gridView(); }
+
+ /*!
+ * \brief the subdomain gridviews
+ */
+ const GridView1 gridView1_() const
+ { return problem_().subProblem1().gridView(); }
+ const GridView2 gridView2_() const
+ { return problem_().subProblem2().gridView(); }
+
+ /*!
+ * \brief the coupled problem
+ */
+ Problem &problem_()
+ { return method_->problem(); }
+ const Problem &problem_() const
+ { return method_->problem(); }
+
+ // returns the actual implementation for the controller we do
+ // it this way in order to allow "poor man's virtual methods",
+ // i.e. methods of subclasses which can be called by the base
+ // class.
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+
+ bool verbose_;
+
+ std::ostringstream endIterMsgStream_;
+ NewtonMethod *method_;
+
+ Scalar tolerance_;
+ Scalar error_;
+ Scalar lastError_;
+
+ // optimal number of iterations we want to achieve
+ int targetSteps_;
+ // maximum number of iterations we do before giving up
+ int maxSteps_;
+ // actual number of steps done so far
+ int numSteps_;
+
+ // the linear solver
+ LinearSolver linearSolver_;
+
+ ConvergenceWriter convergenceWriter_;
+};
+
+} // namespace Dumux
+
+#endif
diff --git a/dumux/multidomain/common/multidomainproblem.hh b/dumux/multidomain/common/multidomainproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..424796aec22f0b8cf1ece0c68d22c0fea1a8d4fa
--- /dev/null
+++ b/dumux/multidomain/common/multidomainproblem.hh
@@ -0,0 +1,387 @@
+// -*- 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 problems which involve two sub problems
+ */
+#ifndef DUMUX_COUPLED_PROBLEM_HH
+#define DUMUX_COUPLED_PROBLEM_HH
+
+#include "coupledmodel.hh"
+#include "couplednewtoncontroller.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup ModelCoupling
+ * \brief Base class for problems which involve two sub problems
+ *
+ * \todo Please doc me more!
+ */
+template<class TypeTag>
+class CoupledProblem
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Implementation;
+
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+ typedef typename GET_PROP_TYPE(TypeTag, NewtonMethod) NewtonMethod;
+ typedef typename GET_PROP_TYPE(TypeTag, NewtonController) NewtonController;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+
+ typedef typename GET_PROP_TYPE(TypeTag, SubProblem1TypeTag) SubTypeTag1;
+ typedef typename GET_PROP_TYPE(TypeTag, SubProblem2TypeTag) SubTypeTag2;
+
+ typedef typename GET_PROP_TYPE(SubTypeTag1, Problem) SubProblem1;
+ typedef typename GET_PROP_TYPE(SubTypeTag2, Problem) SubProblem2;
+
+public:
+ CoupledProblem(TimeManager &timeManager)
+ : timeManager_(timeManager),
+ newtonMethod_(asImp_()),
+ newtonCtl_(asImp_())
+ { };
+
+ /*!
+ * \brief Called by the Dumux::TimeManager in order to
+ * initialize the problem and the sub-problems.
+ *
+ * If you overload this method don't forget to call
+ * ParentType::init()
+ */
+ void init()
+ {
+ // initialize the sub-problems
+ asImp_().subProblem1().init();
+ asImp_().subProblem2().init();
+
+ // specify the elements which couple
+ asImp_().addMetaElements();
+
+ // set the initial condition of the model
+ model().init(asImp_());
+
+ // intialize lagrange multipliers
+ asImp_().initMortarElements();
+ }
+
+ /*!
+ * \brief This method writes the complete state of the problem
+ * to the harddisk.
+ *
+ * \param res docme
+ *
+ * The file will start with the prefix returned by the name()
+ * method, has the current time of the simulation clock in it's
+ * name and uses the extension <tt>.drs</tt>. (Dumux ReStart
+ * file.) See Dumux::Restart for details.
+ */
+ template <class Restarter>
+ void serialize(Restarter &res)
+ {
+ }
+
+ /*!
+ * \brief This method restores the complete state of the problem
+ * from disk.
+ *
+ * \param res docme
+ *
+ * It is the inverse of the serialize() method.
+ */
+ template <class Restarter>
+ void deserialize(Restarter &res)
+ {
+ }
+
+ /*!
+ * \name Simulation control
+ */
+ // \{
+
+ /*!
+ * \brief Start the simulation procedure.
+ *
+ * \param dtInitial docme
+ * \param tEnd docme
+ *
+ * This method is usually called by the main() function and simply
+ * uses Dumux::TimeManager::runSimulation() to do the actual
+ * work.
+ */
+ bool simulate(Scalar dtInitial, Scalar tEnd)
+ {
+ // set the initial time step and the time where the simulation ends
+ timeManager_.setEndTime(tEnd);
+ timeManager_.setTimeStepSize(dtInitial);
+ timeManager_.runSimulation(asImp_());
+ return true;
+ };
+
+
+ /*!
+ * \brief Called by the time manager before the time integration. Calls preTimeStep()
+ * of the subproblems.
+ */
+ void preTimeStep()
+ {
+ asImp_().subProblem1().preTimeStep();
+ asImp_().subProblem2().preTimeStep();
+ }
+
+ /*!
+ * \brief Called by Dumux::TimeManager in order to do a time
+ * integration on the model.
+ */
+ void timeIntegration()
+ {
+ // TODO: should be called from the group Implicit
+ const int maxFails =
+ GET_PARAM_FROM_GROUP(TypeTag, int, Newton, MaxTimeStepDivisions);
+ for (int i = 0; i < maxFails; ++i)
+ {
+ if (model_.update(newtonMethod_, newtonCtl_))
+ return;
+
+ // update failed
+ Scalar dt = timeManager().timeStepSize();
+ Scalar nextDt = dt / 2;
+ timeManager().setTimeStepSize(nextDt);
+
+ std::cout << "Newton solver did not converge. Retrying with time step of "
+ << timeManager().timeStepSize() << "sec\n";
+ }
+
+ DUNE_THROW(Dune::MathError,
+ "Newton solver didn't converge after "
+ << maxFails
+ << " timestep divisions. dt="
+ << timeManager().timeStepSize());
+ }
+
+ /*!
+ * \brief Called by the time manager after the time integration to
+ * do some post processing on the solution. Calls postTimeStep()
+ * of the subproblems.
+ */
+ void postTimeStep()
+ {
+ asImp_().subProblem1().postTimeStep();
+ asImp_().subProblem2().postTimeStep();
+ }
+
+ /*!
+ * \brief Called by Dumux::TimeManager whenever a solution for a
+ * timestep has been computed and the simulation time has
+ * been updated.
+ *
+ * \param dt docme
+ *
+ */
+ Scalar nextTimeStepSize(const Scalar dt)
+ {
+ return newtonCtl_.suggestTimeStepSize(dt);
+ };
+
+ /*!
+ * \brief This method is called by the model if the update to the
+ * next time step failed completely.
+ */
+ void updateSuccessful()
+ {
+ model_.updateSuccessful();
+ };
+
+ /*!
+ * \brief Returns true if the current solution should be written to
+ * disk (i.e. as a VTK file)
+ *
+ * The default behaviour is to write out every the solution for
+ * very time step. This file is intented to be overwritten by the
+ * implementation.
+ */
+ bool shouldWriteOutput() const
+ { return true; }
+
+ /*!
+ * \brief Returns true if the current state of the simulation should be written to
+ * disk
+ */
+ bool shouldWriteRestartFile() const
+ { return false; }
+
+ /*!
+ * \brief Called by the time manager after the end of an episode.
+ */
+ void episodeEnd()
+ {
+ std::cerr << "The end of an episode is reached, but the problem "
+ << "does not override the episodeEnd() method. "
+ << "Doing nothing!\n";
+ }
+
+ // \}
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ * It could be either overwritten by the problem files, or simply
+ * declared over the setName() function in the application file.
+ */
+ const char *name() const
+ {
+ return simname_.c_str();
+ }
+
+ /*!
+ * \brief Set the problem name.
+ *
+ * \param newName docme
+ *
+ * This function sets the simulation name, which should be called before
+ * the application porblem is declared! If not, the default name "sim"
+ * will be used.
+ */
+ static void setName(const char *newName)
+ {
+ simname_ = newName;
+ }
+
+ /*!
+ * \brief Returns TimeManager object used by the simulation
+ */
+ TimeManager &timeManager()
+ { return timeManager_; }
+
+ /*!
+ * \copydoc timeManager()
+ */
+ const TimeManager &timeManager() const
+ { return timeManager_; }
+
+ /*!
+ * \brief Returns NewtonControler object used by the simulation
+ */
+ NewtonController &newtonController()
+ { return newtonCtl_; }
+
+ /*!
+ * \brief Returns NewtonControler object used by the simulation
+ */
+ const NewtonController &newtonController() const
+ { return newtonCtl_; }
+
+ /*!
+ * \brief Returns numerical model used for the problem.
+ */
+ Model &model()
+ { return model_; }
+
+ /*!
+ * \copydoc model()
+ */
+ const Model &model() const
+ { return model_; }
+ // \}
+
+ void computeCouplingIndices() const
+ { DUNE_THROW(Dune::NotImplemented, "Coupled problems need to implement the computeCouplingIndices method!"); }
+
+ SubProblem1 &subProblem1()
+ { DUNE_THROW(Dune::NotImplemented, "Coupled problems need to implement the subProblem1 method!"); }
+ const SubProblem1 &subProblem1() const
+ { DUNE_THROW(Dune::NotImplemented, "Coupled problems need to implement the subProblem1 method!"); }
+
+ SubProblem2 &subProblem2()
+ { DUNE_THROW(Dune::NotImplemented, "Coupled problems need to implement the subProblem2 method!"); }
+ const SubProblem2 &subProblem2() const
+ { DUNE_THROW(Dune::NotImplemented, "Coupled problems need to implement the subProblem2 method!"); }
+
+
+ /*!
+ * \brief Called by the time manager after everything which can be
+ * done about the current time step is finished and the
+ * model should be prepared to do the next time integration.
+ */
+ void advanceTimeLevel()
+ {
+ asImp_().subProblem1().advanceTimeLevel();
+ asImp_().subProblem2().advanceTimeLevel();
+
+ model_.advanceTimeLevel();
+ }
+
+ /*!
+ * \brief Write the relevant quantities of the current solution into an VTK output file.
+ */
+ void writeOutput()
+ {
+ // write the current result to disk
+ if (asImp_().shouldWriteOutput()) {
+ asImp_().subProblem1().writeOutput();
+ asImp_().subProblem2().writeOutput();
+ }
+ }
+
+ /*!
+ * \brief Serialize the simulation's state to disk
+ */
+ void serialize()
+ {
+ DUNE_THROW(Dune::NotImplemented,
+ "Dumux::CoupledProblem::serialize");
+ }
+
+protected:
+ void addMetaElements()
+ {}
+
+ void initMortarElements()
+ {}
+
+ //! 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); }
+
+private:
+ // a string for the name of the current simulation, which could be
+ // set by means of an program argument, for example.
+ static std::string simname_;
+
+ TimeManager &timeManager_;
+ NewtonMethod newtonMethod_;
+ NewtonController newtonCtl_;
+ Model model_;
+};
+// definition of the static class member simname_,
+// which is necessary because it is of type string.
+template <class TypeTag>
+std::string CoupledProblem<TypeTag>::simname_="simCoupled"; //initialized with default "sim"
+
+}
+
+#endif
diff --git a/dumux/multidomain/common/multidomainproperties.hh b/dumux/multidomain/common/multidomainproperties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..154ea37252861c1f8b8d183741a025f66f348f36
--- /dev/null
+++ b/dumux/multidomain/common/multidomainproperties.hh
@@ -0,0 +1,117 @@
+// -*- 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/>. *
+ *****************************************************************************/
+#ifndef DUMUX_COUPLED_PROPERTIES_HH
+#define DUMUX_COUPLED_PROPERTIES_HH
+
+#include <dune/istl/bvector.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dumux/nonlinear/newtonmethod.hh>
+
+#include <dumux/implicit/common/implicitproperties.hh>
+#include <dumux/linear/linearsolverproperties.hh>
+
+#include <dumux/common/timemanager.hh>
+
+/*!
+ * \file
+ * \brief Specify properties required for the coupled model
+ */
+namespace Dumux
+{
+namespace Properties
+{
+/*!
+ * \addtogroup ModelCoupling
+ */
+// \{
+
+//////////////////////////////////////////////////////////////////
+// Type tags tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for problems which utilize the coupling approach
+NEW_TYPE_TAG(CoupledModel, INHERITS_FROM(LinearSolverTypeTag, ImplicitBase));
+
+//! The type tag from which sub-problems of coupling models inherit
+NEW_TYPE_TAG(CoupledSubProblem);
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+//! Specifies the type tag of the first sub-problem
+NEW_PROP_TAG(SubProblem1TypeTag);
+
+//! Specifies the type tag of the second sub-problem
+NEW_PROP_TAG(SubProblem2TypeTag);
+
+//! Specifies the type tag of the other sub-problem
+NEW_PROP_TAG(OtherSubProblemTypeTag);
+
+//! Specifies the type tag of coupled problem
+NEW_PROP_TAG(CoupledProblemTypeTag);
+
+//! Specifies the local jacobian of a meta element
+NEW_PROP_TAG(CoupledLocalJacobian);
+
+//! Specifies the jacobian assembler
+NEW_PROP_TAG(JacobianAssembler);
+
+//! Specifies a meta element
+NEW_PROP_TAG(MetaElement);
+
+//! Specifies a list of meta elements
+NEW_PROP_TAG(MetaElementList);
+
+//! Specifies the type of the jacobian matrix as used for the linear
+//! solver
+NEW_PROP_TAG(JacobianMatrix);
+
+//! specifies whether the convergence rate and the global residual
+//! gets written out to disk for every newton iteration
+NEW_PROP_TAG(NewtonWriteConvergence);
+
+//! the maximum allowed number of timestep divisions for the
+//! Newton solver
+NEW_PROP_TAG(NewtonMaxTimeStepDivisions);
+
+//! Specifies the model
+NEW_PROP_TAG(Model);
+
+//! Specifies the time manager
+NEW_PROP_TAG(TimeManager);
+
+//! Specifies the number of equations in the coupled model
+NEW_PROP_TAG(NumEq);
+
+//! Specifies the number of equations in submodel 1
+NEW_PROP_TAG(NumEq1);
+
+//! Specifies the number of equations in submodel 2
+NEW_PROP_TAG(NumEq2);
+
+//! Specifies the fluidsystem that is used in the subdomains
+NEW_PROP_TAG(FluidSystem);
+
+//! Specifies whether the enriched(mortar) coupling is used (set to false by default)
+NEW_PROP_TAG(DoEnrichedCoupling);
+
+}
+}
+#endif
diff --git a/dumux/multidomain/common/multidomainpropertydefaults.hh b/dumux/multidomain/common/multidomainpropertydefaults.hh
new file mode 100644
index 0000000000000000000000000000000000000000..15ea99f7794972acf568c2db33e389d500c932bd
--- /dev/null
+++ b/dumux/multidomain/common/multidomainpropertydefaults.hh
@@ -0,0 +1,163 @@
+// -*- 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/>. *
+ *****************************************************************************/
+#ifndef DUMUX_COUPLED_PROPERTY_DEFAULTS_HH
+#define DUMUX_COUPLED_PROPERTY_DEFAULTS_HH
+
+#include "coupledcommon.hh"
+#include "couplednewtoncontroller.hh"
+
+#include "coupledproperties.hh"
+
+namespace Dumux
+{
+template <class TypeTag> class CoupledModel;
+template <class TypeTag> class CoupledJacobianAssembler;
+template <class TypeTag> class CoupledNewtonController;
+
+namespace Properties
+{
+
+/////////////////////////////////////7
+// Set property values for the coupled model
+SET_BOOL_PROP(CoupledModel, DoEnrichedCoupling, false);
+SET_TYPE_PROP(CoupledModel, Model, Dumux::CoupledModel<TypeTag>);
+SET_TYPE_PROP(CoupledModel, JacobianAssembler, Dumux::CoupledJacobianAssembler<TypeTag>);
+SET_PROP(CoupledModel, SolutionVector)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+public:
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar, numEq> > type;
+};
+
+//! use the plain newton method for the coupled problems by default
+SET_TYPE_PROP(CoupledModel, NewtonMethod, Dumux::NewtonMethod<TypeTag>);
+
+//! use the plain newton controller for coupled problems by default
+SET_TYPE_PROP(CoupledModel, NewtonController, Dumux::CoupledNewtonController<TypeTag>);
+
+//! Set the default type of the time manager for coupled models
+SET_TYPE_PROP(CoupledModel, TimeManager, Dumux::TimeManager<TypeTag>);
+
+
+SET_PROP(CoupledModel, JacobianMatrix)
+{private:
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+public:
+ typedef Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, numEq, numEq> > type;
+};
+
+SET_PROP(CoupledModel, NumEq)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, SubProblem1TypeTag) TypeTag1;
+ typedef typename GET_PROP_TYPE(TypeTag, SubProblem2TypeTag) TypeTag2;
+
+ enum {
+ numEq1 = GET_PROP_VALUE(TypeTag1, NumEq),
+ numEq2 = GET_PROP_VALUE(TypeTag2, NumEq)
+ };
+
+public:
+ static const int value = numEq1; //TODO: why??
+};
+
+SET_PROP(CoupledModel, NumEq1)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, SubProblem1TypeTag) TypeTag1;
+ enum {numEq = GET_PROP_VALUE(TypeTag1, NumEq)};
+public:
+ static const int value = numEq;
+};
+
+SET_PROP(CoupledModel, NumEq2)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, SubProblem2TypeTag) TypeTag2;
+ enum {numEq = GET_PROP_VALUE(TypeTag2, NumEq)};
+public:
+ static const int value = numEq;
+};
+
+SET_PROP(CoupledModel, MetaElementList)
+{private:
+ typedef typename GET_PROP_TYPE(TypeTag, MetaElement) MetaElement;
+public:
+ typedef std::vector<MetaElement*> type;
+};
+
+SET_PROP(CoupledModel, LinearSolver)
+{public:
+ typedef Dumux::BoxBiCGStabILU0Solver<TypeTag> type;
+};
+
+SET_PROP(CoupledModel, LinearSolverResidualReduction)
+{public:
+ static constexpr double value = 1e-6;
+};
+
+//! set the default number of maximum iterations for the linear solver
+SET_PROP(CoupledModel, LinearSolverMaxIterations)
+{public:
+ static constexpr int value = 250;
+};
+
+SET_INT_PROP(CoupledModel, NewtonMaxTimeStepDivisions, 10);
+
+// use the time manager for the coupled problem in the sub problems
+SET_PROP(CoupledSubProblem, TimeManager)
+{private:
+ typedef typename GET_PROP_TYPE(TypeTag, CoupledProblemTypeTag) CoupledTypeTag;
+public:
+ typedef typename GET_PROP_TYPE(CoupledTypeTag, TimeManager) type;
+};
+
+// use the time manager for the coupled problem in the sub problems
+SET_PROP(CoupledSubProblem, ParameterTree)
+{private:
+ typedef typename GET_PROP_TYPE(TypeTag, CoupledProblemTypeTag) CoupledTypeTag;
+ typedef typename GET_PROP(CoupledTypeTag, ParameterTree) ParameterTree;
+public:
+ typedef typename ParameterTree::type type;
+
+ static type &tree()
+ { return ParameterTree::tree(); }
+
+ static type &compileTimeParams()
+ { return ParameterTree::compileTimeParams(); }
+
+ static type &runTimeParams()
+ { return ParameterTree::runTimeParams(); }
+
+ static type &deprecatedRunTimeParams()
+ { return ParameterTree::deprecatedRunTimeParams(); }
+
+ static type &unusedNewRunTimeParams()
+ { return ParameterTree::unusedNewRunTimeParams(); }
+
+};
+
+// \}
+}
+}
+#endif
+
diff --git a/dumux/multidomain/common/pdelablocaloperator.hh b/dumux/multidomain/common/pdelablocaloperator.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1c061e5170f39ac42ded4509a51f1ff375997af0
--- /dev/null
+++ b/dumux/multidomain/common/pdelablocaloperator.hh
@@ -0,0 +1,147 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2010 by Bernd Flemisch *
+ * Institute for Modelling Hydraulic and Environmental Systems *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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 local operator for PDELab which wraps the box models.
+ */
+#ifndef DUMUX_PDELAB_BOX_LOCAL_OPERATOR_HH
+#define DUMUX_PDELAB_BOX_LOCAL_OPERATOR_HH
+
+#if ! HAVE_DUNE_PDELAB
+#error "DUNE-PDELab must be available in order to include this file!"
+#endif
+
+#include<dune/pdelab/localoperator/pattern.hh>
+#include<dune/pdelab/localoperator/flags.hh>
+
+#include <dumux/implicit/box/boxproperties.hh>
+
+namespace Dumux {
+
+namespace PDELab {
+
+/*!
+ * \brief A local operator for PDELab which wraps the box models.
+ */
+template<class TypeTag>
+class BoxLocalOperator
+ :
+// : public Dune::PDELab::NumericalJacobianApplyVolume<BoxLocalOperatorPDELab<TypeTag> >,
+//public Dune::PDELab::NumericalJacobianVolume<BoxLocalOperatorPDELab<TypeTag> >,
+ public Dune::PDELab::FullVolumePattern,
+ public Dune::PDELab::LocalOperatorDefaultFlags
+{
+ // copying the local operator for PDELab is not a good idea
+ BoxLocalOperator(const BoxLocalOperator &);
+
+ typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+ enum{numEq = GET_PROP_VALUE(TypeTag, NumEq)};
+
+public:
+ // pattern assembly flags
+ enum { doPatternVolume = true };
+
+ // residual assembly flags
+ enum { doAlphaVolume = true };
+
+
+ /*!
+ * \param model The physical model for the box scheme.
+ */
+ BoxLocalOperator(Model &model)
+ : model_(model)
+ {}
+
+ /*!
+ * \brief Volume integral depending on test and ansatz functions
+ *
+ * \tparam EG The entity geometry type from PDELab
+ * \tparam LFSU The type of the local function space of the ansatz functions
+ * \tparam X The type of the container for the coefficients for the ansatz functions
+ * \tparam LFSV The type of the local function space of the test functions
+ * \tparam R The range type (usually FieldVector<double>)
+ *
+ * \param eg The entity geometry object
+ * \param lfsu The local function space object of the ansatz functions
+ * \param x The object of the container for the coefficients for the ansatz functions
+ * \param lfsv The local function space object of the test functions
+ * \param r The object storing the volume integral
+ */
+ template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
+ void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x,
+ const LFSV& lfsv, R& r) const
+ {
+ typedef typename LFSU::Traits::SizeType size_type;
+
+ //enum { Green = JacobianAssembler::Green };
+ //if (model_.jacobianAssembler().elementColor(eg.entity()) == Green)
+ // Green elements don't need to be reassembled
+ // return;
+
+ model_.localResidual().eval(eg.entity());
+
+ int numVertices = x.size()/numEq;
+ for (size_type comp = 0; comp < r.size(); comp++)
+ r[comp] = model_.localResidual().residual(comp%numVertices)[comp/numVertices];
+ }
+
+ /*!
+ * \brief Jacobian of volume term
+ *
+ * \tparam EG The entity geometry type from PDELab
+ * \tparam LFSU The type of the local function space of the ansatz functions
+ * \tparam X The type of the container for the coefficients for the ansatz functions
+ * \tparam LFSV The type of the local function space of the test functions
+ * \tparam R The range type (usually FieldVector<double>)
+ *
+ * \param eg The entity geometry object
+ * \param lfsu The local function space object of the ansatz functions
+ * \param x The object of the container for the coefficients for the ansatz functions
+ * \param lfsv The local function space object of the test functions
+ * \param mat The object containing the local jacobian matrix
+ */
+ template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
+ void jacobian_volume (const EG& eg,
+ const LFSU& lfsu,
+ const X& x,
+ const LFSV& lfsv,
+ Dune::PDELab::LocalMatrix<R>& mat) const
+ {
+ typedef typename LFSU::Traits::SizeType size_type;
+
+ model_.localJacobian().assemble(eg.entity());
+
+ int numVertices = x.size()/numEq;
+ for (size_type j=0; j<lfsu.size(); j++) {
+ for (size_type i=0; i<lfsu.size(); i++) {
+ mat(i,j) = (model_.localJacobian().mat(i%numVertices,j%numVertices))[i/numVertices][j/numVertices];
+ }
+ }
+ }
+
+private:
+ Model& model_;
+};
+
+} // namespace PDELab
+} // namespace Dumux
+
+#endif
diff --git a/dumux/multidomain/common/subdomainpropertydefaults.hh b/dumux/multidomain/common/subdomainpropertydefaults.hh
new file mode 100644
index 0000000000000000000000000000000000000000..80037ba46ab4c8d0d75f17210aae330bed25bc8e
--- /dev/null
+++ b/dumux/multidomain/common/subdomainpropertydefaults.hh
@@ -0,0 +1,80 @@
+/*****************************************************************************
+ * 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/>. *
+ *****************************************************************************/
+#ifndef DUMUX_SUBDOMAIN_PROPERTY_DEFAULTS_HH
+#define DUMUX_SUBDOMAIN_PROPERTY_DEFAULTS_HH
+
+#include <dune/grid/multidomaingrid.hh>
+#include <dumux/implicit/pdelab/pdelabadapter.hh>
+#include <dumux/implicit/common/boxcouplinglocalresidual.hh>
+#include <dumux/modelcoupling/common/multidomainboxlocaloperator.hh>
+#include "coupledproperties.hh"
+
+/*!
+ * \file
+ * \brief Specify default properties required in the subdomains of dune-multidomain
+ */
+namespace Dumux
+{
+namespace Properties
+{
+//! The type tag for problems which use dune-multidomain
+NEW_TYPE_TAG(SubDomain, INHERITS_FROM(BoxPDELab, CoupledSubProblem));
+
+//////////////////////////////////////
+// Set property values
+//////////////////////////////////////
+
+// Specifies the grid type for the subdomains
+SET_PROP(SubDomain, Grid)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, CoupledProblemTypeTag) CoupledTypeTag;
+ typedef typename GET_PROP_TYPE(CoupledTypeTag, Grid) HostGrid;
+ typedef typename Dune::mdgrid::FewSubDomainsTraits<HostGrid::dimension,4> MDGridTraits;
+ typedef typename Dune::MultiDomainGrid<HostGrid, MDGridTraits> Grid;
+
+public:
+ typedef typename Grid::SubDomainGrid type;
+};
+
+// Set the default BaseLocalResidual to BoxCouplingLocalResidual
+SET_TYPE_PROP(SubDomain, BaseLocalResidual, BoxCouplingLocalResidual<TypeTag>);
+
+// set the local operator used for submodels
+SET_TYPE_PROP(SubDomain, LocalOperator,
+ Dumux::PDELab::MultiDomainBoxLocalOperator<TypeTag>);
+
+
+// set the grid functions space for the sub-models
+SET_PROP(SubDomain, ScalarGridFunctionSpace)
+{
+ private:
+ typedef typename GET_PROP_TYPE(TypeTag, LocalFEMSpace) FEM;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Constraints) Constraints;
+ enum{numEq = GET_PROP_VALUE(TypeTag, NumEq)};
+ public:
+ typedef Dune::PDELab::GridFunctionSpace<GridView, FEM, Constraints,
+ Dune::PDELab::ISTLVectorBackend<1> > type;
+};
+
+// \}
+
+}
+}
+
+#endif