[stokesdarcy] Add first draft for iterative coupling

appl/coupling-ff-pm/CMakeLists.txt

@@ -2,4 +2,4 @@
 add_custom_target(test_exercises)
 
 add_subdirectory(monolithic)
-add_subdirectory(solution)
+add_subdirectory(iterative)

appl/coupling-ff-pm/iterative/1pspatialparams.hh

+// -*- 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
+ * \ingroup OnePTests
+ * \brief The spatial parameters class for the test problem using the 1p cc model
+ */
+#ifndef DUMUX_1P_TEST_SPATIALPARAMS_HH
+#define DUMUX_1P_TEST_SPATIALPARAMS_HH
+
+#include <dumux/material/spatialparams/fv1p.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup OnePModel
+ *
+ * \brief The spatial parameters class for the test problem using the
+ *        1p cc model
+ */
+template<class FVGridGeometry, class Scalar>
+class OnePSpatialParams
+: public FVSpatialParamsOneP<FVGridGeometry, Scalar, OnePSpatialParams<FVGridGeometry, Scalar>>
+{
+    using GridView = typename FVGridGeometry::GridView;
+    using ParentType = FVSpatialParamsOneP<FVGridGeometry, Scalar, OnePSpatialParams<FVGridGeometry, Scalar>>;
+
+    using Element = typename GridView::template Codim<0>::Entity;
+    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+    // export permeability type
+    using PermeabilityType = Scalar;
+
+    OnePSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    : ParentType(fvGridGeometry)
+    {
+        permeability_ = getParam<Scalar>("Darcy.SpatialParams.Permeability");
+        porosity_ = getParam<Scalar>("Darcy.SpatialParams.Porosity");
+        alphaBJ_ = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph");
+    }
+
+    /*!
+     * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
+     *
+     * \param globalPos The global position
+     * \return the intrinsic permeability
+     */
+    PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
+    { return permeability_; }
+
+    /*! \brief Define the porosity in [-].
+     *
+     * \param globalPos The global position
+     */
+    Scalar porosityAtPos(const GlobalPosition& globalPos) const
+    { return porosity_; }
+
+    /*! \brief Define the Beavers-Joseph coefficient in [-].
+     *
+     * \param globalPos The global position
+     */
+    Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const
+    { return alphaBJ_; }
+
+
+private:
+    Scalar permeability_;
+    Scalar porosity_;
+    Scalar alphaBJ_;
+};
+
+} // end namespace
+
+#endif

appl/coupling-ff-pm/iterative/CMakeLists.txt

+# executables for the iterative case
+dune_add_test(NAME test_ff
+              SOURCES main_ff.cc
+              COMPILE_DEFINITIONS ENABLEMONOLITHIC=0)
+
+dune_add_test(NAME test_pm
+              SOURCES main_pm.cc
+              COMPILE_DEFINITIONS ENABLEMONOLITHIC=0)
+
+# add a symlink for each input file
+add_input_file_links()

