diff --git a/test/mixeddimension/embedded/1p2c_richards2c/richardstestproblem.hh b/test/mixeddimension/embedded/1p2c_richards2c/richardstestproblem.hh
index f1fdc125c7645f59ee00d8d1068a8046354a665d..b8787e408dea74494316654874aa4731e78507d6 100644
--- a/test/mixeddimension/embedded/1p2c_richards2c/richardstestproblem.hh
+++ b/test/mixeddimension/embedded/1p2c_richards2c/richardstestproblem.hh
@@ -149,7 +149,7 @@ public:
contaminantMoleFraction_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, ContaminantMoleFraction);
// for initial conditions
- const Scalar sw = 0.2; // start with 20% saturation on top
+ const Scalar sw = 0.3; // start with 30% saturation on top
pcTop_ = MaterialLaw::pc(this->spatialParams().materialLawParamsAtPos(this->bBoxMax()), sw);
}
@@ -338,12 +338,13 @@ public:
const auto xTracer = [&,this]()
{
auto contaminationPos = this->bBoxMax()-this->bBoxMin();
- contaminationPos[0] *= 0.2;
- contaminationPos[1] *= 0.5;
- contaminationPos[2] *= 0.2;
+ contaminationPos[0] *= 0.25;
+ contaminationPos[1] *= 0.55;
+ contaminationPos[2] *= 0.25;
contaminationPos += this->bBoxMin();
- if ((globalPos - contaminationPos).two_norm() < 0.1*(this->bBoxMax()-this->bBoxMin()).two_norm() + eps_)
+ static const Scalar extend = 0.15*(this->bBoxMax()[0]-this->bBoxMin()[0]);
+ if ((globalPos - contaminationPos).infinity_norm() < extend + eps_)
return contaminantMoleFraction_;
else
return 0.0;
diff --git a/test/mixeddimension/embedded/1p2c_richards2c/rootsystemtestproblem.hh b/test/mixeddimension/embedded/1p2c_richards2c/rootsystemtestproblem.hh
index 6809d892aa2e710c4cc95fc800d5d9577cea5cb6..40ee2464d658585ee587604111083e18ccffb3ae 100644
--- a/test/mixeddimension/embedded/1p2c_richards2c/rootsystemtestproblem.hh
+++ b/test/mixeddimension/embedded/1p2c_richards2c/rootsystemtestproblem.hh
@@ -26,6 +26,8 @@
#include <cmath>
+#include <dune/grid/common/scsgmapper.hh>
+
#include <dumux/implicit/cellcentered/tpfa/properties.hh>
#include <dumux/porousmediumflow/implicit/problem.hh>
#include <dumux/porousmediumflow/1p2c/implicit/model.hh>
diff --git a/test/mixeddimension/embedded/1p2c_richards2c/rositestproblem.hh b/test/mixeddimension/embedded/1p2c_richards2c/rositestproblem.hh
index 0b449ac614e191c29ddf7d62020d5769e57063d7..8859c6b927086cc8a707426cdafb773e9667432e 100644
--- a/test/mixeddimension/embedded/1p2c_richards2c/rositestproblem.hh
+++ b/test/mixeddimension/embedded/1p2c_richards2c/rositestproblem.hh
@@ -32,6 +32,8 @@
#include <dumux/mixeddimension/embedded/cellcentered/bboxtreecouplingmanagersimple.hh>
#include <dumux/mixeddimension/integrationpointsource.hh>
+#include "schursolver.hh"
+
namespace Dumux
{
template <class TypeTag>
@@ -66,7 +68,7 @@ SET_TYPE_PROP(RootsystemTestCCProblem, PointSourceHelper, Dumux::IntegrationPoin
SET_TYPE_PROP(RichardsTestCCProblem, PointSource, Dumux::IntegrationPointSource<TTAG(RichardsTestCCProblem), unsigned int>);
SET_TYPE_PROP(RichardsTestCCProblem, PointSourceHelper, Dumux::IntegrationPointSourceHelper<TTAG(RichardsTestCCProblem)>);
-SET_TYPE_PROP(RosiTestProblem, LinearSolver, ILU0BiCGSTABBackend<TypeTag>);
+SET_TYPE_PROP(RosiTestProblem, LinearSolver, BlockILU0BiCGSTABSolver<TypeTag>);
}//end namespace properties
diff --git a/test/mixeddimension/embedded/1p2c_richards2c/schursolver.hh b/test/mixeddimension/embedded/1p2c_richards2c/schursolver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a30cd54d62dea640febe17f6eef378816b781d8d
--- /dev/null
+++ b/test/mixeddimension/embedded/1p2c_richards2c/schursolver.hh
@@ -0,0 +1,376 @@
+// -*- 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 Schur complement solver for the embedded problem, the idea is to solve
+ * only linear systems of the size of the much smaller embedded domain by
+ * forming a Schur complement preconditioner. If the embedded system is much
+ * smaller it might be feasible to just use a direct solver.
+ */
+#ifndef DUMUX_EMBEDDED_SCHUR_COMPLEMENT_SOLVER_BACKEND_HH
+#define DUMUX_EMBEDDED_SCHUR_COMPLEMENT_SOLVER_BACKEND_HH
+
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/superlu.hh>
+#include <dune/istl/umfpack.hh>
+
+#include <dumux/common/parameters.hh>
+#include <dumux/common/basicproperties.hh>
+#include <dumux/linear/linearsolverproperties.hh>
+
+#include <dune/istl/paamg/amg.hh>
+#include <dune/istl/paamg/pinfo.hh>
+#include <dumux/linear/amgproperties.hh>
+
+namespace Dumux {
+
+/*!
+ \brief Schur complement linear operator
+ this applies the Schur complement to a vector and can be used in an inverse operator
+*/
+template<class AType, class BType, class CType, class DType, class InvOpDType, class X, class Y>
+class SchurComplement : public Dune::LinearOperator<X,Y>
+{
+public:
+ // export types
+ // typedef DType matrix_type;
+ typedef X domain_type;
+ typedef Y range_type;
+ typedef typename X::field_type field_type;
+
+ //! constructor: just store a reference to a matrix
+ explicit SchurComplement (const AType& A, const BType& B,
+ const CType& C, const DType& D,
+ const std::shared_ptr<InvOpDType>& invD)
+ : A_(A), B_(B), C_(C), D_(D), invD_(invD) {}
+
+ //! apply operator to x: \f$ y = A(x) \f$
+ virtual void apply (const X& x, Y& y) const
+ {
+ // apply C, note x and aTmp1 have different size (Cx)
+ X aTmp1(D_.N());
+ aTmp1 = 0.0;
+ C_.mv(x, aTmp1);
+
+ // apply D^-1 (D^-1Cx)
+ auto aTmp2 = aTmp1;
+ Dune::InverseOperatorResult result;
+ invD_->apply(aTmp2, aTmp1, result);
+
+ // apply B (BD^-1Cx)
+ B_.mv(aTmp2, y);
+
+ y *= -1.0;
+
+ // add Ax (Ax - BD^-1Cx) -> the full Schur complement of S_D := D/M
+ A_.umv(x, y);
+ }
+
+ //! apply operator to x, scale and add: \f$ y = y + \alpha A(x) \f$
+ virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
+ {
+ auto tmp = y;
+ this->apply(x, tmp);
+ y.axpy(alpha, tmp);
+ }
+
+ //! Category of the preconditioner (see SolverCategory::Category)
+ virtual Dune::SolverCategory::Category category() const
+ {
+ return Dune::SolverCategory::sequential;
+ }
+
+ private:
+ const AType& A_;
+ const BType& B_;
+ const CType& C_;
+ const DType& D_;
+ std::shared_ptr<InvOpDType> invD_;
+};
+
+/*! \brief The schur complement preconditioner
+
+ \tparam M The matrix type to operate on
+ \tparam X Type of the update
+ \tparam Y Type of the defect
+ \tparam l The block level to invert. Default is 1
+*/
+template<class TypeTag, class AType, class BType, class CType, class DType, class InvOpDType, class X, class Y>
+class SchurComplementPreconditioner : public Dune::Preconditioner<X, Y>
+{
+public:
+ //! \brief The matrix type the preconditioner is for.
+ // typedef M matrix_type;
+ //! \brief The domain type of the preconditioner.
+ typedef X domain_type;
+ //! \brief The range type of the preconditioner.
+ typedef Y range_type;
+ //! \brief The field type of the preconditioner.
+ typedef typename X::field_type field_type;
+
+ using BulkVector = typename GET_PROP(TypeTag, SolutionVector)::SolutionVectorBulk;
+ using EmbeddedVector = typename GET_PROP(TypeTag, SolutionVector)::SolutionVectorLowDim;
+
+ /*! \brief Approximates the Schur complement
+ */
+ SchurComplementPreconditioner (const AType& A, const BType& B,
+ const CType& C, const DType& D,
+ const std::shared_ptr<InvOpDType>& invD)
+ : A_(A), B_(B), C_(C), D_(D), invD_(invD) {}
+
+
+ virtual void pre (X& x, Y& b)
+ {
+ DUNE_UNUSED_PARAMETER(x);
+ DUNE_UNUSED_PARAMETER(b);
+ }
+
+ /*!
+ \brief Apply the preconditioner.
+
+ \copydoc Preconditioner::apply(X&,const Y&)
+ */
+ virtual void apply (X& x, const Y& b)
+ {
+ using namespace Dune::Indices;
+ const BulkVector& f = b[_0];
+ const EmbeddedVector& g = b[_1];
+
+ static const int verbosity = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, PreconditionerVerbosity);
+ static const int schurMaxIter = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, SchurIterations);
+ static const double schurReduction = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, SchurResidualReduction);
+
+ Dune::Richardson<BulkVector, BulkVector> identity(1.0);
+ SchurComplement<AType, BType, CType, DType, InvOpDType, BulkVector, BulkVector> schurOp(A_, B_, C_, D_, invD_);
+ Dune::BiCGSTABSolver<BulkVector> invOpS(schurOp, identity, schurReduction, schurMaxIter, verbosity);
+
+ // apply S^-1 to solve the bulk block
+ BulkVector fTmp(f);
+ Dune::InverseOperatorResult res;
+ invOpS.apply(x[_0], fTmp, res);
+
+ // prepare rhs for application of D^-1 (g - Ca)
+ EmbeddedVector gTmp(g);
+ C_.mmv(x[_0], gTmp);
+
+ // then solve for embedded block b using the action of D^-1
+ // use direct solver
+ invD_->apply(x[_1], gTmp, res);
+ }
+
+ virtual void post (X& x)
+ { DUNE_UNUSED_PARAMETER(x); }
+
+ //! Category of the preconditioner (see SolverCategory::Category)
+ virtual Dune::SolverCategory::Category category() const
+ { return Dune::SolverCategory::sequential; }
+
+private:
+ const AType& A_;
+ const BType& B_;
+ const CType& C_;
+ const DType& D_;
+ std::shared_ptr<InvOpDType> invD_;
+};
+
+template <class TypeTag>
+class SchurComplementSolver
+{
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+
+public:
+ SchurComplementSolver() = default;
+ SchurComplementSolver(const Problem& problem) {}
+
+ // expects a system as a multi-type matrix
+ // | A B |
+ // | C D |
+
+ // Solve saddle-point problem using a Schur complement based preconditioner
+ template<class Matrix, class Vector>
+ bool solve(const Matrix& M, Vector& x, const Vector& b)
+ {
+ int verbosity = GET_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, Verbosity);
+ static const double reduction = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, ResidualReduction);
+ static const double maxIter = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, MaxIterations);
+ static const int restartGMRes = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, GMResRestart);
+
+ using AType = typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockBulk;
+ using BType = typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockBulkCoupling;
+ using CType = typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockLowDimCoupling;
+ using DType = typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockLowDim;
+
+ using namespace Dune::Indices;
+ const AType& A = M[_0][_0];
+ const BType& B = M[_0][_1];
+ const CType& C = M[_1][_0];
+ const DType& D = M[_1][_1];
+
+ // Solver to compute applications of D^-1 used in the preconditioner
+ // We construct it here so e.g. LU decomposition is only done once
+ using InnerSolver = Dune::UMFPack<DType>;
+ auto invD = std::make_shared<InnerSolver>(D, false);
+
+ SchurComplementPreconditioner<TypeTag, AType, BType, CType, DType, InnerSolver, Vector, Vector> preconditioner(A, B, C, D, invD);
+
+ Dune::MatrixAdapter<Matrix, Vector, Vector> op(M);
+ Dune::RestartedfGMResSolver<Vector> solver(op, preconditioner, reduction, restartGMRes, maxIter, verbosity);
+ auto bTmp(b);
+ solver.apply(x, bTmp, result_);
+
+ return true;
+ }
+
+ const Dune::InverseOperatorResult& result() const
+ {
+ return result_;
+ }
+
+private:
+ Dune::InverseOperatorResult result_;
+};
+
+
+template<class FirstBlock, class FirstX, class FirstY,
+ class SecondBlock, class SecondX, class SecondY,
+ class X, class Y>
+class BlockILU0Preconditioner : public Dune::Preconditioner<X, Y>
+{
+public:
+ //! \brief The matrix type the preconditioner is for.
+ // typedef typename std::remove_const<M>::type matrix_type;
+ //! \brief The domain type of the preconditioner.
+ typedef X domain_type;
+ //! \brief The range type of the preconditioner.
+ typedef Y range_type;
+ //! \brief The field type of the preconditioner.
+ typedef typename X::field_type field_type;
+
+ /*! \brief Constructor.
+
+ Constructor gets all parameters to operate the prec.
+ \param A The matrix to operate on.
+ \param w The relaxation factor.
+ */
+ BlockILU0Preconditioner (const FirstBlock& A, const SecondBlock& B)
+ : firstILU_(A, 1.0), secondILU_(B, 1.0) // does ILU decomposition
+ {}
+
+ /*!
+ \brief Prepare the preconditioner.
+
+ \copydoc Preconditioner::pre(X&,Y&)
+ */
+ virtual void pre (X& x, Y& b)
+ {
+ DUNE_UNUSED_PARAMETER(x);
+ DUNE_UNUSED_PARAMETER(b);
+ }
+
+ /*!
+ \brief Apply the preconditoner.
+
+ \copydoc Preconditioner::apply(X&,const Y&)
+ */
+ virtual void apply (X& v, const Y& d)
+ {
+ using namespace Dune::Indices;
+ firstILU_.apply(v[_0], d[_0]);
+ secondILU_.apply(v[_1], d[_1]);
+ }
+
+ /*!
+ \brief Clean up.
+
+ \copydoc Preconditioner::post(X&)
+ */
+ virtual void post (X& x)
+ {
+ DUNE_UNUSED_PARAMETER(x);
+ }
+
+ //! Category of the preconditioner (see SolverCategory::Category)
+ virtual Dune::SolverCategory::Category category() const
+ {
+ return Dune::SolverCategory::sequential;
+ }
+
+ private:
+ Dune::SeqILU0<FirstBlock, FirstX, FirstY> firstILU_;
+ Dune::SeqILU0<SecondBlock, SecondX, SecondY> secondILU_;
+};
+
+template <class TypeTag>
+class BlockILU0BiCGSTABSolver
+{
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+
+public:
+ BlockILU0BiCGSTABSolver() = default;
+ BlockILU0BiCGSTABSolver(const Problem& problem) {}
+
+ // expects a system as a multi-type matrix
+ // | A B |
+ // | C D |
+
+ // Solve saddle-point problem using a Schur complement based preconditioner
+ template<class Matrix, class Vector>
+ bool solve(const Matrix& M, Vector& x, const Vector& b)
+ {
+ int verbosity = GET_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, Verbosity);
+ static const double reduction = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, ResidualReduction);
+ static const double maxIter = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, MaxIterations);
+
+ using AType = typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockBulk;
+ using DType = typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockLowDim;
+
+ using BulkVector = typename GET_PROP(TypeTag, SolutionVector)::SolutionVectorBulk;
+ using EmbeddedVector = typename GET_PROP(TypeTag, SolutionVector)::SolutionVectorLowDim;
+
+ using namespace Dune::Indices;
+ const AType& A = M[_0][_0];
+ const DType& D = M[_1][_1];
+
+ // Solver to compute applications of D^-1 used in the preconditioner
+ // We construct it here so e.g. LU decomposition is only done once
+ using InnerSolver = Dune::UMFPack<DType>;
+ auto invD = std::make_shared<InnerSolver>(D, false);
+
+ BlockILU0Preconditioner<AType, BulkVector, BulkVector, DType, EmbeddedVector, EmbeddedVector, Vector, Vector> preconditioner(A, D);
+ Dune::MatrixAdapter<Matrix, Vector, Vector> op(M);
+ Dune::BiCGSTABSolver<Vector> solver(op, preconditioner, reduction, maxIter, verbosity);
+ auto bTmp(b);
+ solver.apply(x, bTmp, result_);
+
+ return result_.converged;
+ }
+
+ const Dune::InverseOperatorResult& result() const
+ {
+ return result_;
+ }
+
+private:
+ Dune::InverseOperatorResult result_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/mixeddimension/embedded/1p2c_richards2c/start.hh b/test/mixeddimension/embedded/1p2c_richards2c/start.hh
index 3bebc455759c21f0b55ce3ebff3a20ab87eb3b16..8ed49d6955440fef0b653cdd7d8d600e80e3e0b9 100644
--- a/test/mixeddimension/embedded/1p2c_richards2c/start.hh
+++ b/test/mixeddimension/embedded/1p2c_richards2c/start.hh
@@ -150,8 +150,8 @@ int start_(int argc,
TimeManager timeManager;
Problem problem(timeManager, bulkGridView, lowDimGridView);
- // locally refine levels deep around the embedded grid
- int levels = 2;
+ //locally refine levels deep around the embedded grid
+ int levels = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, SoilGrid, LocalRefinement);
for (int i = 0; i < levels; ++i)
{
CCMultiDimensionGlue<TypeTag> glue(problem.bulkProblem(), problem.lowDimProblem());
@@ -159,8 +159,49 @@ int start_(int argc,
// refine all 3D cells intersected
for (const auto& is : intersections(glue))
+ {
for (unsigned int outsideIdx = 0; outsideIdx < is.neighbor(0); ++outsideIdx)
- BulkGridCreator::grid().mark(1, is.outside(outsideIdx));
+ {
+ const auto cutElement = is.outside(outsideIdx);
+
+ // mark the cut element and all it's neighbors
+ BulkGridCreator::grid().mark(1, cutElement);
+ for (const auto& intersection : intersections(bulkGridView, cutElement))
+ if (intersection.neighbor())
+ BulkGridCreator::grid().mark(1, intersection.outside());
+ }
+
+ }
+
+ // refine all 3D cells that are where the contamination is
+ for (const auto& element : elements(bulkGridView))
+ {
+ const auto globalPos = element.geometry().center();
+ auto contaminationPos = problem.bulkProblem().bBoxMax()-problem.bulkProblem().bBoxMin();
+ contaminationPos[0] *= 0.25;
+ contaminationPos[1] *= 0.55;
+ contaminationPos[2] *= 0.25;
+ contaminationPos += problem.bulkProblem().bBoxMin();
+
+ static const Scalar extend = 0.15*(problem.bulkProblem().bBoxMax()[0]-problem.bulkProblem().bBoxMin()[0]);
+ if ((globalPos - contaminationPos).infinity_norm() < extend + 1e-7)
+ BulkGridCreator::grid().mark(1, element);
+ }
+
+ BulkGridCreator::grid().preAdapt();
+ BulkGridCreator::grid().adapt();
+ BulkGridCreator::grid().postAdapt();
+
+ // make sure there is only one level difference
+ for (const auto& element : elements(bulkGridView))
+ {
+ for (const auto& intersection : intersections(bulkGridView, element))
+ {
+ if (intersection.neighbor())
+ if (intersection.outside().level()-1 > element.level())
+ BulkGridCreator::grid().mark(1, element);
+ }
+ }
BulkGridCreator::grid().preAdapt();
BulkGridCreator::grid().adapt();
diff --git a/test/mixeddimension/embedded/1p2c_richards2c/test_rosi2c.input b/test/mixeddimension/embedded/1p2c_richards2c/test_rosi2c.input
index 07c683890fccd812d81cc9a086bdd8e4f599c539..a2016460e23ac89ef95148009a1a8e87cdde9618 100644
--- a/test/mixeddimension/embedded/1p2c_richards2c/test_rosi2c.input
+++ b/test/mixeddimension/embedded/1p2c_richards2c/test_rosi2c.input
@@ -8,21 +8,23 @@ EpisodeTime = 8640 # [s]
[Grid]
File = ./grids/lupine.dgf
-Refinement = 0
+Refinement = 1
[SoilGrid]
LowerLeft = -0.05 -0.05 -0.1
UpperRight = 0.05 0.05 0
-Cells = 15 15 15
+Cells = 10 10 10
ClosureType = None
+Refinement = 0
+LocalRefinement = 0
[Problem]
Name = rosi2c
-ContaminantMoleFraction = 1e-6
+ContaminantMoleFraction = 3e-7
MolarMass = 20e-3 # in kg/mol, molar mass heavy water D2O
[Component]
-LiquidDiffusionCoefficient = 1e-9
+LiquidDiffusionCoefficient = 2.3e-9
Name = D2O
[SpatialParams]