From c8a947857caec88d6d3d10b185e2425b6a694bd7 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Mon, 9 Jul 2018 18:50:10 +0200
Subject: [PATCH] [test][boundary] Add 1p1p darcy-darcy test with tpfa (domain
decomposition)
---
test/multidomain/CMakeLists.txt | 1 +
test/multidomain/boundary/CMakeLists.txt | 1 +
.../boundary/darcydarcy/1p_1p/CMakeLists.txt | 22 ++
.../boundary/darcydarcy/1p_1p/main.cc | 353 ++++++++++++++++++
.../boundary/darcydarcy/1p_1p/problem.hh | 163 ++++++++
.../darcydarcy/1p_1p/spatialparams.hh | 114 ++++++
.../1p_1p/test_boundary_equaldim_1p_1p.input | 34 ++
.../boundary/darcydarcy/CMakeLists.txt | 1 +
8 files changed, 689 insertions(+)
create mode 100644 test/multidomain/boundary/CMakeLists.txt
create mode 100644 test/multidomain/boundary/darcydarcy/1p_1p/CMakeLists.txt
create mode 100644 test/multidomain/boundary/darcydarcy/1p_1p/main.cc
create mode 100644 test/multidomain/boundary/darcydarcy/1p_1p/problem.hh
create mode 100644 test/multidomain/boundary/darcydarcy/1p_1p/spatialparams.hh
create mode 100644 test/multidomain/boundary/darcydarcy/1p_1p/test_boundary_equaldim_1p_1p.input
create mode 100644 test/multidomain/boundary/darcydarcy/CMakeLists.txt
diff --git a/test/multidomain/CMakeLists.txt b/test/multidomain/CMakeLists.txt
index fe1923bdc7..d6f0a545fe 100644
--- a/test/multidomain/CMakeLists.txt
+++ b/test/multidomain/CMakeLists.txt
@@ -1,2 +1,3 @@
add_subdirectory(embedded)
add_subdirectory(poromechanics)
+add_subdirectory(boundary)
diff --git a/test/multidomain/boundary/CMakeLists.txt b/test/multidomain/boundary/CMakeLists.txt
new file mode 100644
index 0000000000..51f43731d8
--- /dev/null
+++ b/test/multidomain/boundary/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory(darcydarcy)
diff --git a/test/multidomain/boundary/darcydarcy/1p_1p/CMakeLists.txt b/test/multidomain/boundary/darcydarcy/1p_1p/CMakeLists.txt
new file mode 100644
index 0000000000..510844ae3d
--- /dev/null
+++ b/test/multidomain/boundary/darcydarcy/1p_1p/CMakeLists.txt
@@ -0,0 +1,22 @@
+dune_add_test(NAME test_boundary_equaldim_1p_1p_half
+ SOURCES main.cc
+ COMPILE_DEFINITIONS DOMAINSPLIT=0
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/1ptestcc-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/md_1p_half_combined.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_boundary_equaldim_1p_1p_half test_boundary_equaldim_1p_1p.input \
+ -Problem.Name md_1p_half")
+
+dune_add_test(NAME test_boundary_equaldim_1p_1p_lens
+ SOURCES main.cc
+ COMPILE_DEFINITIONS DOMAINSPLIT=1
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMAKE_GUARD dune-subgrid_FOUND
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/1ptestcc-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/md_1p_lens_combined.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_boundary_equaldim_1p_1p_lens test_boundary_equaldim_1p_1p.input \
+ -Problem.Name md_1p_lens")
+
+dune_symlink_to_source_files(FILES "test_boundary_equaldim_1p_1p.input")
diff --git a/test/multidomain/boundary/darcydarcy/1p_1p/main.cc b/test/multidomain/boundary/darcydarcy/1p_1p/main.cc
new file mode 100644
index 0000000000..3d7be9226f
--- /dev/null
+++ b/test/multidomain/boundary/darcydarcy/1p_1p/main.cc
@@ -0,0 +1,353 @@
+// -*- 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 Test for the equal dimension boundary coupling model
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+#include <fstream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+#include <dumux/io/grid/subgridgridcreator.hh>
+
+#include <dumux/porousmediumflow/1p/model.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/fluidsystems/1pliquid.hh>
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+#include <dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh>
+
+#include <dumux/multidomain/traits.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/newtonsolver.hh>
+#include <dumux/multidomain/boundary/darcydarcy/couplingmanager.hh>
+
+#include "problem.hh"
+
+#ifndef DOMAINSPLIT
+#define DOMAINSPLIT 0
+#endif
+
+namespace Dumux {
+namespace Properties {
+
+NEW_TYPE_TAG(OnePSubTypeTag, INHERITS_FROM(CCTpfaModel, OneP));
+// differentiate between the two subproblems
+NEW_TYPE_TAG(OnePSubTypeTag0, INHERITS_FROM(OnePSubTypeTag));
+NEW_TYPE_TAG(OnePSubTypeTag1, INHERITS_FROM(OnePSubTypeTag));
+
+// the coupling manager
+SET_TYPE_PROP(OnePSubTypeTag, CouplingManager,
+ DarcyDarcyBoundaryCouplingManager<MultiDomainTraits<TTAG(OnePSubTypeTag0), TTAG(OnePSubTypeTag1)>>);
+
+// Set the grid type
+#if DOMAINSPLIT==1
+SET_PROP(OnePSubTypeTag, Grid)
+{
+ using FullDomainGrid = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double, 2>>;
+ using type = Dune::SubGrid<FullDomainGrid::dimension, FullDomainGrid>;
+};
+#elif DOMAINSPLIT==0
+SET_PROP(OnePSubTypeTag, Grid)
+{
+ using FullDomainGrid = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double, 2>>;
+ using type = FullDomainGrid;
+};
+#endif
+
+// set the spatial params
+SET_TYPE_PROP(OnePSubTypeTag, SpatialParams, OnePTestSpatialParams<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
+ typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+// the fluid system
+SET_PROP(OnePSubTypeTag, FluidSystem)
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using type = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
+};
+
+// differentiate between the two subproblems
+SET_TYPE_PROP(OnePSubTypeTag0, Problem, OnePTestProblem<TypeTag, 0>);
+SET_TYPE_PROP(OnePSubTypeTag1, Problem, OnePTestProblem<TypeTag, 1>);
+
+} // end namespace Properties
+} // end namespace Dumux
+
+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 TypeTag = TTAG(OnePSubTypeTag);
+ using SubTypeTag0 = TTAG(OnePSubTypeTag0);
+ using SubTypeTag1 = TTAG(OnePSubTypeTag1);
+
+ // create the full grid that we are gonna split for the output
+ using FullDomainGrid = typename GET_PROP(TypeTag, Grid)::FullDomainGrid;
+ GridManager<FullDomainGrid> gridManager;
+ gridManager.init();
+
+ // get the lens
+ using GlobalPosition = typename FullDomainGrid::template Codim<0>::Geometry::GlobalCoordinate;
+ const auto lensLowerLeft = getParam<GlobalPosition>("SpatialParams.LensLowerLeft");
+ const auto lensUpperRight = getParam<GlobalPosition>("SpatialParams.LensUpperRight");
+
+///////////////////////////////////////////////////////////////////////////////////////////
+// split the domains by creating two separate grids for lens and the rest (using sub-grid)
+///////////////////////////////////////////////////////////////////////////////////////////
+#if DOMAINSPLIT==1
+
+ // create subgrids coinciding with the lens
+ auto elementSelector0 = [&lensLowerLeft, &lensUpperRight](const auto& element)
+ { return LensSpatialParams::pointInLens(element.geometry().center(), lensLowerLeft, lensUpperRight); };
+ auto elementSelector1 = [&lensLowerLeft, &lensUpperRight](const auto& element)
+ { return !LensSpatialParams::pointInLens(element.geometry().center(), lensLowerLeft, lensUpperRight); };
+
+ auto subGrid0 = SubgridGridCreator<FullDomainGrid>::makeGrid(gridManager.grid(), elementSelector0);
+ auto subGrid1 = SubgridGridCreator<FullDomainGrid>::makeGrid(gridManager.grid(), elementSelector1);
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid views
+ const auto& gridView0 = subGrid0->leafGridView();
+ const auto& gridView1 = subGrid1->leafGridView();
+
+//////////////////////////////////////////////////////////////////////////////////
+// split the domains by creating two separate grids for the lower and upper half
+//////////////////////////////////////////////////////////////////////////////////
+#elif DOMAINSPLIT==0
+
+ // create an upper half and lower half grid
+ GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager0, gridManager1;
+ gridManager0.init("1");
+ gridManager1.init("2");
+
+ // we compute on the leaf grid views
+ const auto& gridView0 = gridManager0.grid().leafGridView();
+ const auto& gridView1 = gridManager1.grid().leafGridView();
+
+#endif
+
+ ////////////////////////////////////////////////
+ // run the multidomain simulation on two grids
+ ////////////////////////////////////////////////
+
+ // create the finite volume grid geometries
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry0 = std::make_shared<FVGridGeometry>(gridView0);
+ auto fvGridGeometry1 = std::make_shared<FVGridGeometry>(gridView1);
+ fvGridGeometry0->update();
+ fvGridGeometry1->update();
+
+ // the mixed dimension type traits
+ using Traits = MultiDomainTraits<SubTypeTag0, SubTypeTag1>;
+ constexpr auto domain0Idx = Traits::template DomainIdx<0>();
+ constexpr auto domain1Idx = Traits::template DomainIdx<1>();
+
+ // the coupling manager
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ auto couplingManager = std::make_shared<CouplingManager>();
+
+ // the problem (initial and boundary conditions)
+ using Problem0 = typename GET_PROP_TYPE(SubTypeTag0, Problem);
+ auto problem0 = std::make_shared<Problem0>(fvGridGeometry0, couplingManager, "1");
+ problem0->spatialParams().setLens(lensLowerLeft, lensUpperRight);
+ using Problem1 = typename GET_PROP_TYPE(SubTypeTag1, Problem);
+ auto problem1 = std::make_shared<Problem1>(fvGridGeometry1, couplingManager, "2");
+ problem1->spatialParams().setLens(lensLowerLeft, lensUpperRight);
+
+ // the solution vector
+ Traits::SolutionVector sol;
+ sol[domain0Idx].resize(fvGridGeometry0->numDofs());
+ sol[domain1Idx].resize(fvGridGeometry1->numDofs());
+ problem0->applyInitialSolution(sol[domain0Idx]);
+ problem1->applyInitialSolution(sol[domain1Idx]);
+ auto oldSol = sol;
+
+ // compute the coupling map
+ couplingManager->init(problem0, problem1, sol);
+
+ // the grid variables
+ using GridVariables0 = typename GET_PROP_TYPE(SubTypeTag0, GridVariables);
+ auto gridVariables0 = std::make_shared<GridVariables0>(problem0, fvGridGeometry0);
+ using GridVariables1 = typename GET_PROP_TYPE(SubTypeTag1, GridVariables);
+ auto gridVariables1 = std::make_shared<GridVariables1>(problem1, fvGridGeometry1);
+ gridVariables0->init(sol[domain0Idx], oldSol[domain0Idx]);
+ gridVariables1->init(sol[domain1Idx], oldSol[domain1Idx]);
+
+ // get some time loop parameters
+ using Scalar = Traits::Scalar;
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // intialize the vtk output module
+ VtkOutputModule<SubTypeTag0> vtkWriter0(*problem0, *fvGridGeometry0, *gridVariables0, sol[domain0Idx], problem0->name());
+ GET_PROP_TYPE(SubTypeTag0, VtkOutputFields)::init(vtkWriter0);
+ vtkWriter0.write(0.0);
+
+ VtkOutputModule<SubTypeTag1> vtkWriter1(*problem1, *fvGridGeometry1, *gridVariables1, sol[domain1Idx], problem1->name());
+ GET_PROP_TYPE(SubTypeTag1, VtkOutputFields)::init(vtkWriter1);
+ vtkWriter1.write(0.0);
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(std::make_tuple(problem0, problem1),
+ std::make_tuple(fvGridGeometry0, fvGridGeometry1),
+ std::make_tuple(gridVariables0, gridVariables1),
+ couplingManager, timeLoop);
+
+ // the linear solver
+ using LinearSolver = ILU0BiCGSTABBackend;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the non-linear solver
+ using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+
+ // time loop
+ timeLoop->start();
+ while (!timeLoop->finished())
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(oldSol);
+
+ // solve the non-linear system with time step control
+ nonLinearSolver.solve(sol, *timeLoop);
+
+ // make the new solution the old solution
+ oldSol = sol;
+ gridVariables0->advanceTimeStep();
+ gridVariables1->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // write vtk output
+ vtkWriter0.write(timeLoop->time());
+ vtkWriter1.write(timeLoop->time());
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by newton controller
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+ }
+
+ timeLoop->finalize(mpiHelper.getCollectiveCommunication());
+
+ //////////////////////////////////////////////////////////////////////////
+ // write out a combined vtu for vtu comparison in the testing framework
+ //////////////////////////////////////////////////////////////////////////
+
+ const auto& gridView = gridManager.grid().leafGridView();
+ CCTpfaFVGridGeometry<typename FullDomainGrid::LeafGridView> fvGridGeometry(gridView);
+ const auto& bBoxTree = fvGridGeometry.boundingBoxTree();
+ // copy data from the subdomains to full domain data vectors
+ std::vector<int> processRank(gridView.size(0), 0); // sequential simulation
+ std::vector<int> pressure(gridView.size(0));
+
+ for (const auto& element : elements(gridView0))
+ {
+ const auto eIdx = fvGridGeometry0->elementMapper().index(element);
+ const auto eIdxHost = intersectingEntities(element.geometry().center(), bBoxTree)[0];
+ pressure[eIdxHost] = sol[domain0Idx][eIdx][0];
+ }
+
+ for (const auto& element : elements(gridView1))
+ {
+ const auto eIdx = fvGridGeometry1->elementMapper().index(element);
+ const auto eIdxHost = intersectingEntities(element.geometry().center(), bBoxTree)[0];
+ pressure[eIdxHost] = sol[domain1Idx][eIdx][0];
+ }
+
+ Dune::VTKWriter<typename FullDomainGrid::LeafGridView> vtkWriter(gridView);
+ vtkWriter.addCellData(processRank, "process rank");
+ vtkWriter.addCellData(pressure, "pressure");
+ const auto filename = getParam<std::string>("Problem.Name") + "_combined";
+ vtkWriter.write(filename);
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+}
+
+// error handler
+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/test/multidomain/boundary/darcydarcy/1p_1p/problem.hh b/test/multidomain/boundary/darcydarcy/1p_1p/problem.hh
new file mode 100644
index 0000000000..35f82fa8d2
--- /dev/null
+++ b/test/multidomain/boundary/darcydarcy/1p_1p/problem.hh
@@ -0,0 +1,163 @@
+// -*- 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 properties for the incompressible test
+ */
+#ifndef DUMUX_ONEP_SUB_TEST_PROBLEM_HH
+#define DUMUX_ONEP_SUB_TEST_PROBLEM_HH
+
+#include <dumux/porousmediumflow/problem.hh>
+#include "spatialparams.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup OnePTests
+ * \brief Multidomain test problem for the incompressible one-phase model
+ * \todo doc me!
+ */
+template<class TypeTag, std::size_t tag>
+class OnePTestProblem
+: public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
+ using GridView = typename FVGridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ static constexpr int dimWorld = GridView::dimensionworld;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ static constexpr auto domainIdx = Dune::index_constant<tag>{};
+
+public:
+ OnePTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager,
+ const std::string& paramGroup = "")
+ : ParentType(fvGridGeometry, paramGroup)
+ , couplingManager_(couplingManager)
+ {}
+
+ /*!
+ * \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.ipGlobal();
+
+ if (globalPos[dimWorld-1] < this->fvGridGeometry().bBoxMin()[dimWorld-1] + eps_
+ || globalPos[dimWorld-1] > this->fvGridGeometry().bBoxMax()[dimWorld-1] - eps_)
+ values.setAllDirichlet();
+ else
+ values.setAllNeumann();
+
+ if (couplingManager_->isCoupled(domainIdx, scvf))
+ values.setAllCouplingNeumann();
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * This is the method for the case where the Neumann condition is
+ * potentially solution dependent
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry
+ * \param elemVolVars All volume variables for the element
+ * \param scvf The sub control volume face
+ *
+ * Negative values mean influx.
+ * E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
+ */
+ template<class ElementVolumeVariables>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+ const auto bcTypes = boundaryTypes(element, scvf);
+ if (bcTypes.hasCouplingNeumann())
+ values[Indices::conti0EqIdx] = couplingManager_->advectiveFluxCoupling(domainIdx, element, fvGeometry, elemVolVars, scvf);
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * control volume.
+ *
+ * \param values The dirichlet values for the primary variables
+ * \param globalPos The center of the finite volume which ought to be set.
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+ values[0] = 1.0e+5*(2.0 - globalPos[dimWorld-1]);
+ return values;
+ }
+
+ /*!
+ * \brief Returns the temperature \f$\mathrm{[K]}\f$ for an isothermal problem.
+ *
+ * This is not specific to the discretization. By default it just
+ * throws an exception so it must be overloaded by the problem if
+ * no energy equation is used.
+ */
+ Scalar temperature() const
+ {
+ return 283.15; // 10°C
+ }
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ { return dirichletAtPos(globalPos); }
+
+
+private:
+ std::shared_ptr<CouplingManager> couplingManager_;
+ static constexpr Scalar eps_ = 1e-7;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/multidomain/boundary/darcydarcy/1p_1p/spatialparams.hh b/test/multidomain/boundary/darcydarcy/1p_1p/spatialparams.hh
new file mode 100644
index 0000000000..9700bad520
--- /dev/null
+++ b/test/multidomain/boundary/darcydarcy/1p_1p/spatialparams.hh
@@ -0,0 +1,114 @@
+// -*- 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 params the incompressible test
+ */
+#ifndef DUMUX_INCOMPRESSIBLE_ONEP_TEST_SPATIAL_PARAMS_HH
+#define DUMUX_INCOMPRESSIBLE_ONEP_TEST_SPATIAL_PARAMS_HH
+
+#include <limits>
+#include <dumux/material/spatialparams/fv1p.hh>
+
+namespace Dumux {
+namespace LensSpatialParams {
+/*!
+ * \brief If a point is in a lens with a given bounding box
+ *
+ * \param globalPos the position of the point
+ * \param lensLowerLeft the lower left corner of the lens
+ * \param lensUpperRight the upper right corner of the lens
+ */
+template<class GlobalPosition>
+bool pointInLens(const GlobalPosition& globalPos,
+ const GlobalPosition& lensLowerLeft,
+ const GlobalPosition& lensUpperRight)
+{
+ const auto eps = 1e-8*(lensUpperRight - lensLowerLeft).two_norm();
+ for (int i = 0; i < GlobalPosition::size(); ++i)
+ if (globalPos[i] < lensLowerLeft[i] + eps || globalPos[i] > lensUpperRight[i] - eps)
+ return false;
+
+ return true;
+}
+} // end namespace LensSpatialParams
+
+/*!
+ * \ingroup OnePTests
+ * \brief The spatial parameters class for the test problem using the
+ * incompressible 1p model
+ */
+template<class FVGridGeometry, class Scalar>
+class OnePTestSpatialParams
+: public FVSpatialParamsOneP<FVGridGeometry, Scalar, OnePTestSpatialParams<FVGridGeometry, Scalar>>
+{
+ using ThisType = OnePTestSpatialParams<FVGridGeometry, Scalar>;
+ using ParentType = FVSpatialParamsOneP<FVGridGeometry, Scalar, ThisType>;
+ using GridView = typename FVGridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+
+ static constexpr int dimWorld = GridView::dimensionworld;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+ using PermeabilityType = Scalar;
+ OnePTestSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ , lensLowerLeft_(std::numeric_limits<Scalar>::max())
+ , lensUpperRight_(std::numeric_limits<Scalar>::lowest())
+ {
+ permeability_ = getParam<Scalar>("SpatialParams.Permeability");
+ permeabilityLens_ = getParam<Scalar>("SpatialParams.PermeabilityLens");
+ }
+
+ /*!
+ * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
+ * \return the intrinsic permeability
+ */
+ PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
+ { return isInLens_(globalPos) ? permeabilityLens_ : permeability_; }
+
+ /*!
+ * \brief Define the porosity \f$\mathrm{[-]}\f$.
+ *
+ * \param globalPos The global position
+ */
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
+ { return 0.4; }
+
+ /*!
+ * \brief Optionally set a lens
+ */
+ void setLens(const GlobalPosition& lowerLeft, const GlobalPosition& upperRight)
+ { lensLowerLeft_ = lowerLeft; lensUpperRight_ = upperRight; }
+
+private:
+ bool isInLens_(const GlobalPosition &globalPos) const
+ { return LensSpatialParams::pointInLens(globalPos, lensLowerLeft_, lensUpperRight_); }
+
+ GlobalPosition lensLowerLeft_, lensUpperRight_;
+ Scalar permeability_, permeabilityLens_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/multidomain/boundary/darcydarcy/1p_1p/test_boundary_equaldim_1p_1p.input b/test/multidomain/boundary/darcydarcy/1p_1p/test_boundary_equaldim_1p_1p.input
new file mode 100644
index 0000000000..4db37b49c2
--- /dev/null
+++ b/test/multidomain/boundary/darcydarcy/1p_1p/test_boundary_equaldim_1p_1p.input
@@ -0,0 +1,34 @@
+[TimeLoop]
+TEnd = 1.0
+DtInitial = 1.0
+
+[Problem]
+EnableGravity = true
+
+[1.Problem]
+Name = 1p_0
+
+[2.Problem]
+Name = 1p_1
+
+[Grid]
+LowerLeft = 0 0
+UpperRight = 1 1
+Cells = 10 10
+
+[1.Grid]
+LowerLeft = 0 0
+UpperRight = 1 0.5
+Cells = 10 5
+
+[2.Grid]
+LowerLeft = 0 0.5
+UpperRight = 1 1
+Cells = 10 5
+
+[SpatialParams]
+LensLowerLeft = 0.2 0.2
+LensUpperRight = 0.8 0.8
+
+Permeability = 1e-10 # [m^2]
+PermeabilityLens = 1e-12 # [m^2]
diff --git a/test/multidomain/boundary/darcydarcy/CMakeLists.txt b/test/multidomain/boundary/darcydarcy/CMakeLists.txt
new file mode 100644
index 0000000000..f62f49df14
--- /dev/null
+++ b/test/multidomain/boundary/darcydarcy/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory(1p_1p)
--
GitLab