appl/coupling-ff-pm/iterative/main_ff.cc

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief A test problem for the coupled Stokes/Darcy problem (1p)
+ */
+ #include <config.h>
+
+ #include <ctime>
+ #include <iostream>
+
+ #include <dune/common/parallel/mpihelper.hh>
+ #include <dune/common/timer.hh>
+ #include <dune/istl/io.hh>
+
+ #include <dumux/common/properties.hh>
+ #include <dumux/common/parameters.hh>
+ #include <dumux/common/partial.hh>
+ #include <dumux/common/dumuxmessage.hh>
+ #include <dumux/linear/seqsolverbackend.hh>
+ #include <dumux/assembly/fvassembler.hh>
+ #include <dumux/assembly/diffmethod.hh>
+ #include <dumux/discretization/method.hh>
+ #include <dumux/io/vtkoutputmodule.hh>
+ #include <dumux/io/staggeredvtkoutputmodule.hh>
+ #include <dumux/io/grid/gridmanager.hh>
+
+ #include <dumux/assembly/staggeredfvassembler.hh>
+ #include <dumux/nonlinear/newtonsolver.hh>
+
+#include "../monolithic/ffproblem.hh"
+
+int main(int argc, char** argv) try
+{
+    using namespace Dumux;
+
+    // initialize MPI, finalize is done automatically on exit
+    const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+    // print dumux start message
+    if (mpiHelper.rank() == 0)
+        DumuxMessage::print(/*firstCall=*/true);
+
+    // parse command line arguments and input file
+    Parameters::init(argc, argv);
+
+    // Define the sub problem type tags
+    using FreeFlowTypeTag = Properties::TTag::FreeFlowModel;
+
+    // try to create a grid (from the given grid file or the input file)
+    using FreeFlowGridManager = Dumux::GridManager<GetPropType<FreeFlowTypeTag, Properties::Grid>>;
+    FreeFlowGridManager freeFlowGridManager;
+    freeFlowGridManager.init("FreeFlow"); // pass parameter group
+
+    // we compute on the leaf grid view
+    const auto& freeFlowGridView = freeFlowGridManager.grid().leafGridView();
+
+    // create the finite volume grid geometry
+    using FreeFlowFVGridGeometry = GetPropType<FreeFlowTypeTag, Properties::FVGridGeometry>;
+    auto freeFlowFvGridGeometry = std::make_shared<FreeFlowFVGridGeometry>(freeFlowGridView);
+    freeFlowFvGridGeometry->update();
+
+    // the problem (initial and boundary conditions)
+    using FreeFlowProblem = GetPropType<FreeFlowTypeTag, Properties::Problem>;
+    auto freeFlowProblem = std::make_shared<FreeFlowProblem>(freeFlowFvGridGeometry);
+
+    // the solution vector
+    GetPropType<FreeFlowTypeTag, Properties::SolutionVector> sol;
+    sol[FreeFlowFVGridGeometry::cellCenterIdx()].resize(freeFlowFvGridGeometry->numCellCenterDofs());
+    sol[FreeFlowFVGridGeometry::faceIdx()].resize(freeFlowFvGridGeometry->numFaceDofs());
+
+
+    // apply initial solution for instationary problems
+    freeFlowProblem->applyInitialSolution(sol);
+
+    // the grid variables
+    using FreeFlowGridVariables = GetPropType<FreeFlowTypeTag, Properties::GridVariables>;
+    auto freeFlowGridVariables = std::make_shared<FreeFlowGridVariables>(freeFlowProblem, freeFlowFvGridGeometry);
+    freeFlowGridVariables->init(sol);
+
+    // intialize the vtk output module
+    StaggeredVtkOutputModule<FreeFlowGridVariables, decltype(sol)> freeFlowVtkWriter(*freeFlowGridVariables, sol, freeFlowProblem->name());
+    GetPropType<FreeFlowTypeTag, Properties::IOFields>::initOutputModule(freeFlowVtkWriter);
+    freeFlowVtkWriter.addField(freeFlowProblem->getAnalyticalVelocityX(), "analyticalV_x");
+    freeFlowVtkWriter.write(0.0);
+
+    // the assembler for a stationary problem
+    using Assembler = StaggeredFVAssembler<FreeFlowTypeTag, DiffMethod::numeric>;
+    auto assembler = std::make_shared<Assembler>(freeFlowProblem, freeFlowFvGridGeometry, freeFlowGridVariables);
+
+    // the linear solver
+    using LinearSolver = UMFPackBackend;
+    auto linearSolver = std::make_shared<LinearSolver>();
+
+    // the non-linear solver
+    using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
+    NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+    // solve the non-linear system
+    nonLinearSolver.solve(sol);
+
+    // write vtk output
+    freeFlowVtkWriter.write(1.0);
+
+    ////////////////////////////////////////////////////////////
+    // finalize, print dumux message to say goodbye
+    ////////////////////////////////////////////////////////////
+
+    // print dumux end message
+    if (mpiHelper.rank() == 0)
+    {
+        Parameters::print();
+        DumuxMessage::print(/*firstCall=*/false);
+    }
+
+    return 0;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+    std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+    return 1;
+}
+catch (Dune::DGFException & e)
+{
+    std::cerr << "DGF exception thrown (" << e <<
+                 "). Most likely, the DGF file name is wrong "
+                 "or the DGF file is corrupted, "
+                 "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+                 << " ---> Abort!" << std::endl;
+    return 2;
+}
+catch (Dune::Exception &e)
+{
+    std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+    return 3;
+}
+catch (...)
+{
+    std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+    return 4;
+}

