From 741cce5b479affc358b08b270083b9652cbc02db Mon Sep 17 00:00:00 2001
From: Alexander Jaust <alexander.jaust@ipvs.uni-stuttgart.de>
Date: Fri, 5 Jul 2019 16:13:19 +0200
Subject: [PATCH] files for reverse coupling of ff-pm coupling
Change coupling to:
pm receives pressure and communicates velocity
WIP
---
.../iterative-reversed/CMakeLists.txt | 28 ++
.../iterative-reversed/ffproblem-reversed.hh | 415 ++++++++++++++++++
.../iterative-reversed/main_ff-reversed.cc | 317 +++++++++++++
.../iterative-reversed/main_pm-reversed.cc | 340 ++++++++++++++
.../iterative-reversed/params.input | 38 ++
.../iterative-reversed/pmproblem-reversed.hh | 322 ++++++++++++++
6 files changed, 1460 insertions(+)
create mode 100644 appl/coupling-ff-pm/iterative-reversed/CMakeLists.txt
create mode 100644 appl/coupling-ff-pm/iterative-reversed/ffproblem-reversed.hh
create mode 100644 appl/coupling-ff-pm/iterative-reversed/main_ff-reversed.cc
create mode 100644 appl/coupling-ff-pm/iterative-reversed/main_pm-reversed.cc
create mode 100644 appl/coupling-ff-pm/iterative-reversed/params.input
create mode 100644 appl/coupling-ff-pm/iterative-reversed/pmproblem-reversed.hh
diff --git a/appl/coupling-ff-pm/iterative-reversed/CMakeLists.txt b/appl/coupling-ff-pm/iterative-reversed/CMakeLists.txt
new file mode 100644
index 0000000..8b8a730
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed/CMakeLists.txt
@@ -0,0 +1,28 @@
+# executables for the iterative case
+string(REPLACE ":" ";" LIBRARY_DIRS $ENV{LD_LIBRARY_PATH})
+find_library(PRECICE_LIB NAMES precice HINTS ${LIBRARY_DIRS})
+find_package(Boost 1.65.1 REQUIRED COMPONENTS log system) #Require same version as preCICE
+
+#dune_add_test(NAME test_ff_reversed
+# SOURCES main_ff-reversed.cc ../precice/preciceadapter.cc ../precice/dumuxpreciceindexwrapper.cc
+# COMPILE_DEFINITIONS ENABLEMONOLITHIC=0)
+
+#dune_add_test(NAME test_pm_reversed
+# SOURCES main_pm-reversed.cc ../precice/preciceadapter.cc ../precice/dumuxpreciceindexwrapper.cc
+# COMPILE_DEFINITIONS ENABLEMONOLITHIC=0)
+
+add_executable(test_ff_reversed EXCLUDE_FROM_ALL main_ff-reversed.cc ../../precice-adapter/src/preciceadapter.cc ../../precice-adapter/src/dumuxpreciceindexwrapper.cc)
+add_executable(test_pm_reversed EXCLUDE_FROM_ALL main_pm-reversed.cc ../../precice-adapter/src/preciceadapter.cc ../../precice-adapter/src/dumuxpreciceindexwrapper.cc)
+
+target_compile_definitions(test_ff_reversed PUBLIC "ENABLEMONOLITHIC=0")
+target_compile_definitions(test_pm_reversed PUBLIC "ENABLEMONOLITHIC=0")
+
+target_compile_definitions(test_ff_reversed PRIVATE BOOST_ALL_DYN_LINK)
+target_link_libraries(test_ff_reversed PRIVATE ${Boost_LIBRARIES} ${PRECICE_LIB})
+target_include_directories(test_ff_reversed PRIVATE ${Boost_INCLUDE_DIRS})
+
+target_compile_definitions(test_pm_reversed PRIVATE BOOST_ALL_DYN_LINK)
+target_link_libraries(test_pm_reversed PRIVATE ${Boost_LIBRARIES} ${PRECICE_LIB})
+target_include_directories(test_pm_reversed PRIVATE ${Boost_INCLUDE_DIRS})
+# add a symlink for each input file
+add_input_file_links()
diff --git a/appl/coupling-ff-pm/iterative-reversed/ffproblem-reversed.hh b/appl/coupling-ff-pm/iterative-reversed/ffproblem-reversed.hh
new file mode 100644
index 0000000..5c7fd43
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed/ffproblem-reversed.hh
@@ -0,0 +1,415 @@
+// -*- 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();
+
+ // left/right wall
+ if (onRightBoundary_(globalPos) || (onLeftBoundary_(globalPos)))
+ {
+ values.setDirichlet(Indices::pressureIdx);
+ }
+ else
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+
+ // coupling interface
+#if ENABLEMONOLITHIC
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setBJS(Indices::momentumXBalanceIdx);
+ }
+#else
+ // // TODO do preCICE stuff in analogy to heat transfer
+ assert( dataIdsWereSet_ );
+ const auto faceId = scvf.index();
+ if ( couplingInterface_.isCoupledEntity(faceId) )
+ {
+ //TODO What do I want to do here?
+ values.setNeumann(Indices::conti0EqIdx);
+ values.setNeumann(Indices::momentumYBalanceIdx);
+ //values.setCouplingNeumann(Indices::conti0EqIdx);
+ //values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ //values.setBJS(Indices::momentumXBalanceIdx);
+ }
+#endif
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
+ {
+ PrimaryVariables values(0.0);
+ values = initialAtPos(globalPos);
+ 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 velocity = couplingInterface_.getScalarQuantityOnFace( velocityId_, faceId );
+ const auto& volVars = elemVolVars[scvf.insideScvIdx()];
+
+// 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 = Dumux::vtmv( scvf.unitOuterNormal(), K, this->gravity() );
+
+// auto distanceVector = scvf.center() - element.geometry().center();
+// distanceVector /= distanceVector.two_norm2();
+// const Scalar ti = Dumux::vtmv(distanceVector, K, scvf.unitOuterNormal());
+
+// const Scalar pressure = (1/mobility * (scvf.unitOuterNormal() * velocity) + density * alpha)/ti
+// + ccPressure;
+
+ values[Indices::conti0EqIdx] = velocity * density;
+// values[Indices::momentumYBalanceIdx] = scvf.directionSign() * ( pressure - 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/main_ff-reversed.cc b/appl/coupling-ff-pm/iterative-reversed/main_ff-reversed.cc
new file mode 100644
index 0000000..d0d3c6b
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed/main_ff-reversed.cc
@@ -0,0 +1,317 @@
+// -*- 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"
+
+//TODO
+// Helper function to put pressure on interface
+
+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 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 );
+ }
+ }
+ }
+
+}
+
+
+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];
+
+ 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);
+
+ if ( couplingInterface.hasToWriteInitialData() )
+ {
+ //TODO
+// couplingInterface.writeQuantityVector( pressureId );
+ setInterfacePressures( *freeFlowProblem, *freeFlowGridVariables, sol );
+ 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;
+
+ 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( *freeFlowProblem, *freeFlowGridVariables, sol );
+ couplingInterface.writeScalarQuantityToOtherSolver( pressureId );
+
+ const double preciceDt = couplingInterface.advance( dt );
+ dt = std::min( preciceDt, dt );
+
+ if ( couplingInterface.hasToReadIterationCheckpoint() )
+ {
+ //Read checkpoint
+ sol = sol_checkpoint;
+ freeFlowGridVariables->update(sol);
+ freeFlowGridVariables->advanceTimeStep();
+ //freeFlowGridVariables->init(sol);
+ couplingInterface.announceIterationCheckpointRead();
+ }
+ else // coupling successful
+ {
+ // write vtk output
+ freeFlowVtkWriter.write(1.0);
+ }
+
+ }
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ 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/main_pm-reversed.cc b/appl/coupling-ff-pm/iterative-reversed/main_pm-reversed.cc
new file mode 100644
index 0000000..a49cfbe
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed/main_pm-reversed.cc
@@ -0,0 +1,340 @@
+// -*- 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/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"
+
+ /*!
+ * \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);
+
+ // 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 = 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;
+ }
+
+ 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);
+
+ 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 );
+ }
+ }
+ }
+
+ }
+
+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];
+
+ 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);
+
+
+ if ( couplingInterface.hasToWriteInitialData() )
+ {
+ //TODO
+ //couplingInterface.writeQuantityVector(velocityId);
+ setInterfacePressures( *darcyProblem, *darcyGridVariables, sol );
+ couplingInterface.writeScalarQuantityToOtherSolver( pressureId );
+ 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;
+
+ 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 fluxes over boundary to pm: \n" << sum << std::endl;
+ }
+
+ // solve the non-linear system
+ nonLinearSolver.solve(sol);
+ setInterfacePressures( *darcyProblem, *darcyGridVariables, sol );
+ couplingInterface.writeScalarQuantityToOtherSolver( pressureId );
+
+ const double preciceDt = couplingInterface.advance( dt );
+ dt = std::min( preciceDt, dt );
+
+ if ( couplingInterface.hasToReadIterationCheckpoint() )
+ {
+ //Read checkpoint
+ sol = sol_checkpoint;
+ darcyGridVariables->update(sol);
+ darcyGridVariables->advanceTimeStep();
+ //darcyGridVariables->init(sol);
+ couplingInterface.announceIterationCheckpointRead();
+ }
+ else // coupling successful
+ {
+ // write vtk output
+ darcyVtkWriter.write(1.0);
+ }
+
+ }
+ // write vtk output
+ darcyVtkWriter.write(1.0);
+
+ 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/params.input b/appl/coupling-ff-pm/iterative-reversed/params.input
new file mode 100644
index 0000000..c0fa10a
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed/params.input
@@ -0,0 +1,38 @@
+
+
+[FreeFlow.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
+
+[FreeFlow.Problem]
+Name = stokes
+PressureDifference = 1e-9
+EnableInertiaTerms = false
+
+[Darcy.Problem]
+Name = darcy
+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/pmproblem-reversed.hh b/appl/coupling-ff-pm/iterative-reversed/pmproblem-reversed.hh
new file mode 100644
index 0000000..1dd752f
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed/pmproblem-reversed.hh
@@ -0,0 +1,322 @@
+// -*- 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();
+#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);
+
+ 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 );
+ const Scalar facePressure = couplingInterface_.getScalarQuantityOnFace( pressureId_, faceId );
+
+ 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 = Dumux::vtmv( scvf.unitOuterNormal(), K, problem.gravity() );
+
+ values[Indices::conti0EqIdx] = density * ;
+ }
+#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
--
GitLab