diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux/1pspatialparams.hh b/appl/coupling-ff-pm/iterative-reversed-mono-flux/1pspatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a042388f1eac51bb38ef980cb0f183756aaf6362
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux/1pspatialparams.hh
@@ -0,0 +1,91 @@
+// -*- 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
diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux/ffproblem-reversed.hh b/appl/coupling-ff-pm/iterative-reversed-mono-flux/ffproblem-reversed.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0a9ef2af63083faefd885bd4ebe6fee7ab0b73c9
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux/ffproblem-reversed.hh
@@ -0,0 +1,414 @@
+// -*- 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 The free flow sub problem
+ */
+#ifndef DUMUX_STOKES_SUBPROBLEM_HH
+#define DUMUX_STOKES_SUBPROBLEM_HH
+
+#ifndef ENABLEMONOLITHIC
+#define ENABLEMONOLITHIC 0
+#endif
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/material/fluidsystems/1pliquid.hh>
+#include <dumux/material/components/simpleh2o.hh>
+
+#include <dumux/freeflow/navierstokes/problem.hh>
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+#include <dumux/freeflow/navierstokes/model.hh>
+
+#include "../../precice-adapter/include/preciceadapter.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class StokesSubProblem;
+
+namespace Properties
+{
+// Create new type tags
+namespace TTag {
+struct FreeFlowModel { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
+} // end namespace TTag
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::FreeFlowModel>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
+};
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::FreeFlowModel>
+{
+ static constexpr auto dim = 2;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using TensorGrid = Dune::YaspGrid<2, Dune::TensorProductCoordinates<Scalar, dim> >;
+
+//****** comment out for the last exercise *****//
+ using type = TensorGrid;
+
+//****** uncomment for the last exercise *****//
+ // using HostGrid = TensorGrid;
+ // using type = Dune::SubGrid<dim, HostGrid>;
+};
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::FreeFlowModel> { using type = Dumux::StokesSubProblem<TypeTag> ; };
+
+template<class TypeTag>
+struct EnableFVGridGeometryCache<TypeTag, TTag::FreeFlowModel> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::FreeFlowModel> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::FreeFlowModel> { static constexpr bool value = true; };
+}
+
+/*!
+ * \brief The free flow sub problem
+ */
+template <class TypeTag>
+class StokesSubProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ 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, "FreeFlow"),
+ eps_(1e-6),
+ couplingInterface_(precice_adapter::PreciceAdapter::getInstance() ),
+ pressureId_(0),
+ velocityId_(0),
+ dataIdsWereSet_(false)
+#endif
+ {
+ deltaP_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference");
+// pressureId_ = couplingInterface_.getIdFromName( "Pressure" );
+// velocityId_ = couplingInterface_.getIdFromName( "Velocity" );
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature() const
+ { return 273.15 + 10; } // 10°C
+
+ /*!
+ * \brief Return the sources within the domain.
+ *
+ * \param globalPos The global position
+ */
+ NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+ { return NumEqVector(0.0); }
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ *
+ * \param element The finite element
+ * \param scvf The sub control volume face
+ */
+ BoundaryTypes boundaryTypes(const Element& element,
+ const SubControlVolumeFace& scvf) const
+ {
+ BoundaryTypes values;
+
+ const auto& globalPos = scvf.dofPosition();
+ const auto faceId = scvf.index();
+
+ // left/right wall
+ if (onRightBoundary_(globalPos) || (onLeftBoundary_(globalPos)))
+ {
+ values.setDirichlet(Indices::pressureIdx);
+ }
+
+
+ // coupling interface
+#if ENABLEMONOLITHIC
+ else if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setBJS(Indices::momentumXBalanceIdx);
+ }
+#else
+
+ else if ( couplingInterface_.isCoupledEntity(faceId) )
+ {
+ // // TODO do preCICE stuff in analogy to heat transfer
+ assert( dataIdsWereSet_ );
+ //TODO What do I want to do here?
+ // values.setCouplingNeumann(Indices::conti0EqIdx);
+ // values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+
+// values.setNeumann(Indices::conti0EqIdx);
+// values.setNeumann(Indices::momentumYBalanceIdx);
+ values.setBJS(Indices::momentumXBalanceIdx);
+ }
+#endif
+ else
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param globalPos The global position
+ */
+ using ParentType::dirichlet;
+ PrimaryVariables dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
+ {
+ PrimaryVariables values(0.0);
+ values = initialAtPos(scvf.center());
+
+ const auto faceId = scvf.index();
+ if( couplingInterface_.isCoupledEntity( faceId ) )
+ {
+ values[Indices::velocityYIdx] =
+ couplingInterface_.getScalarQuantityOnFace( velocityId_, faceId );
+ }
+
+
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param elemFaceVars The element face variables
+ * \param scvf The boundary sub control volume face
+ */
+ template<class ElementVolumeVariables, class ElementFaceVariables>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ 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
+ assert( dataIdsWereSet_ );
+ const auto faceId = scvf.index();
+ if( couplingInterface_.isCoupledEntity( faceId ) )
+ {
+ const Scalar density = 1000; // TODO how to handle compressible fluids?
+ const auto& volVars = elemVolVars[scvf.insideScvIdx()];
+ const Scalar density_ = volVars.density();
+ values[Indices::conti0EqIdx] = density * elemFaceVars[scvf].velocitySelf() * scvf.directionSign();
+ values[Indices::momentumYBalanceIdx] = scvf.directionSign() * (couplingInterface_.getScalarQuantityOnFace( pressureId_, faceId ) - initialAtPos(scvf.center())[Indices::pressureIdx]) ;
+ }
+#endif
+
+ return values;
+ }
+
+ // \}
+
+#if ENABLEMONOLITHIC
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+#endif
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+ //values[Indices::velocityYIdx] = -1e-6 * globalPos[0] * (this->fvGridGeometry().bBoxMax()[0] - globalPos[0]);
+ if(onLeftBoundary_(globalPos))
+ values[Indices::pressureIdx] = deltaP_;
+ if(onRightBoundary_(globalPos))
+ values[Indices::pressureIdx] = 0.0;
+
+ return values;
+ }
+
+ /*!
+ * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ 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
+ }
+
+ /*!
+ * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ 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
+ }
+
+ /*!
+ * \brief calculate the analytical velocity in x direction based on Beavers & Joseph (1967)
+ */
+ void calculateAnalyticalVelocityX() const
+ {
+ analyticalVelocityX_.resize(this->fvGridGeometry().gridView().size(0));
+
+ using std::sqrt;
+ const Scalar dPdX = -deltaP_ / (this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0]);
+ static const Scalar mu = FluidSystem::viscosity(temperature(), 1e5);
+ static const Scalar alpha = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph");
+ static const Scalar K = getParam<Scalar>("Darcy.SpatialParams.Permeability");
+ static const Scalar sqrtK = sqrt(K);
+ const Scalar sigma = (this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1])/sqrtK;
+
+ const Scalar uB = -K/(2.0*mu) * ((sigma*sigma + 2.0*alpha*sigma) / (1.0 + alpha*sigma)) * dPdX;
+
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ const auto eIdx = this->fvGridGeometry().gridView().indexSet().index(element);
+ const Scalar y = element.geometry().center()[1] - this->fvGridGeometry().bBoxMin()[1];
+
+ const Scalar u = uB*(1.0 + alpha/sqrtK*y) + 1.0/(2.0*mu) * (y*y + 2*alpha*y*sqrtK) * dPdX;
+ analyticalVelocityX_[eIdx] = u;
+ }
+ }
+
+ /*!
+ * \brief Get the analytical velocity in x direction
+ */
+ const std::vector<Scalar>& getAnalyticalVelocityX() const
+ {
+ if(analyticalVelocityX_.empty())
+ calculateAnalyticalVelocityX();
+ return analyticalVelocityX_;
+ }
+
+#if !ENABLEMONOLITHIC
+ void updatePreciceDataIds()
+ {
+ pressureId_ = couplingInterface_.getIdFromName( "Pressure" );
+ velocityId_ = couplingInterface_.getIdFromName( "Velocity" );
+ dataIdsWereSet_ = true;
+ }
+#endif
+
+ // \}
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ Scalar eps_;
+ Scalar deltaP_;
+
+#if ENABLEMONOLITHIC
+ std::shared_ptr<CouplingManager> couplingManager_;
+#else
+ precice_adapter::PreciceAdapter& couplingInterface_;
+ size_t pressureId_;
+ size_t velocityId_;
+ bool dataIdsWereSet_;
+#endif
+
+ mutable std::vector<Scalar> analyticalVelocityX_;
+};
+} //end namespace
+
+#endif // DUMUX_STOKES_SUBPROBLEM_HH
diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux/main_ff-reversed.cc b/appl/coupling-ff-pm/iterative-reversed-mono-flux/main_ff-reversed.cc
new file mode 100644
index 0000000000000000000000000000000000000000..3ffa6979e73354a367ce5a98526a690bd7b08a4e
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux/main_ff-reversed.cc
@@ -0,0 +1,635 @@
+// -*- 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 "ffproblem-reversed.hh"
+
+#include "../../precice-adapter/include/preciceadapter.hh"
+
+#include "monolithicdata.hh"
+
+//TODO
+// Helper function to put pressure on interface
+
+//namespace MonolithicSolution {
+// const std::array<double, 20> velocity = {
+// -2.14294737477106e-12,
+// -1.414265166923949e-12,
+// -1.04852776201853e-12,
+// -8.135033312884269e-13,
+// -6.37657086620555e-13,
+// -4.93766950668703e-13,
+// -3.690552916814545e-13,
+// -2.563759402825702e-13,
+// -1.511247369413653e-13,
+// -4.994226002810334e-14,
+// 4.994226002813352e-14,
+// 1.511247369413988e-13,
+// 2.563759402826087e-13,
+// 3.690552916815079e-13,
+// 4.937669506687764e-13,
+// 6.376570866206665e-13,
+// 8.13503331288573e-13,
+// 1.048527762018745e-12,
+// 1.414265166924232e-12,
+// 2.142947374771422e-12 };
+
+// const std::array<double, 20> pressure = {
+// 9.749674884723169e-10,
+// 9.249954142532818e-10,
+// 8.749977731157065e-10,
+// 8.24998517020602e-10,
+// 7.749988528034125e-10,
+// 7.249990852158898e-10,
+// 6.749992924234883e-10,
+// 6.249994944809285e-10,
+// 5.749996962877403e-10,
+// 5.249998986713713e-10,
+// 4.750001013286287e-10,
+// 4.250003037122596e-10,
+// 3.750005055190715e-10,
+// 3.250007075765118e-10,
+// 2.750009147841103e-10,
+// 2.250011471965874e-10,
+// 1.750014829793979e-10,
+// 1.250022268842933e-10,
+// 7.500458574671815e-11,
+// 2.503251152768306e-11 };
+//}
+
+template<class ElementFaceVariables, class SubControlVolumeFace>
+auto velocityAtInterface(const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf)
+{
+ const double scalarVelocity = elemFaceVars[scvf].velocitySelf();
+ auto velocity = scvf.unitOuterNormal();
+ velocity[scvf.directionIndex()] = scalarVelocity;
+ return velocity;
+ }
+
+template<class FluxVariables,
+ class Problem,
+ class Element,
+ class SubControlVolumeFace,
+ class FVElementGeometry,
+ class ElementVolumeVariables,
+ class ElementFaceVariables,
+ class ElementFluxVariablesCache>
+auto pressureAtInterface(const Problem& problem,
+ const Element& element,
+ const SubControlVolumeFace& scvf,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache)
+{
+ FluxVariables fluxVars;
+ return fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache.gridFluxVarsCache()) /scvf.area();
+}
+
+template<class FluxVariables, class Problem, class GridVariables, class SolutionVector>
+void setInterfacePressures(const Problem& problem,
+ const GridVariables& gridVars,
+ const SolutionVector& sol)
+{
+ const auto& fvGridGeometry = problem.fvGridGeometry();
+ auto fvGeometry = localView(fvGridGeometry);
+ auto elemVolVars = localView(gridVars.curGridVolVars());
+ auto elemFaceVars = localView(gridVars.curGridFaceVars());
+ auto elemFluxVarsCache = localView(gridVars.gridFluxVarsCache());
+
+ auto& couplingInterface = precice_adapter::PreciceAdapter::getInstance();
+ const auto pressureId = couplingInterface.getIdFromName( "Pressure" );
+
+ size_t i = 0;
+ for (const auto& element : elements(fvGridGeometry.gridView()))
+ {
+ fvGeometry.bind(element);
+ elemVolVars.bind(element, fvGeometry, sol);
+ elemFaceVars.bind(element, fvGeometry, sol);
+ elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
+
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+
+ if ( couplingInterface.isCoupledEntity( scvf.index() ) )
+ {
+ //TODO: What to do here?
+// const auto p = pressureAtInterface<FluxVariables>(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
+ const auto p = MonolithicSolution::pressure[i];
+ ++i;
+ couplingInterface.writeScalarQuantityOnFace( pressureId, scvf.index(), p );
+ }
+ }
+ }
+
+}
+
+
+template<class Problem, class GridVariables, class SolutionVector>
+void setInterfaceVelocities(const Problem& problem,
+ const GridVariables& gridVars,
+ const SolutionVector& sol)
+{
+ const auto& fvGridGeometry = problem.fvGridGeometry();
+ auto fvGeometry = localView(fvGridGeometry);
+ auto elemVolVars = localView(gridVars.curGridVolVars());
+ auto elemFaceVars = localView(gridVars.curGridFaceVars());
+
+ auto& couplingInterface = precice_adapter::PreciceAdapter::getInstance();
+ const auto velocityId = couplingInterface.getIdFromName( "Velocity" );
+
+ for (const auto& element : elements(fvGridGeometry.gridView()))
+ {
+ fvGeometry.bindElement(element);
+ elemVolVars.bindElement(element, fvGeometry, sol);
+ elemFaceVars.bindElement(element, fvGeometry, sol);
+
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+
+ if ( couplingInterface.isCoupledEntity( scvf.index() ) )
+ {
+ //TODO: What to do here?
+ const auto v = velocityAtInterface(elemFaceVars, scvf)[scvf.directionIndex()];
+ couplingInterface.writeScalarQuantityOnFace( velocityId, scvf.index(), v );
+ }
+ }
+ }
+
+}
+
+template<class Problem, class GridVariables, class SolutionVector>
+std::tuple<double,double,double> writeVelocitiesOnInterfaceToFile( const std::string& filename,
+ const Problem& problem,
+ const GridVariables& gridVars,
+ const SolutionVector& sol)
+{
+ const auto& fvGridGeometry = problem.fvGridGeometry();
+ auto fvGeometry = localView(fvGridGeometry);
+ auto elemVolVars = localView(gridVars.curGridVolVars());
+ auto elemFaceVars = localView(gridVars.curGridFaceVars());
+
+ const auto& couplingInterface = precice_adapter::PreciceAdapter::getInstance();
+
+ std::ofstream ofs( filename+".csv", std::ofstream::out | std::ofstream::trunc);
+ ofs << "x,y,";
+ if ( couplingInterface.getDimensions() == 3 )
+ ofs << "z,";
+ ofs << "velocityY" << "\n";
+
+
+ double min = std::numeric_limits<double>::max();
+ double max = std::numeric_limits<double>::min();
+ double sum = 0.;
+ for (const auto& element : elements(fvGridGeometry.gridView()))
+ {
+ fvGeometry.bind(element);
+ elemVolVars.bind(element, fvGeometry, sol);
+ elemFaceVars.bindElement(element, fvGeometry, sol);
+
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+
+ if ( couplingInterface.isCoupledEntity( scvf.index() ) )
+ {
+ const auto& pos = scvf.center();
+ for (int i = 0; i < couplingInterface.getDimensions(); ++i )
+ {
+ ofs << pos[i] << ",";
+ }
+ const double v = problem.dirichlet( element, scvf )[1];
+ max = std::max( v, max );
+ min = std::min( v, min );
+ sum += v;
+ //velocityAtInterface(elemFaceVars, scvf)[scvf.directionIndex()];
+ const int prec = ofs.precision();
+ ofs << std::setprecision(std::numeric_limits<double>::digits10 + 1) << v << "\n";
+ ofs.precision( prec );
+ }
+ }
+ }
+
+ ofs.close();
+
+ return std::make_tuple(min, max, sum);
+}
+
+
+
+template<class FluxVariables, class Problem, class GridVariables, class SolutionVector>
+void writePressuresOnInterfaceToFile( const std::string& filename,
+ const Problem& problem,
+ const GridVariables& gridVars,
+ const SolutionVector& sol)
+{
+ const auto& fvGridGeometry = problem.fvGridGeometry();
+ auto fvGeometry = localView(fvGridGeometry);
+ auto elemVolVars = localView(gridVars.curGridVolVars());
+ auto elemFaceVars = localView(gridVars.curGridFaceVars());
+ auto elemFluxVarsCache = localView(gridVars.gridFluxVarsCache());
+
+ const auto& couplingInterface = precice_adapter::PreciceAdapter::getInstance();
+
+ std::ofstream ofs( filename+".csv", std::ofstream::out | std::ofstream::trunc);
+ ofs << "x,y,";
+ if ( couplingInterface.getDimensions() == 3 )
+ ofs << "z,";
+ ofs << "pressure" << "\n";
+ for (const auto& element : elements(fvGridGeometry.gridView()))
+ {
+ fvGeometry.bind(element);
+ elemVolVars.bind(element, fvGeometry, sol);
+ elemFaceVars.bind(element, fvGeometry, sol);
+ elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
+
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+
+ if ( couplingInterface.isCoupledEntity( scvf.index() ) )
+ {
+ const auto& pos = scvf.center();
+ for (int i = 0; i < couplingInterface.getDimensions(); ++i )
+ {
+ ofs << pos[i] << ",";
+ }
+ const double p = pressureAtInterface<FluxVariables>(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
+ ofs << p << "\n";
+ }
+ }
+ }
+
+ ofs.close();
+}
+
+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());
+
+ // Initialize preCICE.Tell preCICE about:
+ // - Name of solver
+ // - What rank of how many ranks this instance is
+ // Configure preCICE. For now the config file is hardcoded.
+ //couplingInterface.createInstance( "FreeFlow", mpiHelper.rank(), mpiHelper.size() );
+ std::string preciceConfigFilename = "precice-config.xml";
+// if (argc == 3)
+// preciceConfigFilename = argv[2];
+ if (argc > 2)
+ preciceConfigFilename = argv[argc-1];
+
+ auto& couplingInterface =
+ precice_adapter::PreciceAdapter::getInstance();
+ couplingInterface.announceSolver( "FreeFlow", preciceConfigFilename,
+ mpiHelper.rank(), mpiHelper.size() );
+
+ const int dim = couplingInterface.getDimensions();
+ std::cout << dim << " " << int(FreeFlowFVGridGeometry::GridView::dimension) << std::endl;
+ if (dim != int(FreeFlowFVGridGeometry::GridView::dimension))
+ DUNE_THROW(Dune::InvalidStateException, "Dimensions do not match");
+
+
+ // GET mesh corodinates
+ const double xMin = getParamFromGroup<std::vector<double>>("Darcy", "Grid.Positions0")[0];
+ const double xMax = getParamFromGroup<std::vector<double>>("Darcy", "Grid.Positions0").back();
+ std::vector<double> coords; //( dim * vertexSize );
+ std::vector<int> coupledScvfIndices;
+
+ for (const auto& element : elements(freeFlowGridView))
+ {
+ auto fvGeometry = localView(*freeFlowFvGridGeometry);
+ fvGeometry.bindElement(element);
+
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+ static constexpr auto eps = 1e-7;
+ const auto& pos = scvf.center();
+ if (pos[1] < freeFlowFvGridGeometry->bBoxMin()[1] + eps)
+ {
+ if (pos[0] > xMin - eps && pos[0] < xMax + eps)
+ {
+ coupledScvfIndices.push_back(scvf.index());
+ for (const auto p : pos)
+ coords.push_back(p);
+ }
+ }
+ }
+ }
+
+ const auto numberOfPoints = coords.size() / dim;
+ const double preciceDt = couplingInterface.setMeshAndInitialize( "FreeFlowMesh",
+ numberOfPoints,
+ coords);
+ couplingInterface.createIndexMapping( coupledScvfIndices );
+
+ const auto velocityId = couplingInterface.announceQuantity( "Velocity" );
+ const auto pressureId = couplingInterface.announceQuantity( "Pressure" );
+
+ freeFlowProblem->updatePreciceDataIds();
+
+ // 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);
+
+ using FluxVariables = GetPropType<FreeFlowTypeTag, Properties::FluxVariables>;
+
+ if ( couplingInterface.hasToWriteInitialData() )
+ {
+ //TODO
+// couplingInterface.writeQuantityVector( pressureId );
+
+ setInterfacePressures<FluxVariables>( *freeFlowProblem, *freeFlowGridVariables, sol );
+ //For testing
+// {
+// std::cout << "Pressures to be sent to pm" << std::endl;
+// const auto p = couplingInterface.getQuantityVector( pressureId );
+// for (size_t i = 0; i < p.size(); ++i) {
+// std::cout << "p[" << i << "]=" <<p[i] << std::endl;
+// }
+// }
+ couplingInterface.writeScalarQuantityToOtherSolver( pressureId );
+ couplingInterface.announceInitialDataWritten();
+ }
+ couplingInterface.initializeData();
+
+ // 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);
+
+
+ auto dt = preciceDt;
+ auto sol_checkpoint = sol;
+
+ double vtkTime = 1.0;
+ size_t iter = 0;
+
+ while ( couplingInterface.isCouplingOngoing() )
+ {
+ if ( couplingInterface.hasToWriteIterationCheckpoint() )
+ {
+ //DO CHECKPOINTING
+ sol_checkpoint = sol;
+ couplingInterface.announceIterationCheckpointWritten();
+ }
+
+ // TODO
+ couplingInterface.readScalarQuantityFromOtherSolver( velocityId );
+// // For testing
+// {
+// const auto v = couplingInterface.getQuantityVector( velocityId );
+// const double sum = std::accumulate( v.begin(), v.end(), 0. );
+// std::cout << "Sum of velocities over boundary to pm: \n" << sum << std::endl;
+// }
+
+ // solve the non-linear system
+ nonLinearSolver.solve(sol);
+
+ // TODO
+ setInterfacePressures<FluxVariables>( *freeFlowProblem, *freeFlowGridVariables, sol );
+ // For testing
+// {
+// const auto p = couplingInterface.getQuantityVector( pressureId );
+// const double sum = std::accumulate( p.begin(), p.end(), 0. );
+// std::cout << "Pressures to be sent to pm" << std::endl;
+//// for (size_t i = 0; i < p.size(); ++i) {
+//// std::cout << "p[" << i << "]=" << p[i] << std::endl;
+//// }
+// std::cout << "Sum of pressures over boundary to pm: \n" << sum << std::endl;
+// }
+ couplingInterface.writeScalarQuantityToOtherSolver( pressureId );
+
+
+ //Read checkpoint
+ freeFlowVtkWriter.write(vtkTime);
+ vtkTime += 1.;
+ const double preciceDt = couplingInterface.advance( dt );
+ dt = std::min( preciceDt, dt );
+
+ //
+ {
+ double min = std::numeric_limits<double>::max();
+ double max = std::numeric_limits<double>::min();
+ double sum = 0.;
+ const std::string filename = getParam<std::string>("Problem.Name") + "-" + freeFlowProblem->name() + "-interface-velocity-" + std::to_string(iter);
+ std::tie(min, max, sum) = writeVelocitiesOnInterfaceToFile( filename,
+ *freeFlowProblem,
+ *freeFlowGridVariables,
+ sol );
+ const int prec = std::cout.precision();
+ std::cout << "Velocity statistics:" << std::endl
+ << std::setprecision(std::numeric_limits<double>::digits10 + 1)
+ << " min: " << min << std::endl
+ << " max: " << max << std::endl
+ << " sum: " << sum << std::endl;
+ std::cout.precision( prec );
+ {
+ const std::string filenameFlow="free-flow-statistics-" + std::to_string(iter);
+ std::ofstream ofs( filenameFlow+".txt", std::ofstream::out | std::ofstream::trunc);
+ const auto prec = ofs.precision();
+ ofs << "Velocity statistics (free flow):" << std::endl
+ << std::setprecision(std::numeric_limits<double>::digits10 + 1)
+ << " min: " << min << std::endl
+ << " max: " << max << std::endl
+ << " sum: " << sum << std::endl;
+ ofs.precision( prec );
+ ofs.close();
+ }
+ }
+ {
+ const std::string filename = getParam<std::string>("Problem.Name") + "-" + freeFlowProblem->name() + "-interface-pressure-" + std::to_string(iter);
+ writePressuresOnInterfaceToFile<FluxVariables>( filename,
+ *freeFlowProblem,
+ *freeFlowGridVariables,
+ sol );
+ }
+ ++iter;
+
+ if ( couplingInterface.hasToReadIterationCheckpoint() )
+ {
+// //Read checkpoint
+// freeFlowVtkWriter.write(vtkTime);
+// vtkTime += 1.;
+ sol = sol_checkpoint;
+ freeFlowGridVariables->update(sol);
+ freeFlowGridVariables->advanceTimeStep();
+ //freeFlowGridVariables->init(sol);
+ couplingInterface.announceIterationCheckpointRead();
+ }
+ else // coupling successful
+ {
+ // write vtk output
+ freeFlowVtkWriter.write(vtkTime);
+ }
+
+ }
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ {
+ double min = std::numeric_limits<double>::max();
+ double max = std::numeric_limits<double>::min();
+ double sum = 0.;
+ const std::string filename = getParam<std::string>("Problem.Name") + "-" + freeFlowProblem->name() + "-interface-velocity";
+ std::tie(min, max, sum) = writeVelocitiesOnInterfaceToFile( filename,
+ *freeFlowProblem,
+ *freeFlowGridVariables,
+ sol );
+ const int prec = std::cout.precision();
+ std::cout << "Velocity statistics:" << std::endl
+ << std::setprecision(std::numeric_limits<double>::digits10 + 1)
+ << " min: " << min << std::endl
+ << " max: " << max << std::endl
+ << " sum: " << sum << std::endl;
+ std::cout.precision( prec );
+ {
+ const std::string filenameFlow="free-flow-statistics";
+ std::ofstream ofs( filenameFlow+".txt", std::ofstream::out | std::ofstream::trunc);
+ const auto prec = ofs.precision();
+ ofs << "Velocity statistics (free flow):" << std::endl
+ << std::setprecision(std::numeric_limits<double>::digits10 + 1)
+ << " min: " << min << std::endl
+ << " max: " << max << std::endl
+ << " sum: " << sum << std::endl;
+ ofs.precision( prec );
+ ofs.close();
+ }
+ }
+ {
+ const std::string filename = getParam<std::string>("Problem.Name") + "-" + freeFlowProblem->name() + "-interface-pressure";
+ writePressuresOnInterfaceToFile<FluxVariables>( filename,
+ *freeFlowProblem,
+ *freeFlowGridVariables,
+ sol );
+ }
+
+ couplingInterface.finalize();
+
+ // 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;
+}
diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux/main_pm-reversed.cc b/appl/coupling-ff-pm/iterative-reversed-mono-flux/main_pm-reversed.cc
new file mode 100644
index 0000000000000000000000000000000000000000..e383d41997c0315f83b291b979a97bc2edbf8284
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux/main_pm-reversed.cc
@@ -0,0 +1,637 @@
+// -*- 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 <fstream>
+#include <iostream>
+#include <string>
+
+bool printstuff = false;
+
+#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/math.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 "pmproblem-reversed.hh"
+
+#include "../../precice-adapter/include/preciceadapter.hh"
+
+#include "monolithicdata.hh"
+
+
+ /*!
+ * \brief Returns the pressure at the interface using Darcy's law for reconstruction
+ */
+ template<class Problem, class Element, class FVElementGeometry, class ElementVolumeVariables, class SubControlVolumeFace>
+ auto pressureAtInterface(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf)
+ {
+ using Scalar = typename ElementVolumeVariables::VolumeVariables::PrimaryVariables::value_type;
+ const auto& volVars = elemVolVars[scvf.insideScvIdx()];
+
+ const Scalar boundaryFlux = problem.neumann(element, fvGeometry, elemVolVars, scvf);
+
+ const auto K = volVars.permeability();
+ const Scalar ccPressure = volVars.pressure();
+ const Scalar mobility = volVars.mobility();
+ const Scalar density = volVars.density();
+
+ // v = -K/mu * (gradP + rho*g)
+ auto velocity = scvf.unitOuterNormal();
+ velocity *= boundaryFlux; // TODO check sign
+ velocity /= density;
+
+ // v = -kr/mu*K * (gradP + rho*g) = -mobility*K * (gradP + rho*g)
+ const auto alpha = Dumux::vtmv( scvf.unitOuterNormal(), K, problem.gravity() );
+
+ auto distanceVector = scvf.center() - element.geometry().center();
+ distanceVector /= distanceVector.two_norm2();
+ const Scalar ti = Dumux::vtmv(distanceVector, K, scvf.unitOuterNormal());
+
+// return (1/mobility * (scvf.unitOuterNormal() * velocity) + density * alpha)/ti
+// + ccPressure;
+
+ const Scalar p = problem.dirichlet( element, scvf );
+ return p;
+ }
+
+ template<class Problem, class GridVariables, class SolutionVector>
+ void setInterfacePressures(const Problem& problem,
+ const GridVariables& gridVars,
+ const SolutionVector& sol)
+ {
+ const auto& fvGridGeometry = problem.fvGridGeometry();
+ auto fvGeometry = localView(fvGridGeometry);
+ auto elemVolVars = localView(gridVars.curGridVolVars());
+
+ auto& couplingInterface = precice_adapter::PreciceAdapter::getInstance();
+ const auto pressureId = couplingInterface.getIdFromName( "Pressure" );
+
+ for (const auto& element : elements(fvGridGeometry.gridView()))
+ {
+ fvGeometry.bindElement(element);
+ elemVolVars.bindElement(element, fvGeometry, sol);
+
+ //sstd::cout << "Pressure by reconstruction" << std::endl;
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+
+ if ( couplingInterface.isCoupledEntity( scvf.index() ) )
+ {
+ //TODO: What to do here?
+ const double p = pressureAtInterface(problem, element, fvGridGeometry, elemVolVars, scvf);
+ couplingInterface.writeScalarQuantityOnFace( pressureId, scvf.index(), p );
+ }
+ }
+ }
+ }
+
+ /*!
+ * \brief Returns the velocity at the interface using Darcy's law for reconstruction
+ */
+ template<class FluxVariables, class Problem, class Element, class FVElementGeometry,
+ class ElementVolumeVariables, class SubControlVolumeFace,
+ class ElementFluxVariablesCache>
+ auto velocityAtInterface(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const ElementFluxVariablesCache& elemFluxVarsCache)
+ {
+ const int phaseIdx = 0;
+ FluxVariables fluxVars;
+ fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+ auto upwindTerm = [phaseIdx](const auto& volVars) { return volVars.mobility(phaseIdx); };
+ const auto scalarVelocity = fluxVars.advectiveFlux(phaseIdx, upwindTerm)/scvf.area();
+ return scalarVelocity;
+ }
+
+ template<class FluxVariables, class Problem, class GridVariables, class SolutionVector>
+ void setInterfaceVelocities(const Problem& problem,
+ const GridVariables& gridVars,
+ const SolutionVector& sol)
+ {
+ const auto& fvGridGeometry = problem.fvGridGeometry();
+ auto fvGeometry = localView(fvGridGeometry);
+ auto elemVolVars = localView(gridVars.curGridVolVars());
+ auto elemFluxVarsCache = localView(gridVars.gridFluxVarsCache());
+
+ auto& couplingInterface = precice_adapter::PreciceAdapter::getInstance();
+ const auto velocityId = couplingInterface.getIdFromName( "Velocity" );
+
+
+ size_t i = 0;
+ for (const auto& element : elements(fvGridGeometry.gridView()))
+ {
+ fvGeometry.bind(element);
+ elemVolVars.bind(element, fvGeometry, sol);
+ elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
+
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+
+ if ( couplingInterface.isCoupledEntity( scvf.index() ) )
+ {
+ //TODO: What to do here?
+ //const double v = velocityAtInterface<FluxVariables>(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+ const double v = MonolithicSolution::velocity[i];
+ ++i;
+ couplingInterface.writeScalarQuantityOnFace( velocityId, scvf.index(), v );
+ }
+ }
+ }
+
+ }
+
+
+ template<class FluxVariables, class Problem, class GridVariables, class SolutionVector>
+ std::tuple<double,double,double> writeVelocitiesOnInterfaceToFile( const std::string& filename,
+ const Problem& problem,
+ const GridVariables& gridVars,
+ const SolutionVector& sol)
+ {
+ const auto& fvGridGeometry = problem.fvGridGeometry();
+ auto fvGeometry = localView(fvGridGeometry);
+ auto elemVolVars = localView(gridVars.curGridVolVars());
+ auto elemFluxVarsCache = localView(gridVars.gridFluxVarsCache());
+
+ const auto& couplingInterface = precice_adapter::PreciceAdapter::getInstance();
+
+ std::ofstream ofs( filename+".csv", std::ofstream::out | std::ofstream::trunc);
+ ofs << "x,y,";
+ if ( couplingInterface.getDimensions() == 3 )
+ ofs << "z,";
+ ofs << "velocityY" << "\n";
+
+ double min = std::numeric_limits<double>::max();
+ double max = std::numeric_limits<double>::min();
+ double sum = 0.;
+ for (const auto& element : elements(fvGridGeometry.gridView()))
+ {
+ fvGeometry.bind(element);
+ elemVolVars.bind(element, fvGeometry, sol);
+ elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
+
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+
+ if ( couplingInterface.isCoupledEntity( scvf.index() ) )
+ {
+ const auto& pos = scvf.center();
+ for (int i = 0; i < couplingInterface.getDimensions(); ++i )
+ {
+ ofs << pos[i] << ",";
+ }
+ const double v = velocityAtInterface<FluxVariables>(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+ max = std::max( v, max );
+ min = std::min( v, min );
+ sum += v;
+ const int prec = ofs.precision();
+ ofs << std::setprecision(std::numeric_limits<double>::digits10 + 1) << v << "\n";
+ ofs.precision( prec );
+ }
+ }
+ }
+
+ ofs.close();
+ return std::make_tuple(min, max, sum);
+ }
+
+
+
+ template<class Problem, class GridVariables, class SolutionVector>
+ void writePressuresOnInterfaceToFile( const std::string& filename,
+ const Problem& problem,
+ const GridVariables& gridVars,
+ const SolutionVector& sol)
+ {
+ const auto& fvGridGeometry = problem.fvGridGeometry();
+ auto fvGeometry = localView(fvGridGeometry);
+ auto elemVolVars = localView(gridVars.curGridVolVars());
+ auto elemFluxVarsCache = localView(gridVars.gridFluxVarsCache());
+
+ const auto& couplingInterface = precice_adapter::PreciceAdapter::getInstance();
+
+ std::ofstream ofs( filename+".csv", std::ofstream::out | std::ofstream::trunc);
+ ofs << "x,y,";
+ if ( couplingInterface.getDimensions() == 3 )
+ ofs << "z,";
+ ofs << "pressure" << "\n";
+ for (const auto& element : elements(fvGridGeometry.gridView()))
+ {
+ fvGeometry.bind(element);
+ elemVolVars.bind(element, fvGeometry, sol);
+ elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
+
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+
+ if ( couplingInterface.isCoupledEntity( scvf.index() ) )
+ {
+ const auto& pos = scvf.center();
+ for (int i = 0; i < couplingInterface.getDimensions(); ++i )
+ {
+ ofs << pos[i] << ",";
+ }
+ const double p = pressureAtInterface(problem, element, fvGridGeometry, elemVolVars, scvf);
+ const int prec = ofs.precision();
+ ofs << std::setprecision(std::numeric_limits<double>::digits10 + 1);
+ ofs << p << "\n";
+ ofs.precision( prec );
+ }
+ }
+ }
+
+ ofs.close();
+ }
+
+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());
+
+ // Initialize preCICE.Tell preCICE about:
+ // - Name of solver
+ // - What rank of how many ranks this instance is
+ // Configure preCICE. For now the config file is hardcoded.
+ //couplingInterface.createInstance( "darcy", mpiHelper.rank(), mpiHelper.size() );
+ std::string preciceConfigFilename = "precice-config.xml";
+// if (argc == 3)
+// preciceConfigFilename = argv[2];
+ if (argc > 2)
+ preciceConfigFilename = argv[argc-1];
+
+ auto& couplingInterface =
+ precice_adapter::PreciceAdapter::getInstance();
+ couplingInterface.announceSolver( "Darcy", preciceConfigFilename,
+ mpiHelper.rank(), mpiHelper.size() );
+
+ const int dim = couplingInterface.getDimensions();
+ std::cout << dim << " " << int(DarcyFVGridGeometry::GridView::dimension) << std::endl;
+ if (dim != int(DarcyFVGridGeometry::GridView::dimension))
+ DUNE_THROW(Dune::InvalidStateException, "Dimensions do not match");
+
+ // GET mesh corodinates
+ const double xMin = getParamFromGroup<std::vector<double>>("Darcy", "Grid.Positions0")[0];
+ const double xMax = getParamFromGroup<std::vector<double>>("Darcy", "Grid.Positions0").back();
+ std::vector<double> coords; //( dim * vertexSize );
+ std::vector<int> coupledScvfIndices;
+
+ for (const auto& element : elements(darcyGridView))
+ {
+ auto fvGeometry = localView(*darcyFvGridGeometry);
+ fvGeometry.bindElement(element);
+
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+ static constexpr auto eps = 1e-7;
+ const auto& pos = scvf.center();
+ if (pos[1] > darcyFvGridGeometry->bBoxMax()[1] - eps)
+ {
+ if (pos[0] > xMin - eps && pos[0] < xMax + eps)
+ {
+ coupledScvfIndices.push_back(scvf.index());
+ for (const auto p : pos)
+ coords.push_back(p);
+ }
+ }
+ }
+ }
+
+ const auto numberOfPoints = coords.size() / dim;
+ const double preciceDt = couplingInterface.setMeshAndInitialize( "DarcyMesh",
+ numberOfPoints,
+ coords);
+ couplingInterface.createIndexMapping( coupledScvfIndices );
+
+ const auto velocityId = couplingInterface.announceQuantity( "Velocity" );
+ const auto pressureId = couplingInterface.announceQuantity( "Pressure" );
+
+ darcyProblem->updatePreciceDataIds();
+
+ 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, darcyProblem->name());
+ using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
+ darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
+ GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
+ darcyVtkWriter.write(0.0);
+
+
+
+ using FluxVariables = GetPropType<DarcyTypeTag, Properties::FluxVariables>;
+ if ( couplingInterface.hasToWriteInitialData() )
+ {
+ //TODO
+ //couplingInterface.writeQuantityVector(velocityId);
+ setInterfaceVelocities<FluxVariables>( *darcyProblem, *darcyGridVariables, sol );
+ // For testing
+ {
+ const auto v = couplingInterface.getQuantityVector( velocityId );
+ std::cout << "velocities to be sent to ff" << std::endl;
+ for (size_t i = 0; i < v.size(); ++i) {
+ for (size_t d = 0; d < dim; ++d )
+ {
+ std::cout << coords[i*dim+d] << " ";
+ }
+ std::cout << "| v[" << i << "]=" << v[i] << std::endl;
+ }
+ }
+ couplingInterface.writeScalarQuantityToOtherSolver( velocityId );
+ couplingInterface.announceInitialDataWritten();
+ }
+ couplingInterface.initializeData();
+
+ // 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);
+
+ auto dt = preciceDt;
+ auto sol_checkpoint = sol;
+
+ double vtkTime = 1.0;
+ size_t iter = 0;
+
+ while ( couplingInterface.isCouplingOngoing() )
+ {
+ if ( couplingInterface.hasToWriteIterationCheckpoint() )
+ {
+ //DO CHECKPOINTING
+ sol_checkpoint = sol;
+ couplingInterface.announceIterationCheckpointWritten();
+ }
+
+ // TODO
+ couplingInterface.readScalarQuantityFromOtherSolver( pressureId );
+ // For testing
+ {
+ const auto p = couplingInterface.getQuantityVector( pressureId );
+ for (size_t i = 0; i < p.size(); ++i) {
+ for (size_t d = 0; d < dim; ++d )
+ {
+ std::cout << coords[i*dim+d] << " ";
+ }
+ std::cout << "| p[" << i << "]=" << p[i] << std::endl;
+ }
+ const double sum = std::accumulate( p.begin(), p.end(), 0. );
+ std::cout << "Sum of pressures over boundary to ff: \n" << sum << std::endl;
+ std::cout << "Pressure received from ff" << std::endl;
+// for (size_t i = 0; i < p.size(); ++i) {
+// std::cout << "p[" << i << "]=" << p[i] << std::endl;
+// }
+ }
+
+ // solve the non-linear system
+ nonLinearSolver.solve(sol);
+ setInterfaceVelocities<FluxVariables>( *darcyProblem, *darcyGridVariables, sol );
+ // For testing
+ {
+ const auto v = couplingInterface.getQuantityVector( velocityId );
+ for (size_t i = 0; i < v.size(); ++i) {
+ for (size_t d = 0; d < dim; ++d )
+ {
+ std::cout << coords[i*dim+d] << " ";
+ }
+ std::cout << "| v[" << i << "]=" << v[i] << std::endl;
+ }
+
+ const double sum = std::accumulate( v.begin(), v.end(), 0. );
+ std::cout << "Velocities to be sent to ff" << std::endl;
+// for (size_t i = 0; i < v.size(); ++i) {
+// std::cout << "v[" << i << "]=" << v[i] << std::endl;
+// }
+ std::cout << "Sum of velocities over boundary to ff: \n" << sum << std::endl;
+ }
+ couplingInterface.writeScalarQuantityToOtherSolver( velocityId );
+
+ const double preciceDt = couplingInterface.advance( dt );
+ dt = std::min( preciceDt, dt );
+
+ {
+ double min = std::numeric_limits<double>::max();
+ double max = std::numeric_limits<double>::min();
+ double sum = 0.;
+ const std::string filename = getParam<std::string>("Problem.Name") + "-" + darcyProblem->name() + "-interface-velocity-" + std::to_string(iter);
+ std::tie(min, max, sum) = writeVelocitiesOnInterfaceToFile<FluxVariables>( filename,
+ *darcyProblem,
+ *darcyGridVariables,
+ sol );
+ const int prec = std::cout.precision();
+ std::cout << "Velocity statistics:" << std::endl
+ << std::setprecision(std::numeric_limits<double>::digits10 + 1)
+ << " min: " << min << std::endl
+ << " max: " << max << std::endl
+ << " sum: " << sum << std::endl;
+ std::cout.precision( prec );
+ {
+ const std::string filenameFlow="darcy-statistics-" + std::to_string(iter);
+ std::ofstream ofs( filenameFlow+".txt", std::ofstream::out | std::ofstream::trunc);
+ const auto prec = ofs.precision();
+ ofs << "Velocity statistics (free flow):" << std::endl
+ << std::setprecision(std::numeric_limits<double>::digits10 + 1)
+ << " min: " << min << std::endl
+ << " max: " << max << std::endl
+ << " sum: " << sum << std::endl;
+ ofs.precision( prec );
+ ofs.close();
+ }
+ }
+ {
+ const std::string filename = getParam<std::string>("Problem.Name") + "-" + darcyProblem->name() + "-interface-pressure-" + std::to_string(iter);
+ writePressuresOnInterfaceToFile( filename,
+ *darcyProblem,
+ *darcyGridVariables,
+ sol );
+ }
+ ++iter;
+
+ if ( couplingInterface.hasToReadIterationCheckpoint() )
+ {
+ //Read checkpoint
+ darcyVtkWriter.write(vtkTime);
+ vtkTime += 1.;
+ sol = sol_checkpoint;
+ darcyGridVariables->update(sol);
+ darcyGridVariables->advanceTimeStep();
+ //darcyGridVariables->init(sol);
+ couplingInterface.announceIterationCheckpointRead();
+ }
+ else // coupling successful
+ {
+ // write vtk output
+ darcyVtkWriter.write(vtkTime);
+ }
+
+ }
+ // write vtk output
+ darcyVtkWriter.write(1.0);
+
+ {
+ double min = std::numeric_limits<double>::max();
+ double max = std::numeric_limits<double>::min();
+ double sum = 0.;
+ const std::string filename = getParam<std::string>("Problem.Name") + "-" + darcyProblem->name() + "-interface-velocity";
+ std::tie(min, max, sum) = writeVelocitiesOnInterfaceToFile<FluxVariables>( filename,
+ *darcyProblem,
+ *darcyGridVariables,
+ sol );
+ const int prec = std::cout.precision();
+ std::cout << "Velocity statistics:" << std::endl
+ << std::setprecision(std::numeric_limits<double>::digits10 + 1)
+ << " min: " << min << std::endl
+ << " max: " << max << std::endl
+ << " sum: " << sum << std::endl;
+ std::cout.precision( prec );
+ {
+ const std::string filenameFlow="darcy-statistics";
+ std::ofstream ofs( filenameFlow+".txt", std::ofstream::out | std::ofstream::trunc);
+ const auto prec = ofs.precision();
+ ofs << "Velocity statistics (free flow):" << std::endl
+ << std::setprecision(std::numeric_limits<double>::digits10 + 1)
+ << " min: " << min << std::endl
+ << " max: " << max << std::endl
+ << " sum: " << sum << std::endl;
+ ofs.precision( prec );
+ ofs.close();
+ }
+ }
+ {
+ const std::string filename = getParam<std::string>("Problem.Name") + "-" + darcyProblem->name() + "-interface-pressure";
+ writePressuresOnInterfaceToFile( filename,
+ *darcyProblem,
+ *darcyGridVariables,
+ sol );
+ }
+
+ couplingInterface.finalize();
+
+ ////////////////////////////////////////////////////////////
+ // 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;
+}
diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux/params.input b/appl/coupling-ff-pm/iterative-reversed-mono-flux/params.input
new file mode 100644
index 0000000000000000000000000000000000000000..31a78808eec638e91edde7449d8fe353ff16c360
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux/params.input
@@ -0,0 +1,41 @@
+
+
+[FreeFlow.Grid]
+Verbosity = true
+Positions0 = 0.0 1.0
+Positions1 = 1.0 2.0
+Cells0 = 20
+Cells1 = 20
+Grading1 = 1
+
+[Darcy.Grid]
+Verbosity = true
+Positions0 = 0.0 1.0
+Positions1 = 0.0 1.0
+Cells0 = 20
+Cells1 = 20
+Grading1 = 1
+
+[FreeFlow.Problem]
+Name = stokes-iterative
+EnableInertiaTerms = false
+#Name = navier-stokes-iterative
+#EnableInertiaTerms = true
+PressureDifference = 1e-9
+
+
+[Darcy.Problem]
+Name = darcy-iterative
+InitialP = 0.0e-9
+
+[Darcy.SpatialParams]
+Permeability = 1e-6 # m^2
+Porosity = 0.4
+AlphaBeaversJoseph = 1.0
+
+[Problem]
+Name = ex_ff-pm-interface
+EnableGravity = false
+
+[Vtk]
+AddVelocity = 1
diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux/pmproblem-reversed.hh b/appl/coupling-ff-pm/iterative-reversed-mono-flux/pmproblem-reversed.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ef1830459d9808f5f53d8c89b77e4f33b29ae960
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux/pmproblem-reversed.hh
@@ -0,0 +1,320 @@
+// -*- 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 The porous medium flow sub problem
+ */
+#ifndef DUMUX_DARCY_SUBPROBLEM_HH
+#define DUMUX_DARCY_SUBPROBLEM_HH
+
+#ifndef ENABLEMONOLITHIC
+#define ENABLEMONOLITHIC 0
+#endif
+
+#include <dune/grid/yaspgrid.hh>
+
+//****** uncomment for the last exercise *****//
+// #include <dumux/io/grid/subgridgridcreator.hh>
+
+#include <dumux/discretization/cctpfa.hh>
+
+#include <dumux/porousmediumflow/1p/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+
+#include "1pspatialparams.hh"
+
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/fluidsystems/1pliquid.hh>
+
+#include "../../precice-adapter/include/preciceadapter.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class DarcySubProblem;
+
+namespace Properties
+{
+// Create new type tags
+namespace TTag {
+struct DarcyOneP { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; };
+} // end namespace TTag
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DarcyOneP> { using type = Dumux::DarcySubProblem<TypeTag>; };
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DarcyOneP>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
+};
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::DarcyOneP>
+{
+ static constexpr auto dim = 2;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using TensorGrid = Dune::YaspGrid<2, Dune::TensorProductCoordinates<Scalar, dim> >;
+
+//****** comment out for the last exercise *****//
+ using type = TensorGrid;
+
+//****** uncomment for the last exercise *****//
+ // using HostGrid = TensorGrid;
+ // using type = Dune::SubGrid<dim, HostGrid>;
+};
+
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyOneP> {
+ using type = OnePSpatialParams<GetPropType<TypeTag, FVGridGeometry>, GetPropType<TypeTag, Scalar>>;
+};
+
+} // end namespace Properties
+
+/*!
+ * \brief The porous medium flow sub problem
+ */
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+ using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+ using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+
+ 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),
+ couplingInterface_(precice_adapter::PreciceAdapter::getInstance() ),
+ pressureId_(0),
+ velocityId_(0),
+ dataIdsWereSet_(false)
+#endif
+ {}
+
+ /*!
+ * \name Simulation steering
+ */
+ // \{
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ *
+ */
+ Scalar temperature() const
+ { return 273.15 + 10; } // 10°C
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary control volume.
+ *
+ * \param element The element
+ * \param scvf The boundary sub control volume face
+ */
+ BoundaryTypes boundaryTypes(const Element &element, const SubControlVolumeFace &scvf) const
+ {
+ BoundaryTypes values;
+
+ // 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();
+#else
+ const auto faceId = scvf.index();
+ if ( couplingInterface_.isCoupledEntity(faceId) )
+ values.setAllDirichlet();
+#endif
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The element for which the Dirichlet boundary condition is set
+ * \param scvf The boundary subcontrolvolumeface
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
+ {
+ // set p = 0 at the bottom
+ PrimaryVariables values(0.0);
+ values = initial(element);
+
+ const auto faceId = scvf.index();
+ if ( couplingInterface_.isCoupledEntity(faceId) )
+ values = couplingInterface_.getScalarQuantityOnFace( pressureId_, faceId );
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scvf The boundary sub control volume face
+ *
+ * For this method, the \a values variable stores primary variables.
+ */
+ template<class ElementVolumeVariables>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ // 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
+// assert( dataIdsWereSet_ );
+// const auto faceId = scvf.index();
+// if ( couplingInterface_.isCoupledEntity(faceId) )
+// {
+// const Scalar density = 1000.;
+// values[Indices::conti0EqIdx] = density * couplingInterface_.getScalarQuantityOnFace( velocityId_, faceId );
+// std::cout << "pm: values[Indices::conti0EqIdx] = " << values << std::endl;
+// }
+#endif
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * \param element The element for which the source term is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scv The subcontrolvolume
+ */
+ template<class ElementVolumeVariables>
+ NumEqVector source(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ { return NumEqVector(0.0); }
+
+ // \}
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param element The element
+ *
+ * For this method, the \a priVars parameter stores primary
+ * variables.
+ */
+ PrimaryVariables initial(const Element &element) const
+ {
+ static const Scalar p = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.InitialP");
+ return PrimaryVariables(p);
+ }
+
+ // \}
+
+#if !ENABLEMONOLITHIC
+ void updatePreciceDataIds()
+ {
+ pressureId_ = couplingInterface_.getIdFromName( "Pressure" );
+ velocityId_ = couplingInterface_.getIdFromName( "Velocity" );
+ dataIdsWereSet_ = true;
+ }
+#endif
+
+#if ENABLEMONOLITHIC
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+#endif
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ Scalar eps_;
+
+#if ENABLEMONOLITHIC
+ std::shared_ptr<CouplingManager> couplingManager_;
+#else
+ precice_adapter::PreciceAdapter& couplingInterface_;
+ size_t pressureId_;
+ size_t velocityId_;
+ bool dataIdsWereSet_;
+
+#endif
+};
+} //end namespace
+
+#endif //DUMUX_DARCY_SUBPROBLEM_HH
diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux/precice-config-parallel-implicit-reversed.xml b/appl/coupling-ff-pm/iterative-reversed-mono-flux/precice-config-parallel-implicit-reversed.xml
new file mode 100644
index 0000000000000000000000000000000000000000..b0ff6c0d2b15e16dc691c1c9b79816cfa9997aee
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux/precice-config-parallel-implicit-reversed.xml
@@ -0,0 +1,79 @@
+<?xml version="1.0"?>
+
+<precice-configuration>
+ <log>
+ <sink type="stream" output="stdout" filter= "(%Severity% > debug) or (%Severity% >= trace and %Module% contains SolverInterfaceImpl)" enabled="false" />
+ <sink type="stream" output="stdout" enabled="false" />
+ </log>
+
+ <solver-interface dimensions="2">
+
+ <data:scalar name="Velocity"/>
+ <data:scalar name="Pressure"/>
+
+ <mesh name="FreeFlowMesh">
+ <use-data name="Velocity" />
+ <use-data name="Pressure" />
+ </mesh>
+
+ <mesh name="DarcyMesh">
+ <use-data name="Velocity" />
+ <use-data name="Pressure" />
+ </mesh>
+
+ <participant name="FreeFlow">
+ <use-mesh name="FreeFlowMesh" provide="yes"/>
+ <use-mesh name="DarcyMesh" from="Darcy"/>
+
+ <read-data name="Velocity" mesh="FreeFlowMesh"/>
+ <write-data name="Pressure" mesh="FreeFlowMesh"/>
+
+ <mapping:nearest-neighbor direction="write" from="FreeFlowMesh" to="DarcyMesh" constraint="consistent"/>
+ <mapping:nearest-neighbor direction="read" from="DarcyMesh" to="FreeFlowMesh" constraint="consistent"/>
+
+ </participant>
+
+ <participant name="Darcy">
+ <use-mesh name="DarcyMesh" provide="yes"/>
+
+ <read-data name="Pressure" mesh="DarcyMesh"/>
+ <write-data name="Velocity" mesh="DarcyMesh"/>
+ </participant>
+
+ <m2n:sockets from="FreeFlow" to="Darcy" distribution-type="gather-scatter" network="lo" />
+
+
+ <coupling-scheme:parallel-implicit>
+ <max-time value="1"/>
+ <timestep-length value="1" />
+ <max-iterations value="100"/>
+
+
+ <participants first="FreeFlow" second="Darcy"/>
+ <exchange data="Pressure" mesh="DarcyMesh" from="FreeFlow" to="Darcy" initialize="true" />
+ <exchange data="Velocity" mesh="DarcyMesh" from="Darcy" to="FreeFlow" initialize="true" />
+
+ <relative-convergence-measure limit="1.0e-5" data="Pressure" mesh="DarcyMesh"/>
+ <relative-convergence-measure limit="1.0e-5" data="Velocity" mesh="DarcyMesh"/>
+
+
+ <!--
+ <relative-convergence-measure limit="1.0e-2" data="Velocity" mesh="FreeFlowMesh"/>
+ -->
+
+ <extrapolation-order value="0"/>
+
+ <post-processing:IQN-ILS>
+ <data mesh="DarcyMesh" name="Velocity" />
+ <data mesh="DarcyMesh" name="Pressure" />
+ <preconditioner type="residual-sum" />
+ <initial-relaxation value="0.1" />
+ <max-used-iterations value="9" />
+ <timesteps-reused value="10" />
+ <filter type="QR2" limit="1e-3" />
+ </post-processing:IQN-ILS>
+
+ </coupling-scheme:parallel-implicit>
+ </solver-interface>
+</precice-configuration>
+
diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux/precice-config-serial-explicit-reversed.xml b/appl/coupling-ff-pm/iterative-reversed-mono-flux/precice-config-serial-explicit-reversed.xml
new file mode 100644
index 0000000000000000000000000000000000000000..d78a37837f8b6f765a5b691c8ca0caf1dbff79b7
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux/precice-config-serial-explicit-reversed.xml
@@ -0,0 +1,57 @@
+<?xml version="1.0"?>
+
+<precice-configuration>
+ <log>
+ <sink type="stream" output="stdout" filter= "(%Severity% > debug) or (%Severity% >= trace and %Module% contains SolverInterfaceImpl)" enabled="false" />
+ <sink type="stream" output="stdout" enabled="false" />
+ </log>
+
+ <solver-interface dimensions="2">
+
+ <data:scalar name="Velocity"/>
+ <data:scalar name="Pressure"/>
+
+ <mesh name="FreeFlowMesh">
+ <use-data name="Velocity" />
+ <use-data name="Pressure" />
+ </mesh>
+
+ <mesh name="DarcyMesh">
+ <use-data name="Velocity" />
+ <use-data name="Pressure" />
+ </mesh>
+
+ <participant name="FreeFlow">
+ <use-mesh name="FreeFlowMesh" provide="yes"/>
+ <use-mesh name="DarcyMesh" from="Darcy"/>
+
+ <read-data name="Velocity" mesh="FreeFlowMesh"/>
+ <write-data name="Pressure" mesh="FreeFlowMesh"/>
+
+ <mapping:nearest-neighbor direction="write" from="FreeFlowMesh" to="DarcyMesh" constraint="consistent"/>
+ <mapping:nearest-neighbor direction="read" from="DarcyMesh" to="FreeFlowMesh" constraint="consistent"/>
+
+ </participant>
+
+ <participant name="Darcy">
+ <use-mesh name="DarcyMesh" provide="yes"/>
+
+ <read-data name="Pressure" mesh="DarcyMesh"/>
+ <write-data name="Velocity" mesh="DarcyMesh"/>
+ </participant>
+
+ <m2n:sockets from="FreeFlow" to="Darcy" distribution-type="gather-scatter" network="lo" />
+
+
+ <coupling-scheme:serial-explicit>
+ <max-time value="1"/>
+ <timestep-length value="1" />
+
+ <participants first="FreeFlow" second="Darcy"/>
+ <exchange data="Pressure" mesh="DarcyMesh" from="FreeFlow" to="Darcy" initialize="false" />
+ <exchange data="Velocity" mesh="DarcyMesh" from="Darcy" to="FreeFlow" initialize="true" />
+
+ </coupling-scheme:serial-explicit>
+ </solver-interface>
+</precice-configuration>
+
diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux/precice-config-serial-implicit-reversed-darcy-first.xml b/appl/coupling-ff-pm/iterative-reversed-mono-flux/precice-config-serial-implicit-reversed-darcy-first.xml
new file mode 100644
index 0000000000000000000000000000000000000000..55a904e532fe7df4c702ff7cc9ed923ec9f14bd2
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux/precice-config-serial-implicit-reversed-darcy-first.xml
@@ -0,0 +1,79 @@
+<?xml version="1.0"?>
+
+<precice-configuration>
+ <log>
+ <sink type="stream" output="stdout" filter= "(%Severity% > debug) or (%Severity% >= trace and %Module% contains SolverInterfaceImpl)" enabled="false" />
+ <sink type="stream" output="stdout" enabled="false" />
+ </log>
+
+ <solver-interface dimensions="2">
+
+ <data:scalar name="Velocity"/>
+ <data:scalar name="Pressure"/>
+
+ <mesh name="FreeFlowMesh">
+ <use-data name="Velocity" />
+ <use-data name="Pressure" />
+ </mesh>
+
+ <mesh name="DarcyMesh">
+ <use-data name="Velocity" />
+ <use-data name="Pressure" />
+ </mesh>
+
+ <participant name="FreeFlow">
+ <use-mesh name="FreeFlowMesh" provide="yes"/>
+ <use-mesh name="DarcyMesh" from="Darcy"/>
+
+ <read-data name="Velocity" mesh="FreeFlowMesh"/>
+ <write-data name="Pressure" mesh="FreeFlowMesh"/>
+
+ <mapping:nearest-neighbor direction="write" from="FreeFlowMesh" to="DarcyMesh" constraint="consistent"/>
+ <mapping:nearest-neighbor direction="read" from="DarcyMesh" to="FreeFlowMesh" constraint="consistent"/>
+
+ </participant>
+
+ <participant name="Darcy">
+ <use-mesh name="DarcyMesh" provide="yes"/>
+
+ <read-data name="Pressure" mesh="DarcyMesh"/>
+ <write-data name="Velocity" mesh="DarcyMesh"/>
+ </participant>
+
+ <m2n:sockets from="FreeFlow" to="Darcy" distribution-type="gather-scatter" network="lo" />
+
+
+ <coupling-scheme:serial-implicit>
+ <max-time value="1"/>
+ <timestep-length value="1" />
+ <max-iterations value="10"/>
+
+
+ <participants first="Darcy" second="FreeFlow"/>
+ <exchange data="Pressure" mesh="DarcyMesh" from="FreeFlow" to="Darcy" initialize="true" />
+ <exchange data="Velocity" mesh="DarcyMesh" from="Darcy" to="FreeFlow" initialize="false" />
+
+ <relative-convergence-measure limit="1.0e-4" data="Pressure" mesh="DarcyMesh"/>
+ <relative-convergence-measure limit="1.0e-4" data="Velocity" mesh="DarcyMesh"/>
+
+
+ <!--
+ <relative-convergence-measure limit="1.0e-2" data="Velocity" mesh="FreeFlowMesh"/>
+ -->
+
+
+ <extrapolation-order value="0"/>
+<!--
+ <post-processing:IQN-ILS>
+ <data mesh="DarcyMesh" name="Velocity" />
+ <initial-relaxation value="0.1" />
+ <max-used-iterations value="9" />
+ <timesteps-reused value="10" />
+ <filter type="QR2" limit="1e-3" />
+ </post-processing:IQN-ILS>
+-->
+
+ </coupling-scheme:serial-implicit>
+ </solver-interface>
+</precice-configuration>
+
diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux/precice-config-serial-implicit-reversed.xml b/appl/coupling-ff-pm/iterative-reversed-mono-flux/precice-config-serial-implicit-reversed.xml
new file mode 100644
index 0000000000000000000000000000000000000000..756b2b07b6a895b55715740a40e5118b9614daa0
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux/precice-config-serial-implicit-reversed.xml
@@ -0,0 +1,79 @@
+<?xml version="1.0"?>
+
+<precice-configuration>
+ <log>
+ <sink type="stream" output="stdout" filter= "(%Severity% > debug) or (%Severity% >= trace and %Module% contains SolverInterfaceImpl)" enabled="false" />
+ <sink type="stream" output="stdout" enabled="false" />
+ </log>
+
+ <solver-interface dimensions="2">
+
+ <data:scalar name="Velocity"/>
+ <data:scalar name="Pressure"/>
+
+ <mesh name="FreeFlowMesh">
+ <use-data name="Velocity" />
+ <use-data name="Pressure" />
+ </mesh>
+
+ <mesh name="DarcyMesh">
+ <use-data name="Velocity" />
+ <use-data name="Pressure" />
+ </mesh>
+
+ <participant name="FreeFlow">
+ <use-mesh name="FreeFlowMesh" provide="yes"/>
+ <use-mesh name="DarcyMesh" from="Darcy"/>
+
+ <read-data name="Velocity" mesh="FreeFlowMesh"/>
+ <write-data name="Pressure" mesh="FreeFlowMesh"/>
+
+ <mapping:nearest-neighbor direction="write" from="FreeFlowMesh" to="DarcyMesh" constraint="consistent"/>
+ <mapping:nearest-neighbor direction="read" from="DarcyMesh" to="FreeFlowMesh" constraint="consistent"/>
+
+ </participant>
+
+ <participant name="Darcy">
+ <use-mesh name="DarcyMesh" provide="yes"/>
+
+ <read-data name="Pressure" mesh="DarcyMesh"/>
+ <write-data name="Velocity" mesh="DarcyMesh"/>
+ </participant>
+
+ <m2n:sockets from="FreeFlow" to="Darcy" distribution-type="gather-scatter" network="lo" />
+
+
+ <coupling-scheme:serial-implicit>
+ <max-time value="1"/>
+ <timestep-length value="1" />
+ <max-iterations value="10"/>
+
+
+ <participants first="FreeFlow" second="Darcy"/>
+ <exchange data="Pressure" mesh="DarcyMesh" from="FreeFlow" to="Darcy" initialize="false" />
+ <exchange data="Velocity" mesh="DarcyMesh" from="Darcy" to="FreeFlow" initialize="true" />
+
+ <relative-convergence-measure limit="1.0e-4" data="Pressure" mesh="DarcyMesh"/>
+ <relative-convergence-measure limit="1.0e-4" data="Velocity" mesh="DarcyMesh"/>
+
+
+ <!--
+ <relative-convergence-measure limit="1.0e-2" data="Velocity" mesh="FreeFlowMesh"/>
+ -->
+
+
+ <extrapolation-order value="0"/>
+<!--
+ <post-processing:IQN-ILS>
+ <data mesh="DarcyMesh" name="Velocity" />
+ <initial-relaxation value="0.1" />
+ <max-used-iterations value="9" />
+ <timesteps-reused value="10" />
+ <filter type="QR2" limit="1e-3" />
+ </post-processing:IQN-ILS>
+-->
+
+ </coupling-scheme:serial-implicit>
+ </solver-interface>
+</precice-configuration>
+