appl/coupling-ff-pm/iterative/main_pm.cc

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief A test problem for the coupled Stokes/Darcy problem (1p)
+ */
+ #include <config.h>
+
+ #include <ctime>
+ #include <iostream>
+
+ #include <dune/common/parallel/mpihelper.hh>
+ #include <dune/common/timer.hh>
+ #include <dune/grid/io/file/dgfparser/dgfexception.hh>
+ #include <dune/grid/io/file/vtk.hh>
+ #include <dune/istl/io.hh>
+
+ #include <dumux/common/properties.hh>
+ #include <dumux/common/parameters.hh>
+ #include <dumux/common/valgrind.hh>
+ #include <dumux/common/dumuxmessage.hh>
+ #include <dumux/common/defaultusagemessage.hh>
+
+ #include <dumux/linear/amgbackend.hh>
+ #include <dumux/nonlinear/newtonsolver.hh>
+
+ #include <dumux/assembly/fvassembler.hh>
+ #include <dumux/assembly/diffmethod.hh>
+
+ #include <dumux/discretization/method.hh>
+
+ #include <dumux/io/vtkoutputmodule.hh>
+ #include <dumux/io/grid/gridmanager.hh>
+
+ #include "../monolithic/pmproblem.hh"
+
+ /*!
+  * \brief Returns the pressure at the interface using Darcy's law for reconstruction
+  */
+ template<class Problem, class Element, class SubControlVolumeFace, class VolumeVariables, class Scalar>
+ Scalar pressureAtInterface(const Problem& problem,
+                            const Element& element,
+                            const SubControlVolumeFace& scvf,
+                            const VolumeVariables& volVars,
+                            const Scalar boundaryFlux)
+ {
+     const Scalar cellCenterPressure = volVars.pressure();
+
+     // v = -K/mu * (gradP + rho*g)
+     auto velocity = scvf.unitOuterNormal();
+     velocity *= boundaryFlux; // TODO check sign
+
+    const Scalar ccPressure = volVars.pressure();
+    const Scalar mobility = volVars.mobility();
+    const Scalar density = volVars.density();
+    const auto K = volVars.permeability();
+
+    // v = -kr/mu*K * (gradP + rho*g) = -mobility*K * (gradP + rho*g)
+
+    const auto alpha = vtmv(scvf.unitOuterNormal(), K, problem().spatialParams.gravity(scvf.center()));
+
+    auto distanceVector = scvf.center() - element.geometry().center();
+    distanceVector /= distanceVector.two_norm2();
+    const Scalar ti = vtmv(distanceVector, K, scvf.unitOuterNormal());
+
+    return (1/mobility * (scvf.unitOuterNormal() * velocity) + density * alpha)/ti
+           + ccPressure;
+ }
+
+int main(int argc, char** argv) try
+{
+    using namespace Dumux;
+
+    // initialize MPI, finalize is done automatically on exit
+    const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+    // print dumux start message
+    if (mpiHelper.rank() == 0)
+        DumuxMessage::print(/*firstCall=*/true);
+
+    // parse command line arguments and input file
+    Parameters::init(argc, argv);
+
+    using DarcyTypeTag = Properties::TTag::DarcyOneP;
+
+    using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
+    DarcyGridManager darcyGridManager;
+    darcyGridManager.init("Darcy"); // pass parameter group
+
+    // we compute on the leaf grid view
+    const auto& darcyGridView = darcyGridManager.grid().leafGridView();
+
+    // create the finite volume grid geometry
+    using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::FVGridGeometry>;
+    auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
+    darcyFvGridGeometry->update();
+
+    using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
+    auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry);
+
+    // the solution vector
+    GetPropType<DarcyTypeTag, Properties::SolutionVector> sol;
+    sol.resize(darcyFvGridGeometry->numDofs());
+
+    darcyProblem->applyInitialSolution(sol);
+
+    // the grid variables
+    using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
+    auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
+    darcyGridVariables->init(sol);
+
+    // intialize the vtk output module
+    const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
+
+    VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol, darcyName);
+    using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
+    darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
+    GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
+    darcyVtkWriter.write(0.0);
+
+    // the assembler for a stationary problem
+    using Assembler = FVAssembler<DarcyTypeTag, DiffMethod::numeric>;
+    auto assembler = std::make_shared<Assembler>(darcyProblem, darcyFvGridGeometry, darcyGridVariables);
+
+    // the linear solver
+    using LinearSolver = UMFPackBackend;
+    auto linearSolver = std::make_shared<LinearSolver>();
+
+    // the non-linear solver
+    using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
+    NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+    // solve the non-linear system
+    nonLinearSolver.solve(sol);
+
+    // write vtk output
+    darcyVtkWriter.write(1.0);
+
+    ////////////////////////////////////////////////////////////
+    // finalize, print dumux message to say goodbye
+    ////////////////////////////////////////////////////////////
+
+    // print dumux end message
+    if (mpiHelper.rank() == 0)
+    {
+        Parameters::print();
+        DumuxMessage::print(/*firstCall=*/false);
+    }
+
+    return 0;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+    std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+    return 1;
+}
+catch (Dune::DGFException & e)
+{
+    std::cerr << "DGF exception thrown (" << e <<
+                 "). Most likely, the DGF file name is wrong "
+                 "or the DGF file is corrupted, "
+                 "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+                 << " ---> Abort!" << std::endl;
+    return 2;
+}
+catch (Dune::Exception &e)
+{
+    std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+    return 3;
+}
+catch (...)
+{
+    std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+    return 4;
+}

