From ad573a24bbb1495361832862948b13e7201d696a Mon Sep 17 00:00:00 2001
From: Alexander Jaust <alexander.jaust@ipvs.uni-stuttgart.de>
Date: Tue, 21 Apr 2020 15:06:08 +0200
Subject: [PATCH] [WIP] initial commit for solver using monolithic flux
---
.../1pspatialparams.hh | 91 +++
.../CMakeLists.txt | 37 +
.../ffproblem-reversed.hh | 414 ++++++++++++
.../main_ff-reversed.cc | 584 ++++++++++++++++
.../main_pm-reversed.cc | 637 ++++++++++++++++++
.../monolithicdata.hh | 55 ++
.../params.input | 41 ++
.../pmproblem-reversed.hh | 320 +++++++++
8 files changed, 2179 insertions(+)
create mode 100644 appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/1pspatialparams.hh
create mode 100644 appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/CMakeLists.txt
create mode 100644 appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/ffproblem-reversed.hh
create mode 100644 appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/main_ff-reversed.cc
create mode 100644 appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/main_pm-reversed.cc
create mode 100644 appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/monolithicdata.hh
create mode 100644 appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/params.input
create mode 100644 appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/pmproblem-reversed.hh
diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/1pspatialparams.hh b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/1pspatialparams.hh
new file mode 100644
index 0000000..a042388
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/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-no-precice/CMakeLists.txt b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/CMakeLists.txt
new file mode 100644
index 0000000..5e92710
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/CMakeLists.txt
@@ -0,0 +1,37 @@
+# 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_mono_flux
+# SOURCES main_ff-reversed.cc ../precice/preciceadapter.cc ../precice/dumuxpreciceindexwrapper.cc
+# COMPILE_DEFINITIONS ENABLEMONOLITHIC=0)
+
+#dune_add_test(NAME test_pm_reversed_mono_flux
+# SOURCES main_pm-reversed.cc ../precice/preciceadapter.cc ../precice/dumuxpreciceindexwrapper.cc
+# COMPILE_DEFINITIONS ENABLEMONOLITHIC=0)
+
+add_executable(test_ff_reversed_mono_flux-no-precice EXCLUDE_FROM_ALL main_ff-reversed.cc ../../precice-adapter/src/preciceadapter.cc ../../precice-adapter/src/dumuxpreciceindexwrapper.cc)
+add_executable(test_pm_reversed_mono_flux-no-precice EXCLUDE_FROM_ALL main_pm-reversed.cc ../../precice-adapter/src/preciceadapter.cc ../../precice-adapter/src/dumuxpreciceindexwrapper.cc)
+
+target_compile_definitions(test_ff_reversed_mono_flux-no-precice PUBLIC "ENABLEMONOLITHIC=0")
+target_compile_definitions(test_pm_reversed_mono_flux-no-precice PUBLIC "ENABLEMONOLITHIC=0")
+
+target_compile_definitions(test_ff_reversed_mono_flux-no-precice PRIVATE BOOST_ALL_DYN_LINK)
+target_link_libraries(test_ff_reversed_mono_flux-no-precice PRIVATE ${Boost_LIBRARIES} ${PRECICE_LIB})
+target_include_directories(test_ff_reversed_mono_flux-no-precice PRIVATE ${Boost_INCLUDE_DIRS})
+
+target_compile_definitions(test_pm_reversed_mono_flux-no-precice PRIVATE BOOST_ALL_DYN_LINK)
+target_link_libraries(test_pm_reversed_mono_flux-no-precice PRIVATE ${Boost_LIBRARIES} ${PRECICE_LIB})
+target_include_directories(test_pm_reversed_mono_flux-no-precice PRIVATE ${Boost_INCLUDE_DIRS})
+# add a symlink for each input file
+add_input_file_links()
+
+macro(add_precice_file_links)
+ FILE(GLOB precice_files *.xml)
+ foreach(VAR ${input_files})
+ get_filename_component(file_name ${VAR} NAME)
+ dune_symlink_to_source_files(FILES ${file_name})
+ endforeach()
+endmacro()
+add_precice_file_links()
diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/ffproblem-reversed.hh b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/ffproblem-reversed.hh
new file mode 100644
index 0000000..0a9ef2a
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/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-no-precice/main_ff-reversed.cc b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/main_ff-reversed.cc
new file mode 100644
index 0000000..847c08d
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/main_ff-reversed.cc
@@ -0,0 +1,584 @@
+// -*- 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 "monolithicdata.hh"
+
+
+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-no-precice/main_pm-reversed.cc b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/main_pm-reversed.cc
new file mode 100644
index 0000000..e383d41
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/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-no-precice/monolithicdata.hh b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/monolithicdata.hh
new file mode 100644
index 0000000..5653a59
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/monolithicdata.hh
@@ -0,0 +1,55 @@
+#ifndef MONOLITHICDATA_HH
+#define MONOLITHICDATA_HH
+
+#include <array>
+
+namespace MonolithicSolution
+{
+static 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
+};
+
+static 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 };
+}
+
+#endif // MONOLITHICDATA_HH
diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/params.input b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/params.input
new file mode 100644
index 0000000..31a7880
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/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-no-precice/pmproblem-reversed.hh b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/pmproblem-reversed.hh
new file mode 100644
index 0000000..ef18304
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/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
--
GitLab