appl/coupling-ff-pm/iterative/params.input

+ # for dune-subgrid
+ [Grid]
+Positions0 = 0 1
+Positions1 = 0 0.2 0.3 0.65
+Cells0 = 100
+Cells1 = 10 50 18
+Baseline = 0.25 # [m]
+Amplitude = 0.04 # [m]
+Offset = 0.5 # [m]
+Scaling = 0.2 #[m]
+
+[Stokes.Grid]
+Verbosity = true
+Positions0 = 0.0 1.0
+Positions1 = 1.0 2.0
+Cells0 = 20
+Cells1 = 100
+Grading1 = 1
+
+[Darcy.Grid]
+Verbosity = true
+Positions0 = 0.0 1.0
+Positions1 = 0.0 1.0
+Cells0 = 20
+Cells1 = 20
+Grading1 = 1
+
+[Stokes.Problem]
+Name = stokes
+PressureDifference = 1e-9
+EnableInertiaTerms = false
+
+[Darcy.Problem]
+Name = darcy
+
+[Darcy.SpatialParams]
+Permeability = 1e-6 # m^2
+Porosity = 0.4
+AlphaBeaversJoseph = 1.0
+
+[Problem]
+Name = ex_ff-pm-interface
+EnableGravity = false
+
+[Vtk]
+AddVelocity = 1

appl/coupling-ff-pm/monolithic/CMakeLists.txt

-add_subdirectory(interface)
-add_subdirectory(models)
-add_subdirectory(turbulence)
+# executables for the monolithic case
+dune_add_test(NAME ff-pm_monolithic
+              SOURCES main.cc)
+
+# add a symlink for each input file
+add_input_file_links()

appl/coupling-ff-pm/monolithic/ffproblem.hh

@@ -23,10 +23,11 @@
 #ifndef DUMUX_STOKES_SUBPROBLEM_HH
 #define DUMUX_STOKES_SUBPROBLEM_HH
 
-#include <dune/grid/yaspgrid.hh>
+#ifndef ENABLEMONOLITHIC
+#define ENABLEMONOLITHIC 1
+#endif
 
-//****** uncomment for the last exercise *****//
-// #include <dumux/io/grid/subgridgridcreator.hh>
+#include <dune/grid/yaspgrid.hh>
 
 #include <dumux/material/fluidsystems/1pliquid.hh>
 #include <dumux/material/components/simpleh2o.hh>
@@ -44,12 +45,12 @@ namespace Properties
 {
 // Create new type tags
 namespace TTag {
-struct StokesOneP { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
+struct FreeFlowModel { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
 } // end namespace TTag
 
 // the fluid system
 template<class TypeTag>
-struct FluidSystem<TypeTag, TTag::StokesOneP>
+struct FluidSystem<TypeTag, TTag::FreeFlowModel>
 {
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
@@ -57,7 +58,7 @@ struct FluidSystem<TypeTag, TTag::StokesOneP>
 
 // Set the grid type
 template<class TypeTag>
-struct Grid<TypeTag, TTag::StokesOneP>
+struct Grid<TypeTag, TTag::FreeFlowModel>
 {
     static constexpr auto dim = 2;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
@@ -73,14 +74,14 @@ struct Grid<TypeTag, TTag::StokesOneP>
 
 // Set the problem property
 template<class TypeTag>
-struct Problem<TypeTag, TTag::StokesOneP> { using type = Dumux::StokesSubProblem<TypeTag> ; };
+struct Problem<TypeTag, TTag::FreeFlowModel> { using type = Dumux::StokesSubProblem<TypeTag> ; };
 
 template<class TypeTag>
-struct EnableFVGridGeometryCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; };
+struct EnableFVGridGeometryCache<TypeTag, TTag::FreeFlowModel> { static constexpr bool value = true; };
 template<class TypeTag>
-struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; };
+struct EnableGridFluxVariablesCache<TypeTag, TTag::FreeFlowModel> { static constexpr bool value = true; };
 template<class TypeTag>
-struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; };
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::FreeFlowModel> { static constexpr bool value = true; };
 }
 
 /*!
@@ -109,11 +110,18 @@ class StokesSubProblem : public NavierStokesProblem<TypeTag>
     using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
 
+#if ENABLEMONOLITHIC
     using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+#endif
 
 public:
+#if ENABLEMONOLITHIC
     StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
     : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+#else
+    StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6)
+#endif
     {
         deltaP_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference");
     }
@@ -162,8 +170,6 @@ public:
         // left/right wall
         if (onRightBoundary_(globalPos) || (onLeftBoundary_(globalPos)))
         {
-//            values.setDirichlet(Indices::velocityXIdx);
-//            values.setDirichlet(Indices::velocityYIdx);
           values.setDirichlet(Indices::pressureIdx);
         }
         else
@@ -173,13 +179,22 @@ public:
         }
 
         // coupling interface
+#if ENABLEMONOLITHIC
         if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
         {
             values.setCouplingNeumann(Indices::conti0EqIdx);
             values.setCouplingNeumann(Indices::momentumYBalanceIdx);
-//            values.setDirichlet(Indices::velocityXIdx); // assume no slip on interface
-            values.setBJS(Indices::momentumXBalanceIdx); // assume no slip on interface
+            values.setBJS(Indices::momentumXBalanceIdx);
         }
+#else
+    // // TODO do preCICE stuff in analogy to heat transfer
+    // if(/*preCice*/)
+    // {
+    //     values.setCouplingNeumann(Indices::conti0EqIdx);
+    //     values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+    //     values.setBJS(Indices::momentumXBalanceIdx);
+    // }
+#endif
 
         return values;
     }
@@ -193,18 +208,6 @@ public:
     {
         PrimaryVariables values(0.0);
         values = initialAtPos(globalPos);
-
-//        if(onUpperBoundary_(globalPos))
-//        {
-//          values[Indices::velocityXIdx] = 0.0;
-//          values[Indices::velocityYIdx] = 0.0;
-//        }
-
-//        if(onLeftBoundary_(globalPos))
-//            values[Indices::pressureIdx] = deltaP_;
-//        if(onRightBoundary_(globalPos))
-//            values[Indices::pressureIdx] = 0.0;
-
         return values;
     }
 
@@ -226,24 +229,31 @@ public:
     {
         NumEqVector values(0.0);
 
+#if ENABLEMONOLITHIC
         if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
         {
             values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
             values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
         }
+#else
+        // if(/*preCICE*/)
+        // {
+        //     const Scalar density = 1000; // TODO how to handle compressible fluids?
+        //     values[Indices::conti0EqIdx] = elemFaceVars[scvf].velocitySelf() * scvf.directionSign() * density;
+        //     values[Indices::momentumYBalanceIdx] = /*pressure from Darcy*/
+        // }
+#endif
 
         return values;
     }
 
     // \}
 
-    //! Set the coupling manager
-    void setCouplingManager(std::shared_ptr<CouplingManager> cm)
-    { couplingManager_ = cm; }
-
+#if ENABLEMONOLITHIC
     //! Get the coupling manager
     const CouplingManager& couplingManager() const
     { return *couplingManager_; }
+#endif
 
    /*!
      * \name Volume terms
@@ -272,7 +282,11 @@ public:
      */
     Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const
     {
+#if ENABLEMONOLITHIC
         return couplingManager().couplingData().darcyPermeability(element, scvf);
+#else
+        return 1e-10; // TODO transfer information or just use constant value
+#endif
     }
 
     /*!
@@ -280,7 +294,11 @@ public:
      */
     Scalar alphaBJ(const SubControlVolumeFace& scvf) const
     {
+#if ENABLEMONOLITHIC
         return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
+#else
+        return 1.0; // TODO transfer information or just use constant value
+#endif
     }
 
     /*!
@@ -338,7 +356,9 @@ private:
     Scalar eps_;
     Scalar deltaP_;
 
+#if ENABLEMONOLITHIC
     std::shared_ptr<CouplingManager> couplingManager_;
+#endif
 
     mutable std::vector<Scalar> analyticalVelocityX_;
 };

appl/coupling-ff-pm/monolithic/main.cc

@@ -46,14 +46,14 @@
 
 #include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
 
-#include "ex_interface_pmproblem.hh"
-#include "ex_interface_ffproblem.hh"
+#include "pmproblem.hh"
+#include "ffproblem.hh"
 
 namespace Dumux {
 namespace Properties {
 
 template<class TypeTag>
-struct CouplingManager<TypeTag, TTag::StokesOneP>
+struct CouplingManager<TypeTag, TTag::FreeFlowModel>
 {
     using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOneP>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
@@ -62,7 +62,7 @@ struct CouplingManager<TypeTag, TTag::StokesOneP>
 template<class TypeTag>
 struct CouplingManager<TypeTag, TTag::DarcyOneP>
 {
-    using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesOneP, Properties::TTag::StokesOneP, TypeTag>;
+    using Traits = StaggeredMultiDomainTraits<Properties::TTag::FreeFlowModel, Properties::TTag::FreeFlowModel, TypeTag>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 
@@ -84,7 +84,7 @@ int main(int argc, char** argv) try
     Parameters::init(argc, argv);
 
     // Define the sub problem type tags
-    using StokesTypeTag = Properties::TTag::StokesOneP;
+    using StokesTypeTag = Properties::TTag::FreeFlowModel;
     using DarcyTypeTag = Properties::TTag::DarcyOneP;
 
 

appl/coupling-ff-pm/monolithic/pmproblem.hh

@@ -24,6 +24,10 @@
 #ifndef DUMUX_DARCY_SUBPROBLEM_HH
 #define DUMUX_DARCY_SUBPROBLEM_HH
 
+#ifndef ENABLEMONOLITHIC
+#define ENABLEMONOLITHIC 1
+#endif
+
 #include <dune/grid/yaspgrid.hh>
 
 //****** uncomment for the last exercise *****//
@@ -34,7 +38,7 @@
 #include <dumux/porousmediumflow/1p/model.hh>
 #include <dumux/porousmediumflow/problem.hh>
 
-#include "../1pspatialparams.hh"
+#include "1pspatialparams.hh"
 
 #include <dumux/material/components/simpleh2o.hh>
 #include <dumux/material/fluidsystems/1pliquid.hh>
@@ -109,12 +113,19 @@ class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
     using Element = typename GridView::template Codim<0>::Entity;
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
 
+#if ENABLEMONOLITHIC
     using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+#endif
 
 public:
+#if ENABLEMONOLITHIC
     DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
                    std::shared_ptr<CouplingManager> couplingManager)
     : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+#else
+DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7)
+#endif
     {}
 
     /*!
@@ -149,14 +160,11 @@ public:
         // set Neumann BCs to all boundaries first
         values.setAllNeumann();
 
+#if ENABLEMONOLITHIC
         // set the coupling boundary condition at the interface
         if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
             values.setAllCouplingNeumann();
-
-        // set a Dirichlet boundary condition at the bottom
-        //if (onLowerBoundary_(scvf.center()))
-        //    values.setAllDirichlet();
-
+#endif
         return values;
     }
 
@@ -196,10 +204,16 @@ public:
         // no-flow everywhere ...
         NumEqVector values(0.0);
 
+#if ENABLEMONOLITHIC
         // ... except at the coupling interface
         if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
             values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf);
-
+#else
+        // if (/*preCICE*/)
+        // {
+        //     values[Indices::conti0EqIdx] = /*mass flux from stokes*/;
+        // }
+#endif
         return values;
     }
 
@@ -242,13 +256,11 @@ public:
 
     // \}
 
-    //! Set the coupling manager
-    void setCouplingManager(std::shared_ptr<CouplingManager> cm)
-    { couplingManager_ = cm; }
-
+#if ENABLEMONOLITHIC
     //! Get the coupling manager
     const CouplingManager& couplingManager() const
     { return *couplingManager_; }
+#endif
 
 private:
     bool onLeftBoundary_(const GlobalPosition &globalPos) const
@@ -264,7 +276,10 @@ private:
     { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
 
     Scalar eps_;
+
+#if ENABLEMONOLITHIC
     std::shared_ptr<CouplingManager> couplingManager_;
+#endif
 };
 } //end namespace