From 9c3ab5e1a43bff526fd6590c71f3b3150d29d4dc Mon Sep 17 00:00:00 2001
From: Klaus Mosthaf <klmos@env.dtu.dk>
Date: Wed, 29 Jan 2014 10:27:05 +0000
Subject: [PATCH] added 2cnistokes2p2cni test and adapted makefiles and
configure
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@12384 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
configure.ac | 1 +
.../2cnistokes2p2cnilocaloperator.hh | 4 +-
.../2cnistokes2p2cniproblem.hh | 758 ++++++++++++++++++
.../2cnistokes2p2cnispatialparams.hh | 371 +++++++++
.../2cnistokes2p2cni/2p2cnisubproblem.hh | 675 ++++++++++++++++
test/multidomain/2cnistokes2p2cni/Makefile.am | 5 +
.../grids/interfacedomain.dgf | 15 +
.../grids/test_2cnistokes2p2cni.dgf | 17 +
.../2cnistokes2p2cni/stokes2cnisubproblem.hh | 592 ++++++++++++++
.../2cnistokes2p2cni/test_2cnistokes2p2cni.cc | 279 +++++++
.../test_2cnistokes2p2cni.input | 193 +++++
test/multidomain/Makefile.am | 3 +-
12 files changed, 2910 insertions(+), 3 deletions(-)
create mode 100644 test/multidomain/2cnistokes2p2cni/2cnistokes2p2cniproblem.hh
create mode 100644 test/multidomain/2cnistokes2p2cni/2cnistokes2p2cnispatialparams.hh
create mode 100644 test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh
create mode 100644 test/multidomain/2cnistokes2p2cni/Makefile.am
create mode 100644 test/multidomain/2cnistokes2p2cni/grids/interfacedomain.dgf
create mode 100644 test/multidomain/2cnistokes2p2cni/grids/test_2cnistokes2p2cni.dgf
create mode 100644 test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh
create mode 100644 test/multidomain/2cnistokes2p2cni/test_2cnistokes2p2cni.cc
create mode 100644 test/multidomain/2cnistokes2p2cni/test_2cnistokes2p2cni.input
diff --git a/configure.ac b/configure.ac
index e311465a31..10c58d3271 100644
--- a/configure.ac
+++ b/configure.ac
@@ -131,6 +131,7 @@ AC_CONFIG_FILES([dumux.pc
test/material/tabulation/Makefile
test/multidomain/Makefile
test/multidomain/2cstokes2p2c/Makefile
+ test/multidomain/2cnistokes2p2cni/Makefile
test/references/Makefile
tutorial/Makefile
])
diff --git a/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh b/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh
index bdc7153fc4..7672bed006 100644
--- a/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh
+++ b/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh
@@ -26,8 +26,8 @@
#include <dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh>
-#include <dumux/freeflow/stokesncni/stokesncnimodel.hh>
-#include <dumux/implicit/2p2cnicoupling/2p2cnicouplingmodel.hh>
+//#include <dumux/freeflow/stokesncni/stokesncnimodel.hh>
+//#include <dumux/implicit/2p2cnicoupling/2p2cnicouplingmodel.hh>
namespace Dumux {
diff --git a/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cniproblem.hh b/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cniproblem.hh
new file mode 100644
index 0000000000..9c3ca6c470
--- /dev/null
+++ b/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cniproblem.hh
@@ -0,0 +1,758 @@
+// -*- 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 docme
+ * @brief docme
+ */
+#ifndef DUMUX_2CNISTOKES2P2CNIPROBLEM_HH
+#define DUMUX_2CNISTOKES2P2CNIPROBLEM_HH
+
+#include <dune/grid/common/gridinfo.hh>
+#include <dune/grid/multidomaingrid.hh>
+#include <dune/grid/io/file/dgfparser.hh>
+
+#include <dumux/multidomain/common/multidomainproblem.hh>
+#include <dumux/multidomain/2cstokes2p2c/2cstokes2p2cnewtoncontroller.hh>
+#include <dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh>
+//#include <dumux/linear/amgbackend.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/linear/pardisobackend.hh>
+
+#include "2cnistokes2p2cnispatialparams.hh"
+#include <dumux/material/fluidsystems/h2oairfluidsystem.hh>
+
+#include "stokes2cnisubproblem.hh"
+#include "2p2cnisubproblem.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class TwoCNIStokesTwoPTwoCNIProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(TwoCNIStokesTwoPTwoCNIProblem, INHERITS_FROM(MultiDomain));
+
+// Set the grid type
+SET_PROP(TwoCNIStokesTwoPTwoCNIProblem, Grid)
+{
+ public:
+#ifdef HAVE_UG
+ typedef typename Dune::UGGrid<2> type;
+#elif HAVE_ALUGRID
+ typedef typename Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming> type;
+#endif
+};
+
+// Specify the multidomain gridview
+SET_PROP(TwoCNIStokesTwoPTwoCNIProblem, GridView)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, MultiDomainGrid) MDGrid;
+ public:
+ typedef typename MDGrid::LeafGridView type;
+};
+
+// Set the global problem
+SET_TYPE_PROP(TwoCNIStokesTwoPTwoCNIProblem, Problem, TwoCNIStokesTwoPTwoCNIProblem<TypeTag>);
+
+// Set the two sub-problems of the global problem
+SET_TYPE_PROP(TwoCNIStokesTwoPTwoCNIProblem, SubDomain1TypeTag, TTAG(Stokes2cniSubProblem));
+SET_TYPE_PROP(TwoCNIStokesTwoPTwoCNIProblem, SubDomain2TypeTag, TTAG(TwoPTwoCNISubProblem));
+
+// Set the global problem in the context of the two sub-problems
+SET_TYPE_PROP(Stokes2cniSubProblem, MultiDomainTypeTag, TTAG(TwoCNIStokesTwoPTwoCNIProblem));
+SET_TYPE_PROP(TwoPTwoCNISubProblem, MultiDomainTypeTag, TTAG(TwoCNIStokesTwoPTwoCNIProblem));
+
+// Set the other sub-problem for each of the two sub-problems
+SET_TYPE_PROP(Stokes2cniSubProblem, OtherSubDomainTypeTag, TTAG(TwoPTwoCNISubProblem));
+SET_TYPE_PROP(TwoPTwoCNISubProblem, OtherSubDomainTypeTag, TTAG(Stokes2cniSubProblem));
+
+// Set the same spatial parameters for both sub-problems
+SET_PROP(Stokes2cniSubProblem, SpatialParams)
+{ typedef Dumux::TwoCNIStokesTwoPTwoCNISpatialParams<TypeTag> type; };
+SET_PROP(TwoPTwoCNISubProblem, SpatialParams)
+{ typedef Dumux::TwoCNIStokesTwoPTwoCNISpatialParams<TypeTag> type; };
+
+// Set the fluid system
+SET_PROP(TwoCNIStokesTwoPTwoCNIProblem, FluidSystem)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef Dumux::FluidSystems::H2OAir<Scalar,
+ Dumux::H2O<Scalar>,
+ /*useComplexrelations=*/true> type;
+};
+
+SET_TYPE_PROP(TwoCNIStokesTwoPTwoCNIProblem, JacobianAssembler, Dumux::MultiDomainAssembler<TypeTag>);
+
+// Set the local coupling operator
+SET_TYPE_PROP(TwoCNIStokesTwoPTwoCNIProblem, MultiDomainCouplingLocalOperator,
+ Dumux::TwoCNIStokesTwoPTwoCNILocalOperator<TypeTag>);
+
+SET_PROP(TwoCNIStokesTwoPTwoCNIProblem, SolutionVector)
+{
+ private:
+ typedef typename GET_PROP_TYPE(TypeTag, MultiDomainGridOperator) MDGridOperator;
+ public:
+ typedef typename MDGridOperator::Traits::Domain type;
+};
+
+// the newton controller from the TwoCStokesTwoPTwoC model is used
+SET_TYPE_PROP(TwoCNIStokesTwoPTwoCNIProblem, NewtonController, TwoCStokesTwoPTwoCNewtonController<TypeTag>);
+
+//SET_TYPE_PROP(TwoCNIStokesTwoPTwoCNIProblem, LinearSolver, SSORRestartedGMResBackend<TypeTag>);
+// if you do not have PARDISO, the SuperLU solver is used:
+#ifdef HAVE_PARDISO
+SET_TYPE_PROP(TwoCNIStokesTwoPTwoCNIProblem, LinearSolver, PardisoBackend<TypeTag>);
+#else
+SET_TYPE_PROP(TwoCNIStokesTwoPTwoCNIProblem, LinearSolver, SuperLUBackend<TypeTag>);
+#endif // HAVE_PARDISO
+
+
+// set this to one here (must fit to the structure of the coupled matrix which has block length 1)
+SET_INT_PROP(TwoCNIStokesTwoPTwoCNIProblem, NumEq, 1);
+
+// Write the solutions of individual newton iterations?
+SET_BOOL_PROP(TwoCNIStokesTwoPTwoCNIProblem, NewtonWriteConvergence, false);
+}
+
+/*!
+ * \brief The problem class for the coupling of a non-isothermal two-component Stokes (stokes2cni)
+ * and a non-isothermal two-phase two-component Darcy model (2p2cni).
+ *
+ * The problem class for the coupling of a non-isothermal two-component Stokes (stokes2cni)
+ * and a non-isothermal two-phase two-component Darcy model (2p2cni).
+ * It uses the 2p2cniCoupling model and the Stokes2cnicoupling model and provides
+ * the problem specifications for common parameters of the two submodels.
+ * The initial and boundary conditions of the submodels are specified in the two subproblems,
+ * 2p2cnisubproblem.hh and stokes2cnisubproblem.hh, which are accessible via the coupled problem.
+ */
+template <class TypeTag = TTAG(TwoCNIStokesTwoPTwoCNIProblem) >
+class TwoCNIStokesTwoPTwoCNIProblem : public MultiDomainProblem<TypeTag>
+{
+ template<int dim>
+ struct VertexLayout
+ {
+ bool contains(Dune::GeometryType geomtype)
+ { return geomtype.dim() == 0; }
+ };
+
+ typedef TwoCNIStokesTwoPTwoCNIProblem<TypeTag> ThisType;
+ typedef MultiDomainProblem<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+ typedef typename GET_PROP_TYPE(TypeTag, SubDomain1TypeTag) Stokes2cniTypeTag;
+ typedef typename GET_PROP_TYPE(TypeTag, SubDomain2TypeTag) TwoPTwoCNITypeTag;
+
+ typedef typename GET_PROP_TYPE(Stokes2cniTypeTag, Problem) Stokes2cniSubProblem;
+ typedef typename GET_PROP_TYPE(TwoPTwoCNITypeTag, Problem) TwoPTwoCNISubProblem;
+
+ typedef typename GET_PROP_TYPE(Stokes2cniTypeTag, GridView) Stokes2cniGridView;
+ typedef typename GET_PROP_TYPE(TwoPTwoCNITypeTag, GridView) TwoPTwoCNIGridView;
+
+ typedef typename GET_PROP_TYPE(Stokes2cniTypeTag, PrimaryVariables) Stokes2cniPrimaryVariables;
+ typedef typename GET_PROP_TYPE(TwoPTwoCNITypeTag, PrimaryVariables) TwoPTwoCNIPrimaryVariables;
+
+ typedef typename GET_PROP_TYPE(Stokes2cniTypeTag, ElementSolutionVector) ElementSolutionVector1;
+ typedef typename GET_PROP_TYPE(TwoPTwoCNITypeTag, ElementSolutionVector) ElementSolutionVector2;
+
+ typedef typename GET_PROP_TYPE(Stokes2cniTypeTag, ElementVolumeVariables) ElementVolumeVariables1;
+ typedef typename GET_PROP_TYPE(TwoPTwoCNITypeTag, ElementVolumeVariables) ElementVolumeVariables2;
+
+ typedef typename GET_PROP_TYPE(Stokes2cniTypeTag, FluxVariables) BoundaryVariables1;
+
+ typedef typename GET_PROP_TYPE(Stokes2cniTypeTag, FVElementGeometry) FVElementGeometry1;
+ typedef typename GET_PROP_TYPE(TwoPTwoCNITypeTag, FVElementGeometry) FVElementGeometry2;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Grid) HostGrid;
+ typedef typename GET_PROP_TYPE(TypeTag, MultiDomainGrid) MDGrid;
+ typedef typename MDGrid::LeafGridView MDGridView;
+ typedef typename MDGrid::SubDomainGrid SDGrid;
+
+ typedef typename MDGrid::Traits::template Codim<0>::Entity MDElement;
+ typedef typename MDGrid::Traits::template Codim<0>::EntityPointer MDElementPointer;
+ typedef typename Stokes2cniGridView::template Codim<0>::Entity SDElement1;
+ typedef typename TwoPTwoCNIGridView::template Codim<0>::Entity SDElement2;
+ typedef typename SDGrid::Traits::template Codim<0>::EntityPointer SDElementPointer;
+
+ typedef typename GET_PROP_TYPE(Stokes2cniTypeTag, Indices) Stokes2cniIndices;
+ typedef typename GET_PROP_TYPE(TwoPTwoCNITypeTag, Indices) TwoPTwoCNIIndices;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+ enum { dim = Stokes2cniGridView::dimension };
+ enum { // indices in the Stokes domain
+ momentumXIdx1 = Stokes2cniIndices::momentumXIdx, //!< Index of the x-component of the momentum balance
+ momentumYIdx1 = Stokes2cniIndices::momentumYIdx, //!< Index of the y-component of the momentum balance
+ momentumZIdx1 = Stokes2cniIndices::momentumZIdx, //!< Index of the z-component of the momentum balance
+ massBalanceIdx1 = Stokes2cniIndices::massBalanceIdx, //!< Index of the mass balance
+ transportEqIdx1 = Stokes2cniIndices::transportEqIdx, //!< Index of the transport equation
+ energyEqIdx1 = Stokes2cniIndices::energyEqIdx //!< Index of the transport equation
+ };
+ enum { // indices of the PVs in the Darcy domain
+ massBalanceIdx2 = TwoPTwoCNIIndices::pressureIdx,
+ switchIdx2 = TwoPTwoCNIIndices::switchIdx,
+ temperatureIdx2 = TwoPTwoCNIIndices::temperatureIdx
+ };
+ enum { // indices of the balance equations
+ contiTotalMassIdx2 = TwoPTwoCNIIndices::contiNEqIdx,
+ contiWEqIdx2 = TwoPTwoCNIIndices::contiWEqIdx,
+ energyEqIdx2 = TwoPTwoCNIIndices::energyEqIdx
+ };
+ enum { transportCompIdx1 = Stokes2cniIndices::transportCompIdx };
+ enum {
+ wCompIdx2 = TwoPTwoCNIIndices::wCompIdx,
+ nCompIdx2 = TwoPTwoCNIIndices::nCompIdx
+ };
+ enum { phaseIdx = GET_PROP_VALUE(Stokes2cniTypeTag, PhaseIdx) };
+ enum {
+ numEq1 = GET_PROP_VALUE(Stokes2cniTypeTag, NumEq),
+ numEq2 = GET_PROP_VALUE(TwoPTwoCNITypeTag, NumEq)
+ };
+
+ typedef Dune::FieldVector<Scalar, dim> GlobalPosition;
+ typedef Dune::FieldVector<Scalar, dim> FieldVector;
+
+ typedef typename MDGrid::template Codim<0>::LeafIterator ElementIterator;
+ typedef typename MDGrid::LeafSubDomainInterfaceIterator SDInterfaceIterator;
+
+ /*!
+ * \brief docme
+ *
+ * \param hostGrid docme
+ * \param timeManager The TimeManager which is used by the simulation
+ *
+ */
+public:
+ TwoCNIStokesTwoPTwoCNIProblem(MDGrid &mdGrid,
+ TimeManager &timeManager)
+ : ParentType(mdGrid, timeManager)
+ {
+ try
+ {
+ // define location of the interface
+ interface_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePos);
+ noDarcyX_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, NoDarcyX);
+ episodeLength_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, EpisodeLength);
+ initializationTime_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, InitTime);
+ dtInit_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, DtInitial);
+
+ // define output options
+ freqRestart_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqRestart);
+ freqOutput_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqOutput);
+ freqFluxOutput_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqFluxOutput);
+ freqVaporFluxOutput_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqVaporFluxOutput);
+ }
+ catch (Dumux::ParameterException &e) {
+ std::cerr << e << ". Abort!\n";
+ exit(1) ;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ exit(1);
+ }
+
+ stokes2cni_ = this->sdID1();
+ twoPtwoCNI_ = this->sdID2();
+
+ initializeGrid();
+
+ // initialize the tables of the fluid system
+ FluidSystem::init(/*tempMin=*/273.15, /*tempMax=*/373.15, /*numTemp=*/200,
+ /*pMin=*/1e3, /*pMax=*/2e5, /*numP=*/200);
+
+ if (initializationTime_ > 0.0)
+ this->timeManager().startNextEpisode(initializationTime_);
+ else
+ this->timeManager().startNextEpisode(episodeLength_);
+ }
+
+ ~TwoCNIStokesTwoPTwoCNIProblem()
+ {
+ fluxFile_.close();
+ };
+
+ //! \copydoc Dumux::CoupledProblem::init()
+ void init()
+ {
+ ParentType::init();
+
+ // plot the Pc-Sw curves, if requested
+ this->sdProblem2().spatialParams().plotMaterialLaw();
+
+ std::cout << "Writing flux data at interface\n";
+ if (this->timeManager().time() == 0)
+ {
+ fluxFile_.open("fluxes.out");
+ fluxFile_ << "time; flux1; advFlux1; diffFlux1; totalFlux1; energyFlux1; "
+ "flux2; wPhaseFlux2; nPhaseFlux2; energyFlux2\n";
+ counter_ = 1;
+ }
+ else
+ fluxFile_.open("fluxes.out", std::ios_base::app);
+ }
+
+ /*!
+ * \brief docme
+ */
+ void initializeGrid()
+ {
+ MDGrid& mdGrid = this->mdGrid();
+ mdGrid.startSubDomainMarking();
+
+ // subdivide grid in two subdomains
+ ElementIterator eendit = mdGrid.template leafend<0>();
+ for (ElementIterator elementIt = mdGrid.template leafbegin<0>();
+ elementIt != eendit; ++elementIt)
+ {
+ // this is required for parallelization
+ // checks if element is within a partition
+ if (elementIt->partitionType() != Dune::InteriorEntity)
+ continue;
+
+ GlobalPosition globalPos = elementIt->geometry().center();
+
+ if (globalPos[1] > interface_)
+ mdGrid.addToSubDomain(stokes2cni_,*elementIt);
+ else
+ if(globalPos[0] > noDarcyX_)
+ mdGrid.addToSubDomain(twoPtwoCNI_,*elementIt);
+ }
+ mdGrid.preUpdateSubDomains();
+ mdGrid.updateSubDomains();
+ mdGrid.postUpdateSubDomains();
+
+ gridinfo(this->sdGrid1());
+ gridinfo(this->sdGrid2());
+ }
+
+ //! \copydoc Dumux::CoupledProblem::postTimeStep()
+ void postTimeStep()
+ {
+ // call the postTimeStep function of the subproblems
+ this->sdProblem1().postTimeStep();
+ this->sdProblem2().postTimeStep();
+
+ if (shouldWriteFluxFile() || shouldWriteVaporFlux())
+ {
+ counter_ = this->sdProblem1().currentVTKFileNumber() + 1;
+
+ calculateFirstInterfaceFluxes();
+ calculateSecondInterfaceFluxes();
+ }
+ }
+
+ //! \copydoc Dumux::CoupledProblem::episodeEnd()
+ void episodeEnd()
+ {
+ this->timeManager().startNextEpisode(episodeLength_);
+ if (this->timeManager().time() <= initializationTime_ + dtInit_)
+ {
+ std::cout << "setting timeStepSize to " << dtInit_ << std::endl;
+ this->timeManager().setTimeStepSize(dtInit_);
+ }
+ }
+
+ /*
+ * \brief Calculates fluxes and coupling terms at the interface for the Stokes model.
+ * Flux output files are created and the summarized flux is written to a file.
+ */
+ void calculateFirstInterfaceFluxes()
+ {
+ const MDGrid& mdGrid = this->mdGrid();
+ ElementVolumeVariables1 elemVolVarsPrev1, elemVolVarsCur1;
+ Scalar sumVaporFluxes = 0.;
+ Scalar advectiveVaporFlux = 0.;
+ Scalar diffusiveVaporFlux = 0.;
+ Scalar sumEnergyFluxes = 0.;
+
+ // count number of elements to determine number of interface nodes
+ int numElements = 0;
+ const SDInterfaceIterator endIfIt = mdGrid.leafSubDomainInterfaceEnd(stokes2cni_, twoPtwoCNI_);
+ for (SDInterfaceIterator ifIt =
+ mdGrid.leafSubDomainInterfaceBegin(stokes2cni_, twoPtwoCNI_); ifIt != endIfIt; ++ifIt)
+ numElements++;
+
+ const int numInterfaceVertices = numElements + 1;
+ std::vector<InterfaceFluxes<numEq1> > outputVector(numInterfaceVertices); // vector for the output of the fluxes
+ FVElementGeometry1 fvGeometry1;
+ int interfaceVertIdx = -1;
+
+ // loop over the element faces on the interface
+ for (SDInterfaceIterator ifIt =
+ mdGrid.leafSubDomainInterfaceBegin(stokes2cni_, twoPtwoCNI_); ifIt != endIfIt; ++ifIt)
+ {
+ const int firstFaceIdx = ifIt->indexInFirstCell();
+ const MDElementPointer mdElementPointer1 = ifIt->firstCell(); // ATTENTION!!!
+ const MDElement& mdElement1 = *mdElementPointer1; // Entity pointer has to be copied before.
+ const SDElementPointer sdElementPointer1 = this->sdElementPointer1(mdElement1);
+ const SDElement1& sdElement1 = *sdElementPointer1;
+ fvGeometry1.update(this->sdGridView1(), sdElement1);
+
+ const Dune::GenericReferenceElement<typename MDGrid::ctype,dim>& referenceElement1 =
+ Dune::GenericReferenceElements<typename MDGrid::ctype,dim>::general(mdElement1.type());
+ const int numVerticesOfFace = referenceElement1.size(firstFaceIdx, 1, dim);
+
+ // evaluate residual of the sub model without boundary conditions (stabilization is removed)
+ // the element volume variables are updated here
+ this->localResidual1().evalNoBoundary(sdElement1, fvGeometry1,
+ elemVolVarsPrev1, elemVolVarsCur1);
+
+ for (int nodeInFace = 0; nodeInFace < numVerticesOfFace; nodeInFace++)
+ {
+ const int vertInElem1 = referenceElement1.subEntity(firstFaceIdx, 1, nodeInFace, dim);
+ const FieldVector& vertexGlobal = mdElement1.geometry().corner(vertInElem1);
+ const unsigned firstGlobalIdx = this->mdVertexMapper().map(stokes2cni_, mdElement1, vertInElem1, dim);
+ const ElementSolutionVector1& firstVertexDefect = this->localResidual1().residual();
+
+ // loop over all interface vertices to check if vertex id is already in stack
+ bool existing = false;
+ for (int interfaceVertex=0; interfaceVertex < numInterfaceVertices; ++interfaceVertex)
+ {
+ if (firstGlobalIdx == outputVector[interfaceVertex].globalIdx)
+ {
+ existing = true;
+ interfaceVertIdx = interfaceVertex;
+ break;
+ }
+ }
+
+ if (!existing)
+ interfaceVertIdx++;
+
+ if (shouldWriteFluxFile()) // compute only if required
+ {
+ outputVector[interfaceVertIdx].interfaceVertex = interfaceVertIdx;
+ outputVector[interfaceVertIdx].globalIdx = firstGlobalIdx;
+ outputVector[interfaceVertIdx].xCoord = vertexGlobal[0];
+ outputVector[interfaceVertIdx].yCoord = vertexGlobal[1];
+ outputVector[interfaceVertIdx].count += 1;
+ for (int eqIdx=0; eqIdx < numEq1; ++eqIdx)
+ outputVector[interfaceVertIdx].defect[eqIdx] +=
+ firstVertexDefect[vertInElem1][eqIdx];
+ }
+
+ // compute summarized fluxes for output
+ if (shouldWriteVaporFlux())
+ {
+ int boundaryFaceIdx =
+ fvGeometry1.boundaryFaceIndex(firstFaceIdx, nodeInFace);
+
+ const BoundaryVariables1 boundaryVars1(this->sdProblem1(),
+ sdElement1,
+ fvGeometry1,
+ boundaryFaceIdx,
+ elemVolVarsCur1,
+ /*onBoundary=*/true);
+
+ advectiveVaporFlux += computeAdvectiveVaporFluxes1(elemVolVarsCur1, boundaryVars1, vertInElem1);
+ diffusiveVaporFlux += computeDiffusiveVaporFluxes1(elemVolVarsCur1, boundaryVars1, vertInElem1);
+ sumVaporFluxes += firstVertexDefect[vertInElem1][transportEqIdx1];
+ sumEnergyFluxes += firstVertexDefect[vertInElem1][energyEqIdx1];
+ }
+ }
+ } // end loop over element faces on interface
+
+ if (shouldWriteFluxFile())
+ {
+ std::cout << "Writing flux file\n";
+ char outputname[20];
+ sprintf(outputname, "%s%05d%s","fluxes1_", counter_,".out");
+ std::ofstream outfile(outputname, std::ios_base::out);
+ outfile << "Xcoord1 "
+ << "totalFlux1 "
+ << "componentFlux1 "
+ << "heatFlux1 "
+ << std::endl;
+ for (int interfaceVertIdx=0; interfaceVertIdx < numInterfaceVertices; interfaceVertIdx++)
+ {
+ if (outputVector[interfaceVertIdx].count > 2)
+ std::cerr << "too often at one node!!";
+
+ if (outputVector[interfaceVertIdx].count==2)
+ outfile << outputVector[interfaceVertIdx].xCoord << " "
+ << outputVector[interfaceVertIdx].defect[massBalanceIdx1] << " " // total mass flux
+ << outputVector[interfaceVertIdx].defect[transportEqIdx1] << " " // total flux of component
+ << outputVector[interfaceVertIdx].defect[energyEqIdx1] << " " // total flux of heat
+ << std::endl;
+ }
+ outfile.close();
+ }
+ if (shouldWriteVaporFlux())
+ fluxFile_ << this->timeManager().time() + this->timeManager().timeStepSize() << "; "
+ << sumVaporFluxes << "; "
+ << advectiveVaporFlux << "; "
+ << diffusiveVaporFlux << "; "
+ << advectiveVaporFlux-diffusiveVaporFlux << "; "
+ << sumEnergyFluxes << "; ";
+ }
+
+ /*
+ * \brief Calculates fluxes and coupling terms at the interface for the Darcy model.
+ * Flux output files are created and the summarized flux is written to a file.
+ */
+ void calculateSecondInterfaceFluxes()
+ {
+ const MDGrid& mdGrid = this->mdGrid();
+ ElementVolumeVariables2 elemVolVarsPrev2, elemVolVarsCur2;
+
+ Scalar sumVaporFluxes = 0.;
+ Scalar sumEnergyFluxes = 0.;
+ Scalar sumWaterFluxInGasPhase = 0.;
+
+ // count number of elements to determine number of interface nodes
+ int numElements = 0;
+ const SDInterfaceIterator endIfIt = mdGrid.leafSubDomainInterfaceEnd(stokes2cni_, twoPtwoCNI_);
+ for (SDInterfaceIterator ifIt =
+ mdGrid.leafSubDomainInterfaceBegin(stokes2cni_, twoPtwoCNI_); ifIt != endIfIt; ++ifIt)
+ numElements++;
+
+ const int numInterfaceVertices = numElements + 1;
+ std::vector<InterfaceFluxes<numEq2> > outputVector(numInterfaceVertices); // vector for the output of the fluxes
+ FVElementGeometry2 fvGeometry2;
+ int interfaceVertIdx = -1;
+
+ // loop over the element faces on the interface
+ for (SDInterfaceIterator ifIt =
+ mdGrid.leafSubDomainInterfaceBegin(stokes2cni_, twoPtwoCNI_); ifIt != endIfIt; ++ifIt)
+ {
+ const int secondFaceIdx = ifIt->indexInSecondCell();
+ const MDElementPointer mdElementPointer2 = ifIt->secondCell(); // ATTENTION!!!
+ const MDElement& mdElement2 = *mdElementPointer2; // Entity pointer has to be copied before.
+ const SDElementPointer sdElementPointer2 = this->sdElementPointer2(mdElement2);
+ const SDElement2& sdElement2 = *sdElementPointer2;
+ fvGeometry2.update(this->sdGridView2(), sdElement2);
+
+ const Dune::GenericReferenceElement<typename MDGrid::ctype,dim>& referenceElement2 =
+ Dune::GenericReferenceElements<typename MDGrid::ctype,dim>::general(mdElement2.type());
+ const int numVerticesOfFace = referenceElement2.size(secondFaceIdx, 1, dim);
+
+ // evaluate residual of the sub model without boundary conditions
+ this->localResidual2().evalNoBoundary(sdElement2, fvGeometry2,
+ elemVolVarsPrev2, elemVolVarsCur2);
+ // evaluate the vapor fluxes within each phase
+ this->localResidual2().evalPhaseFluxes();
+
+ for (int nodeInFace = 0; nodeInFace < numVerticesOfFace; nodeInFace++)
+ {
+ const int vertInElem2 = referenceElement2.subEntity(secondFaceIdx, 1, nodeInFace, dim);
+ const FieldVector& vertexGlobal = mdElement2.geometry().corner(vertInElem2);
+ const unsigned secondGlobalIdx = this->mdVertexMapper().map(twoPtwoCNI_, mdElement2, vertInElem2, dim);
+ const ElementSolutionVector2& secondVertexDefect = this->localResidual2().residual();
+
+ bool existing = false;
+ // loop over all interface vertices to check if vertex id is already in stack
+ for (int interfaceVertex=0; interfaceVertex < numInterfaceVertices; ++interfaceVertex)
+ {
+ if (secondGlobalIdx == outputVector[interfaceVertex].globalIdx)
+ {
+ existing = true;
+ interfaceVertIdx = interfaceVertex;
+ break;
+ }
+ }
+
+ if (!existing)
+ interfaceVertIdx++;
+
+ if (shouldWriteFluxFile())
+ {
+ outputVector[interfaceVertIdx].interfaceVertex = interfaceVertIdx;
+ outputVector[interfaceVertIdx].globalIdx = secondGlobalIdx;
+ outputVector[interfaceVertIdx].xCoord = vertexGlobal[0];
+ outputVector[interfaceVertIdx].yCoord = vertexGlobal[1];
+ for (int eqIdx=0; eqIdx < numEq2; ++eqIdx)
+ outputVector[interfaceVertIdx].defect[eqIdx] += secondVertexDefect[vertInElem2][eqIdx];
+ outputVector[interfaceVertIdx].count += 1;
+ }
+ if (shouldWriteVaporFlux())
+ {
+ if (!existing) // add phase storage only once per vertex
+ sumWaterFluxInGasPhase +=
+ this->localResidual2().evalPhaseStorage(vertInElem2);
+
+ sumVaporFluxes += secondVertexDefect[vertInElem2][contiWEqIdx2];
+ sumWaterFluxInGasPhase +=
+ this->localResidual2().elementFluxes(vertInElem2);
+ sumEnergyFluxes += secondVertexDefect[vertInElem2][energyEqIdx2];
+
+ }
+ }
+ }
+
+ if (shouldWriteFluxFile())
+ {
+ char outputname[20];
+ sprintf(outputname, "%s%05d%s","fluxes2_", counter_,".out");
+ std::ofstream outfile(outputname, std::ios_base::out);
+ outfile << "Xcoord2 "
+ << "totalFlux2 "
+ << "componentFlux2 "
+ << "heatFlux2 "
+ << std::endl;
+
+ for (int interfaceVertIdx=0; interfaceVertIdx < numInterfaceVertices; interfaceVertIdx++)
+ {
+ if (outputVector[interfaceVertIdx].count > 2)
+ std::cerr << "too often at one node!!";
+
+ if (outputVector[interfaceVertIdx].count==2)
+ outfile << outputVector[interfaceVertIdx].xCoord << " "
+ << outputVector[interfaceVertIdx].defect[contiTotalMassIdx2] << " " // total mass flux
+ << outputVector[interfaceVertIdx].defect[contiWEqIdx2] << " " // total flux of component
+ << outputVector[interfaceVertIdx].defect[energyEqIdx2] << " " // total heat flux
+ << std::endl;
+ }
+ outfile.close();
+ }
+ if (shouldWriteVaporFlux()){
+ Scalar sumWaterFluxInLiquidPhase = sumVaporFluxes - sumWaterFluxInGasPhase;
+ fluxFile_ << sumVaporFluxes << "; "
+ << sumWaterFluxInLiquidPhase << "; "
+ << sumWaterFluxInGasPhase << "; "
+ << sumEnergyFluxes
+ << std::endl;
+ }
+ }
+
+ Scalar computeAdvectiveVaporFluxes1(const ElementVolumeVariables1& elemVolVars1,
+ const BoundaryVariables1& boundaryVars1,
+ int vertInElem1)
+ {
+ Scalar advFlux = elemVolVars1[vertInElem1].density() *
+ elemVolVars1[vertInElem1].fluidState().massFraction(phaseIdx, transportCompIdx1) *
+ boundaryVars1.normalVelocity();
+ return advFlux;
+ }
+
+ Scalar computeDiffusiveVaporFluxes1(const ElementVolumeVariables1& elemVolVars1,
+ const BoundaryVariables1& boundaryVars1,
+ int vertInElem1)
+ {
+ Scalar diffFlux = boundaryVars1.moleFractionGrad(transportCompIdx1) *
+ boundaryVars1.face().normal *
+ boundaryVars1.diffusionCoeff(transportCompIdx1) *
+ boundaryVars1.molarDensity() *
+ FluidSystem::molarMass(transportCompIdx1);
+ return diffFlux;
+ }
+
+ //! \copydoc Dumux::CoupledProblem::shouldWriteRestartFile()
+ bool shouldWriteRestartFile() const
+ {
+ if ((this->timeManager().timeStepIndex() > 0 &&
+ (this->timeManager().timeStepIndex() % freqRestart_ == 0))
+ // also write a restart file at the end of each episode
+ || this->timeManager().episodeWillBeOver())
+ return true;
+ return false;
+ }
+
+
+ //! \copydoc Dumux::CoupledProblem::shouldWriteOutput()
+ bool shouldWriteOutput() const
+ {
+ if (this->timeManager().timeStepIndex() % freqOutput_ == 0
+ || this->timeManager().episodeWillBeOver())
+ return true;
+ return false;
+ }
+
+ /*!
+ * \brief Returns true if a file with the fluxes across the
+ * free-flow -- porous-medium interface should be written.
+ */
+ bool shouldWriteFluxFile() const
+ {
+ if (this->timeManager().timeStepIndex() % freqFluxOutput_ == 0
+ || this->timeManager().episodeWillBeOver())
+ return true;
+ return false;
+ }
+
+ /*!
+ * \brief Returns true if the summarized vapor fluxes
+ * across the free-flow -- porous-medium interface,
+ * representing the evaporation rate (related to the
+ * interface area), should be written.
+ */
+ bool shouldWriteVaporFlux() const
+ {
+ if (this->timeManager().timeStepIndex() % freqVaporFluxOutput_ == 0
+ || this->timeManager().episodeWillBeOver())
+ return true;
+ return false;
+
+ }
+
+ Stokes2cniSubProblem& stokes2cniProblem()
+ { return this->sdProblem1(); }
+ const Stokes2cniSubProblem& stokes2cniProblem() const
+ { return this->sdProblem1(); }
+
+ TwoPTwoCNISubProblem& twoPtwoCNIProblem()
+ { return this->sdProblem2(); }
+ const TwoPTwoCNISubProblem& twoPtwoCNIProblem() const
+ { return this->sdProblem2(); }
+
+private:
+ typename MDGrid::SubDomainType stokes2cni_;
+ typename MDGrid::SubDomainType twoPtwoCNI_;
+
+ unsigned counter_;
+ unsigned freqRestart_;
+ unsigned freqOutput_;
+ unsigned freqFluxOutput_;
+ unsigned freqVaporFluxOutput_;
+
+ Scalar interface_;
+ Scalar noDarcyX_;
+ Scalar episodeLength_;
+ Scalar initializationTime_;
+ Scalar dtInit_;
+
+ template <int numEq>
+ struct InterfaceFluxes
+ {
+ unsigned count;
+ unsigned interfaceVertex;
+ unsigned globalIdx;
+ Scalar xCoord;
+ Scalar yCoord;
+ Dune::FieldVector<Scalar, numEq> defect;
+
+ InterfaceFluxes()
+ {
+ count = 0;
+ interfaceVertex = 0;
+ globalIdx = 0;
+ xCoord = 0.0;
+ yCoord = 0.0;
+ defect = 0.0;
+ }
+ };
+ std::ofstream fluxFile_;
+};
+
+} //end namespace
+
+#endif
diff --git a/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cnispatialparams.hh b/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cnispatialparams.hh
new file mode 100644
index 0000000000..3c1973350b
--- /dev/null
+++ b/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cnispatialparams.hh
@@ -0,0 +1,371 @@
+// -*- 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/>. *
+ *****************************************************************************/
+#ifndef DUMUX_TWOCNISTOKES2P2CNISPATIALPARAMS_HH
+#define DUMUX_TWOCNISTOKES2P2CNISPATIALPARAMS_HH
+
+#include <dune/grid/io/file/vtk/common.hh>
+
+#include <dumux/material/spatialparams/implicitspatialparams.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+//#include <dumux/io/plotfluidmatrixlaw.hh>
+
+namespace Dumux
+{
+
+//forward declaration
+template<class TypeTag>
+class TwoCNIStokesTwoPTwoCNISpatialParams;
+
+namespace Properties
+{
+// The spatial parameters TypeTag
+NEW_TYPE_TAG(TwoCNIStokesTwoPTwoCNISpatialParams);
+
+// Set the spatial parameters
+SET_TYPE_PROP(TwoCNIStokesTwoPTwoCNISpatialParams, SpatialParams,
+ TwoCNIStokesTwoPTwoCNISpatialParams<TypeTag>);
+
+// Set the material Law
+SET_PROP(TwoCNIStokesTwoPTwoCNISpatialParams, MaterialLaw)
+{
+private:
+ // define the material law which is parameterized by effective
+ // saturations
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef RegularizedVanGenuchten<Scalar> EffMaterialLaw;
+public:
+ // define the material law parameterized by absolute saturations
+ typedef EffToAbsLaw<EffMaterialLaw> type;
+};
+}
+
+
+/*!
+ * \brief docme
+ */
+template<class TypeTag>
+class TwoCNIStokesTwoPTwoCNISpatialParams : public ImplicitSpatialParams<TypeTag>
+{
+ typedef ImplicitSpatialParams<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GridView::ctype CoordScalar;
+
+ enum {
+ dim=GridView::dimension,
+ dimWorld=GridView::dimensionworld
+ };
+
+ typedef Dune::FieldVector<CoordScalar,dim> LocalPosition;
+ typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<CoordScalar,dimWorld> DimVector;
+
+ typedef typename GridView::IndexSet IndexSet;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::template Codim<0>::Entity Element;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+
+ typedef std::vector<Scalar> PermeabilityType;
+ typedef typename GET_PROP(TypeTag, ParameterTree) ParameterTree;
+ typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
+
+public:
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+ typedef typename MaterialLaw::Params MaterialLawParams;
+ typedef std::vector<MaterialLawParams> MaterialLawParamsVector;
+
+ /*!
+ * \brief docme
+ *
+ * \param gridView The GridView which is used by the problem.
+ *
+ */
+ TwoCNIStokesTwoPTwoCNISpatialParams(const GridView& gridView)
+ : ParentType(gridView),
+ permeability_(gridView.size(dim), 0.0),
+ vanGenuchtenAlpha_(gridView.size(dim), 0.0),
+ indexSet_(gridView.indexSet())
+ {
+ try
+ {
+ soilType_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, SpatialParams, SoilType);
+ xMaterialInterface_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, MaterialInterfaceX);
+
+ // porosities
+ coarsePorosity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Coarse, Porosity1);
+ mediumPorosity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Porosity2);
+ finePorosity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Fine, Porosity3);
+
+ // intrinsic permeabilities
+ coarsePermeability_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Coarse, Permeability1);
+ mediumPermeability_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Permeability2);
+ finePermeability_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Fine, Permeability3);
+
+ // thermal conductivity of the solid material
+ coarseLambdaSolid_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Coarse, LambdaSolid1);
+ mediumLambdaSolid_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, LambdaSolid2);
+ fineLambdaSolid_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Fine, LambdaSolid3);
+
+ if (soilType_ != 0)
+ {
+ // residual saturations
+ coarseParams_.setSwr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Coarse, Swr1));
+ coarseParams_.setSnr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Coarse, Snr1));
+ mediumParams_.setSwr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Swr2));
+ mediumParams_.setSnr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Snr2));
+ fineParams_.setSwr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Fine, Swr3));
+ fineParams_.setSnr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Fine, Snr3));
+
+ // parameters for the vanGenuchten law
+ coarseParams_.setVgAlpha(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Coarse, VgAlpha1));
+ coarseParams_.setVgn(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Coarse, VgN1));
+ mediumParams_.setVgAlpha(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, VgAlpha2));
+ mediumParams_.setVgn(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, VgN2));
+ fineParams_.setVgAlpha(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Fine, VgAlpha3));
+ fineParams_.setVgn(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Fine, VgN3));
+ }
+ }
+ catch (Dumux::ParameterException &e) {
+ std::cerr << e << ". Abort!\n";
+ exit(1) ;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ exit(1);
+ }
+ }
+
+ ~TwoCNIStokesTwoPTwoCNISpatialParams()
+ { }
+
+ /*!
+ * \brief Apply the intrinsic permeability tensor to a pressure
+ * potential gradient.
+ *
+ * \param element The current finite element
+ * \param fvGeometry The current finite volume geometry of the element
+ * \param scvIdx The index sub-control volume face where the
+ * intrinsic velocity ought to be calculated.
+ */
+ const Scalar intrinsicPermeability(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ if (checkSoilType(globalPos) == 1)
+ return coarsePermeability_;
+ if (checkSoilType(globalPos) == 2)
+ return mediumPermeability_;
+ if (checkSoilType(globalPos) == 3)
+ return finePermeability_;
+ if (checkSoilType(globalPos) == 4)
+ return finePermeability_;
+ else
+ return mediumPermeability_;
+ }
+
+ //! \copydoc Dumux::ImplicitSpatialParamsOneP::porosity()
+ Scalar porosity(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ if (checkSoilType(globalPos) == 1)
+ return coarsePorosity_;
+ if (checkSoilType(globalPos) == 2)
+ return mediumPorosity_;
+ if (checkSoilType(globalPos) == 3)
+ return finePorosity_;
+ else
+ return mediumPorosity_;
+ }
+
+
+ //! \copydoc Dumux::ImplicitSpatialParams::materialLawParams()
+ const MaterialLawParams& materialLawParams(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+ if (checkSoilType(globalPos)==1)
+ return coarseParams_;
+ if (checkSoilType(globalPos)==2)
+ return mediumParams_;
+ if (checkSoilType(globalPos)==3)
+ return fineParams_;
+// if (checkSoilType(globalPos)==4)
+// return leverettJParams_;
+ else
+ return mediumParams_;
+ }
+
+
+ /*!
+ * \brief Returns the heat capacity \f$[J/m^3 K]\f$ of the rock matrix.
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume where
+ * the heat capacity needs to be defined
+ */
+ Scalar heatCapacity(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ return
+ // 830 for sand
+ 790 // specific heat capacity of granite [J / (kg K)]
+ * 2700 // density of granite [kg/m^3]
+ * (1 - porosity(element, fvGeometry, scvIdx));
+ }
+
+ /*!
+ * \brief docme
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry
+ * \param scvIdx The local subcontrolvolume index
+ *
+ */
+ Scalar thermalConductivitySolid(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ const GlobalPosition &pos = element.geometry().corner(scvIdx);
+
+ if (checkSoilType(pos) == 1)
+ return coarseLambdaSolid_;
+ if (checkSoilType(pos) == 2)
+ return mediumLambdaSolid_;
+ if (checkSoilType(pos) == 3)
+ return fineLambdaSolid_;
+ else
+ return mediumLambdaSolid_;
+ }
+
+ /*!
+ * \brief docme
+ *
+ * \param pos The position of the boundary face's integration point in global coordinates
+ *
+ */
+ const unsigned checkSoilType(const GlobalPosition &pos) const
+ {
+
+// one soil, which can be chosen as runtime parameter:
+// 1: coarse
+// 2: medium
+// 3: fine
+// 4: LeverettJ (x < xMaterialInterface)
+ return soilType_;
+ }
+
+ // this is called from the coupled problem
+ // and creates a gnuplot output of the Pc-Sw curve
+ /*!
+ * \brief docme
+ */
+ void plotMaterialLaw()
+ {
+// if (soilType_ == 0)
+// {
+// std::cout << "Material law plot not possible for heterogeneous media!\n";
+// return;
+// }
+// if (GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, SpatialParams.Coarse, PlotMaterialLaw1))
+// {
+// PlotFluidMatrixLaw<MaterialLaw> coarseFluidMatrixLaw_;
+// coarseFluidMatrixLaw_.plotpC(coarseParams_,
+// GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Coarse, Swr1),
+// 1.0 - GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Coarse, Snr1),
+// "pcSw_coarse");
+// }
+// if (GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, SpatialParams.Medium, PlotMaterialLaw2))
+// {
+// PlotFluidMatrixLaw<MaterialLaw> mediumFluidMatrixLaw_;
+// mediumFluidMatrixLaw_.plotpC(mediumParams_,
+// GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Swr2),
+// 1.0 - GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Snr2),
+// "pcSw_medium");
+// PlotFluidMatrixLaw<MaterialLaw> mediumFluidMatrixLaw2_;
+// mediumFluidMatrixLaw_.plotkrw(mediumParams_,
+// GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Swr2),
+// 1.0 - GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Snr2));
+// PlotFluidMatrixLaw<MaterialLaw> mediumFluidMatrixLaw3_;
+// mediumFluidMatrixLaw_.plotkrn(mediumParams_,
+// GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Swr2),
+// 1.0 - GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Snr2));
+// }
+// if (GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, SpatialParams.Fine, PlotMaterialLaw3))
+// {
+// PlotFluidMatrixLaw<MaterialLaw> fineFluidMatrixLaw_;
+// fineFluidMatrixLaw_.plotpC(fineParams_,
+// GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Fine, Swr3),
+// 1.0 - GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Fine, Snr3),
+// "pcSw_fine");
+// }
+// if (GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, SpatialParams.LeverettJ, PlotMaterialLawJ))
+// {
+// PlotFluidMatrixLaw<MaterialLaw> leverettJFluidMatrixLaw_;
+// leverettJFluidMatrixLaw_.plotpC(leverettJParams_,
+// GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Swr2),
+// 1.0 - GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Snr2),
+// "pcSw_leverett");
+// }
+ }
+
+private:
+
+ unsigned soilType_;
+
+ Scalar coarsePermeability_;
+ Scalar mediumPermeability_;
+ Scalar finePermeability_;
+
+ Scalar coarsePorosity_;
+ Scalar mediumPorosity_;
+ Scalar finePorosity_;
+
+ Scalar coarseLambdaSolid_;
+ Scalar mediumLambdaSolid_;
+ Scalar fineLambdaSolid_;
+
+ Scalar xMaterialInterface_;
+ MaterialLawParamsVector materialLawParams_;
+ PermeabilityType permeability_;
+ PermeabilityType vanGenuchtenAlpha_;
+ const IndexSet& indexSet_;
+
+ MaterialLawParams coarseParams_;
+ MaterialLawParams mediumParams_;
+ MaterialLawParams fineParams_;
+};
+
+} // end namespace
+#endif
+
diff --git a/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh b/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh
new file mode 100644
index 0000000000..772e1988c8
--- /dev/null
+++ b/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh
@@ -0,0 +1,675 @@
+// -*- 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 Porous-medium subproblem with coupling at the top boundary.
+ */
+#ifndef DUMUX_2P2CNISUB_PROBLEM_HH
+#define DUMUX_2P2CNISUB_PROBLEM_HH
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dumux/implicit/common/implicitporousmediaproblem.hh>
+#include <dumux/implicit/2p2cni/2p2cnimodel.hh>
+#include <dumux/multidomain/couplinglocalresiduals/2p2cnicouplinglocalresidual.hh>
+#include <dumux/multidomain/common/subdomainpropertydefaults.hh>
+#include <dumux/multidomain/common/multidomainlocaloperator.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
+#include <dune/pdelab/finiteelementmap/conformingconstraints.hh>
+
+#include "2cnistokes2p2cniproblem.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class TwoPTwoCNISubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(TwoPTwoCNISubProblem,
+ INHERITS_FROM(BoxTwoPTwoCNI, SubDomain, TwoCNIStokesTwoPTwoCNISpatialParams));
+
+// Set the problem property
+SET_TYPE_PROP(TwoPTwoCNISubProblem, Problem, TwoPTwoCNISubProblem<TTAG(TwoPTwoCNISubProblem)>);
+
+//! Use the 2p2cni local jacobian operator for the 2p2cniCoupling model
+SET_TYPE_PROP(TwoPTwoCNISubProblem, LocalResidual, TwoPTwoCNICouplingLocalResidual<TypeTag>);
+
+// choose pn and Sw as primary variables
+SET_INT_PROP(TwoPTwoCNISubProblem, Formulation, TwoPTwoCFormulation::pnsw);
+
+// the gas component balance (air) is replaced by the total mass balance
+SET_PROP(TwoPTwoCNISubProblem, ReplaceCompEqIdx)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ static const int value = Indices::contiNEqIdx;
+};
+
+// Set the model parameter group for the TypeTag
+//SET_PROP(TwoPTwoCNISubProblem, ModelParameterGroup)
+//{ static const char *value() { return "PM"; }; };
+
+// set the constraints for multidomain / pdelab
+SET_TYPE_PROP(TwoPTwoCNISubProblem, Constraints,
+ Dune::PDELab::NonoverlappingConformingDirichletConstraints);
+
+// set the grid operator
+SET_PROP(TwoPTwoCNISubProblem, GridOperator)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, ConstraintsTrafo) ConstraintsTrafo;
+ typedef typename GET_PROP_TYPE(TypeTag, GridFunctionSpace) GridFunctionSpace;
+ typedef Dumux::PDELab::MultiDomainLocalOperator<TypeTag> LocalOperator;
+
+ enum{numEq = GET_PROP_VALUE(TypeTag, NumEq)};
+
+public:
+ typedef Dune::PDELab::GridOperator<GridFunctionSpace,
+ GridFunctionSpace, LocalOperator,
+ Dune::PDELab::ISTLBCRSMatrixBackend<numEq, numEq>,
+ Scalar, Scalar, Scalar,
+ ConstraintsTrafo,
+ ConstraintsTrafo,
+ true> type;
+};
+
+// the fluidsystem is set in the coupled problem
+SET_PROP(TwoPTwoCNISubProblem, FluidSystem)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, MultiDomainTypeTag) CoupledTypeTag;
+ typedef typename GET_PROP_TYPE(CoupledTypeTag, FluidSystem) FluidSystem;
+public:
+ typedef FluidSystem type;
+};
+
+
+//! Somerton is used as model to compute the effective thermal heat conductivity
+SET_PROP(TwoPTwoCNISubProblem, ThermalConductivityModel)
+{ private :
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ public:
+ typedef ThermalConductivitySomerton<Scalar> type;
+};
+
+// use formulation based on mass fractions
+SET_BOOL_PROP(TwoPTwoCNISubProblem, UseMoles, false);
+
+// enable/disable velocity output
+SET_BOOL_PROP(TwoPTwoCNISubProblem, VtkAddVelocity, true);
+
+// Enable gravity
+SET_BOOL_PROP(TwoPTwoCNISubProblem, ProblemEnableGravity, true);
+}
+
+
+template <class TypeTag = TTAG(TwoPTwoCNISubProblem) >
+class TwoPTwoCNISubProblem : public ImplicitPorousMediaProblem<TypeTag>
+{
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GridView::Grid Grid;
+
+ typedef TwoPTwoCNISubProblem<TypeTag> ThisType;
+ typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
+
+ // copy some indices for convenience
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ // the type tag of the coupled problem
+ typedef typename GET_PROP_TYPE(TypeTag, MultiDomainTypeTag) CoupledTypeTag;
+
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+ enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+ enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+ enum { // the equation indices
+ contiTotalMassIdx = Indices::contiNEqIdx,
+ contiWEqIdx = Indices::contiWEqIdx,
+ energyEqIdx = Indices::energyEqIdx
+ };
+ enum { // the indices of the primary variables
+ pNIdx = Indices::pressureIdx,
+ sWIdx = Indices::switchIdx,
+ temperatureIdx = Indices::temperatureIdx
+ };
+ enum { // the indices for the phase presence
+ wCompIdx = Indices::wCompIdx,
+ nCompIdx = Indices::nCompIdx
+ };
+ enum { // the indices for the phase presence
+ wPhaseOnly = Indices::wPhaseOnly,
+ nPhaseOnly = Indices::nPhaseOnly,
+ bothPhases = Indices::bothPhases
+ };
+ enum {
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx
+ };
+ enum { // grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
+// typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+
+ // only required for radiation
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+ //////////////
+
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+ /*!
+ * \brief docme
+ *
+ * \param timeManager The TimeManager which is used by the simulation
+ * \param gridView The simulation's idea about physical space
+ *
+ */
+public:
+ TwoPTwoCNISubProblem(TimeManager &timeManager, const GridView &gridView)
+ : ParentType(timeManager, gridView)
+ {
+ try
+ {
+ Scalar noDarcyX = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, NoDarcyX);
+ Scalar xMin = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMin);
+
+ bboxMin_[0] = std::max(xMin,noDarcyX);
+ bboxMax_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMax);
+
+ bboxMin_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, YMin);
+ bboxMax_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePos);
+
+ runUpDistanceX_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, RunUpDistanceX); // first part of the interface without coupling
+ xMaterialInterface_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, MaterialInterfaceX);
+ initializationTime_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, InitTime);
+
+ refTemperature_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, PorousMedium, RefTemperaturePM);
+ refPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, PorousMedium, RefPressurePM);
+ initialSw1_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, PorousMedium, InitialSw1);
+ initialSw2_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, PorousMedium, InitialSw2);
+
+// porosity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Porosity2);
+
+ freqMassOutput_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqMassOutput);
+
+ storageLastTimestep_ = Scalar(0);
+ lastMassOutputTime_ = Scalar(0);
+
+ outfile.open("evaporation.out");
+ outfile << "time; massChange; evaporationRate; energyChange; massH2O" << std::endl;
+ }
+ catch (Dumux::ParameterException &e) {
+ std::cerr << e << ". Abort!\n";
+ exit(1) ;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ exit(1);
+ }
+ }
+
+ ~TwoPTwoCNISubProblem()
+ {
+ outfile.close();
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ //! \copydoc Dumux::ImplicitProblem::name()
+ const std::string &name() const
+ { return GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Vtk, NamePM); }
+
+ //! \copydoc Dumux::ImplicitProblem::init()
+ void init()
+ {
+ // set the initial condition of the model
+ ParentType::init();
+
+ this->model().globalStorage(storageLastTimestep_);
+
+ const int blModel = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, FreeFlow, UseBoundaryLayerModel);
+ std::cout << "Using boundary layer model " << blModel << std::endl;
+ }
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ //! \copydoc Dumux::ImplicitProblem::boundaryTypes()
+ void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
+ {
+ const GlobalPosition globalPos = vertex.geometry().center();
+ Scalar time = this->timeManager().time();
+
+ values.setAllNeumann();
+
+ if (onLowerBoundary_(globalPos))
+ {
+// values.setAllDirichlet();
+// values.setDirichlet(pressureIdx, contiTotalMassIdx);
+ values.setDirichlet(temperatureIdx, energyEqIdx);
+ }
+
+ if (onUpperBoundary_(globalPos))
+ {
+ if (time > initializationTime_)
+ {
+ if (globalPos[0] > runUpDistanceX_ - eps_)
+ {
+ values.setAllCouplingInflow();
+ }
+ else
+ values.setAllNeumann();
+ }
+ else
+ {
+// values.setAllDirichlet();
+ values.setAllNeumann();
+ }
+ }
+ }
+
+ //! \copydoc Dumux::ImplicitProblem::dirichlet()
+ void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
+ {
+ const GlobalPosition globalPos = vertex.geometry().center();
+
+ initial_(values, globalPos);
+ }
+
+ //! \copydoc Dumux::ImplicitProblem::neumann()
+ void neumann(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const Intersection &is,
+ const int scvIdx,
+ const int boundaryFaceIdx) const
+ {
+ values = 0.;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+ //! \copydoc Dumux::ImplicitProblem::solDependentSource()
+ void solDependentSource(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx,
+ const ElementVolumeVariables &elemVolVars) const
+ {
+ values = 0;
+// const Scalar time = this->timeManager().time();
+// const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+//
+// const Scalar temperatureOut = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefTemperature);
+//
+//// const Scalar S = variationSun_(0); //Solar irradiance (Yamanaka 1998)
+// const Scalar sigma = 5.671e-8; //Boltzmann-Constant
+// const Scalar e_a = atmosphericEmissivity(); //atmosphere emissivity (Brutsaert 1982)
+// const Scalar e_s = surfaceEmissivity_(elemVolVars, scvIdx);
+// const Scalar a = albedo_(elemVolVars, scvIdx);
+// const Scalar T_a = 298.15;//variationT_a_(6); //atmosphere Temperature (Yamanaka 1998)
+// const Scalar T_s = elemVolVars[scvIdx].temperature(); //surface Temperature
+//
+// const Scalar vertexPosition = 0.249248;//1.0; //GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InterfacePos);
+//
+// //For d depth with b cells!!
+//// Scalar radiation = S*(1.0-a)+
+// Scalar radiation = sigma*e_s*(e_a*std::pow(T_a,4.0)-(std::pow(T_s,4.0))); //Equation from Novak 2010
+//
+// if (time > initializationTime_
+// && globalPos[1] > (vertexPosition - 0.00002)
+// && globalPos[1] < (vertexPosition + 0.00002))
+// values[energyEqIdx] += radiation/cellHeight_();
+ }
+
+
+ //! \copydoc Dumux::ImplicitProblem::initial()
+ void initial(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ values = 0.;
+
+ initial_(values, globalPos);
+ }
+
+ /*!
+ * \brief Return the initial phase state inside a control volume.
+ *
+ * \param vert The vertex
+ * \param globalIdx The index of the global vertex
+ * \param globalPos The global position
+ */
+ int initialPhasePresence(const Vertex &vert,
+ const int &globalIdx,
+ const GlobalPosition &globalPos) const
+ {
+ return bothPhases;
+ }
+
+
+// /*!
+// * \brief docme
+// */
+// Scalar globalRadiation()
+// {
+// FVElementGeometry fvGeometry;
+// ElementVolumeVariables elemVolVars;
+// Scalar source = 0;
+// PrimaryVariables values(0);
+//
+// ElementIterator elemIt = this->gridView().template begin<0>();
+// const ElementIterator elemEndIt = this->gridView().template end<0>();
+// for (elemIt = this->gridView().template begin<0>(); elemIt != elemEndIt; ++elemIt)
+// {
+// if(elemIt->partitionType() == Dune::InteriorEntity)
+// {
+// fvGeometry.update(this->gridView(), *elemIt);
+// elemVolVars.update(*this, *elemIt, fvGeometry, false);
+//
+// for (int scvIdx = 0; scvIdx < elemIt->template count<dim>(); ++scvIdx)
+//
+// {
+// boxSDSource(values,
+// *elemIt,
+// fvGeometry,
+// scvIdx,
+// elemVolVars);
+// Valgrind::CheckDefined(values);
+//
+// // thickness of the domain, for 2D usually 1m
+// Scalar extrusionFactor =
+// elemVolVars[scvIdx].extrusionFactor();
+// values[energyEqIdx] *= fvGeometry.subContVol[scvIdx].volume * extrusionFactor;
+//
+// source += values[energyEqIdx];
+// }
+// }
+// }
+//
+// if (this->gridView().comm().size() > 1)
+// source = this->gridView().comm().sum(source);
+//
+// return source;
+// }
+
+ //! \copydoc Dumux::ImplicitProblem::postTimeStep()
+ void postTimeStep()
+ {
+ // Calculate masses
+ PrimaryVariables storage;
+
+ this->model().globalStorage(storage);
+ const Scalar time = this->timeManager().time() + this->timeManager().timeStepSize();
+
+ // Write mass balance information for rank 0
+ if (this->gridView().comm().rank() == 0)
+ {
+ if (this->timeManager().timeStepIndex() % freqMassOutput_ == 0
+ || this->timeManager().episodeWillBeOver())
+ {
+ PrimaryVariables evaporationRate(0.);
+ evaporationRate = storageLastTimestep_ - storage;
+
+ assert(time - lastMassOutputTime_ != 0);
+ evaporationRate /= (time - lastMassOutputTime_);
+ // 2d: interface length has to be accounted for
+ // in order to obtain kg/m²s
+ evaporationRate /= (bboxMax_[0]-bboxMin_[0]);
+
+ std::cout << "TotalMass: " << storage[contiTotalMassIdx]
+ << " MassH2O: " << storage[contiWEqIdx]
+ << " IntEnergy: " << storage[energyEqIdx]
+ << " EvaporationRate: " << evaporationRate[contiWEqIdx]
+ << " Time: " << time <<std::endl;
+ if (this->timeManager().time() != 0.)
+ outfile << time <<
+ "; " << evaporationRate[contiTotalMassIdx] <<
+ "; " << evaporationRate[contiWEqIdx] <<
+ "; " << evaporationRate[energyEqIdx] <<
+ "; " << storage[contiWEqIdx] <<
+ std::endl;
+
+ storageLastTimestep_ = storage;
+ lastMassOutputTime_ = time;
+ }
+ }
+ }
+
+ /*!
+ * \brief Determine if we are on a corner of the grid
+ *
+ * \param globalPos The global position
+ *
+ */
+ bool isCornerPoint(const GlobalPosition &globalPos)
+ {
+ return ((onLeftBoundary_(globalPos) && onLowerBoundary_(globalPos)) ||
+ (onLeftBoundary_(globalPos) && onUpperBoundary_(globalPos)) ||
+ (onRightBoundary_(globalPos) && onLowerBoundary_(globalPos)) ||
+ (onRightBoundary_(globalPos) && onUpperBoundary_(globalPos)));
+ }
+
+ // required in case of mortar coupling
+ // otherwise it should return false
+ /*!
+ * \brief docme
+ *
+ * \param global Pos The global position
+ *
+ */
+ bool isInterfaceCornerPoint(const GlobalPosition &globalPos) const
+ { return false; }
+
+ // \}
+
+private:
+ // internal method for the initial condition (reused for the
+ // dirichlet conditions!)
+
+ void initial_(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ //TODO: call density from fluidsystem
+// const MaterialLawParams &materialParams =
+// this->spatialParams().materialLawParams(globalPos);
+// Scalar pC = MaterialLaw::pC(materialParams, initialSw1_);
+
+ values[pNIdx] = refPressure_
+ + 1000.*this->gravity()[1]*(globalPos[1]-bboxMax_[1]);
+// + pC;
+
+ if (globalPos[0] < xMaterialInterface_)
+ values[sWIdx] = initialSw1_;
+ else
+ values[sWIdx] = initialSw2_;
+ values[temperatureIdx] = refTemperature_;
+ // if (onUpperBoundary_(globalPos) && onLeftBoundary_(globalPos))
+// {
+// values[sWIdx] = 0.99e-4;
+// values[temperatureIdx] = 283.55;
+// }
+
+// values[pNIdx] = 1e5*(2.0 - globalPos[0]);
+// if (onLeftBoundary_(globalPos))
+// values[sWIdx] = 0.9e-4;
+// else
+// values[sWIdx] = 1e-4;
+// values[temperatureIdx] = 283.15;
+ }
+
+ /// RADIATION stuff
+// const Scalar albedo_(const ElementVolumeVariables &elemVolVars, const int scvIdx) const //albedo
+// {
+// if (elemVolVars[scvIdx].saturation(wPhaseIdx)*porosity_ > 0.3)
+// return 0.075;
+// if (0.04 <= elemVolVars[scvIdx].saturation(wPhaseIdx)*porosity_ &&
+// elemVolVars[scvIdx].saturation(wPhaseIdx)*porosity_ <= 0.3)
+// return 0.1846-0.3654*elemVolVars[scvIdx].saturation(wPhaseIdx)*porosity_;
+// if (elemVolVars[scvIdx].saturation(wPhaseIdx)*porosity_ > 0.4)
+// return 0.17;
+// return 0.17;
+// }
+//
+// const Scalar surfaceEmissivity_(const ElementVolumeVariables &elemVolVars, const int scvIdx) const //surface emissivity
+// {
+// if (elemVolVars[scvIdx].saturation(wPhaseIdx)*porosity_ < 0.3)
+// return 0.93+0.1333*elemVolVars[scvIdx].saturation(wPhaseIdx)*porosity_;
+// if (elemVolVars[scvIdx].saturation(wPhaseIdx)*porosity_ >= 0.3)
+// return 0.97;
+// return 0.;
+// }
+//
+// //to simulate the Solar irradiance for a number of days (definable in the input file)
+// const Scalar variationSun_(const Scalar value) const
+// {
+// const Scalar time = this->timeManager().time() + this->timeManager().timeStepSize();
+// const Scalar days = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, days);
+//
+// Scalar solar = 0;
+// for (int n=0; n<days; ++n)
+// {
+// if (time > 21600 + n*86400 && time < 64800 + n*86400)
+// solar = cos(2*M_PI*((time+43200)/(86400))) * value;
+// }
+// return solar;
+// }
+//
+// //to simulate the atmospheric Temperature for a number of days (definable in the input file)
+// const Scalar variationT_a_(const Scalar value) const
+// {
+// const Scalar time = this->timeManager().time() + this->timeManager().timeStepSize();
+// const Scalar T = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefTemperature);
+// {
+// return T+cos(2*M_PI*((time+36000)/(86400))) * value;
+// }
+// std::cerr << "error with T_a at time " << time << std::endl;
+// }
+//
+// //to calculate the atmospheric Emissivity depending on the Temperature near the Surface
+// const Scalar atmosphericEmissivity() const
+// {
+// const Scalar T_a = variationT_a_(6);
+//
+// return 1.24*std::pow((((Dumux::H2O<Scalar>::vaporPressure(T_a)/100)/T_a)),(1.0/7.0));
+// }
+//
+// const Scalar cellHeight_() const
+// {
+// return 0.0004575;
+// }
+ /////////
+
+// TODO: what should be used as reference mass fraction??
+// void updateFluidStateForBC_(FluidState& fluidState) const
+// {
+// for (unsigned phaseIdx=0; phaseIdx<numPhases; ++phaseIdx)
+// {
+// fluidState.setTemperature(refTemperature_);
+// fluidState.setPressure(phaseIdx, refPressure_);
+//
+// Scalar massFraction[numComponents];
+// massFraction[wCompIdx] = refMassfrac_;
+// massFraction[nCompIdx] = 1 - massFraction[wCompIdx];
+//
+// // calculate average molar mass of the gas phase
+// Scalar M1 = FluidSystem::molarMass(wCompIdx);
+// Scalar M2 = FluidSystem::molarMass(nCompIdx);
+// Scalar X2 = massFraction[nCompIdx];
+// Scalar massToMoleDenominator = M2 + X2*(M1 - M2);
+//
+// fluidState.setMoleFraction(phaseIdx, wCompIdx, massFraction[wCompIdx]*M2/massToMoleDenominator);
+// fluidState.setMoleFraction(phaseIdx, nCompIdx, massFraction[nCompIdx]*M1/massToMoleDenominator);
+// }
+// }
+
+
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < bboxMin_[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > bboxMax_[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < bboxMin_[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > bboxMax_[1] - eps_; }
+
+ bool onBoundary_(const GlobalPosition &globalPos) const
+ {
+ return (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)
+ || onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos));
+ }
+
+ static constexpr Scalar eps_ = 1e-8;
+ GlobalPosition bboxMin_;
+ GlobalPosition bboxMax_;
+ Scalar xMaterialInterface_;
+
+ int freqMassOutput_;
+
+ PrimaryVariables storageLastTimestep_;
+ Scalar lastMassOutputTime_;
+
+ Scalar refTemperature_;
+ Scalar refPressure_;
+ Scalar initialSw1_;
+ Scalar initialSw2_;
+
+// Scalar porosity_;
+ Scalar runUpDistanceX_;
+ Scalar initializationTime_;
+ std::ofstream outfile;
+};
+} //end namespace
+
+#endif
diff --git a/test/multidomain/2cnistokes2p2cni/Makefile.am b/test/multidomain/2cnistokes2p2cni/Makefile.am
new file mode 100644
index 0000000000..e6520df3f3
--- /dev/null
+++ b/test/multidomain/2cnistokes2p2cni/Makefile.am
@@ -0,0 +1,5 @@
+# tests where program to build and program to run are equal
+check_PROGRAMS = test_2cnistokes2p2cni
+test_2cnistokes2p2cni_SOURCES = test_2cnistokes2p2cni.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/multidomain/2cnistokes2p2cni/grids/interfacedomain.dgf b/test/multidomain/2cnistokes2p2cni/grids/interfacedomain.dgf
new file mode 100644
index 0000000000..b0290172cb
--- /dev/null
+++ b/test/multidomain/2cnistokes2p2cni/grids/interfacedomain.dgf
@@ -0,0 +1,15 @@
+DGF
+Interval
+0 0 % first corner
+0.25 0.5 % second corner
+1 1 % 1 cells in x and 1 in y direction
+#
+
+Cube
+0 1 2 3
+#
+
+BOUNDARYDOMAIN
+default 1 % all boundaries have id 1
+#BOUNDARYDOMAIN
+# unitcube.dgf
diff --git a/test/multidomain/2cnistokes2p2cni/grids/test_2cnistokes2p2cni.dgf b/test/multidomain/2cnistokes2p2cni/grids/test_2cnistokes2p2cni.dgf
new file mode 100644
index 0000000000..f13c8f14e2
--- /dev/null
+++ b/test/multidomain/2cnistokes2p2cni/grids/test_2cnistokes2p2cni.dgf
@@ -0,0 +1,17 @@
+DGF
+Interval
+0 0 % first corner
+0.25 0.5 % second corner
+20 64 % cells in x and y direction
+#
+GridParameter
+overlap 0
+periodic
+closure green
+copies none
+heapsize 500
+
+BOUNDARYDOMAIN
+default 1 % all boundaries have id 1
+#BOUNDARYDOMAIN
+# unitcube.dgf
diff --git a/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh b/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh
new file mode 100644
index 0000000000..3602c2837e
--- /dev/null
+++ b/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh
@@ -0,0 +1,592 @@
+// -*- 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 Stokes2cniProblems
+ * @brief Definition of a non-isothermal compositional Stokes problem
+ * @author Klaus Mosthaf
+ */
+#ifndef DUMUX_STOKES2CNI_SUBPROBLEM_HH
+#define DUMUX_STOKES2CNI_SUBPROBLEM_HH
+
+#include <dune/pdelab/finiteelementmap/conformingconstraints.hh>
+
+#include <dumux/freeflow/stokesncni/stokesncnimodel.hh>
+#include <dumux/multidomain/couplinglocalresiduals/stokesncnicouplinglocalresidual.hh>
+#include <dumux/multidomain/common/subdomainpropertydefaults.hh>
+
+namespace Dumux
+{
+
+template <class TypeTag>
+class Stokes2cniSubProblem;
+
+//////////
+// Specify the properties for the Stokes problem
+//////////
+namespace Properties
+{
+NEW_TYPE_TAG(Stokes2cniSubProblem,
+ INHERITS_FROM(BoxStokesncni, SubDomain, TwoCNIStokesTwoPTwoCNISpatialParams));
+
+//SET_PROP(Stokes2cniSubProblem, ModelParameterGroup)
+//{ static const char *value() { return "FF"; }; };
+
+// Set the problem property
+SET_TYPE_PROP(Stokes2cniSubProblem, Problem, Dumux::Stokes2cniSubProblem<TypeTag>);
+
+//! Use the Stokes2cniCouplingLocalResidual for the computation of the local residual in the Stokes domain
+SET_TYPE_PROP(Stokes2cniSubProblem, LocalResidual, StokesncniCouplingLocalResidual<TypeTag>);
+
+//! Set the property for the material parameters by extracting it from the material law
+SET_PROP(Stokes2cniSubProblem, MaterialLawParams)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+public:
+ typedef typename MaterialLaw::Params type;
+};
+
+SET_TYPE_PROP(Stokes2cniSubProblem, Constraints,
+ Dune::PDELab::NonoverlappingConformingDirichletConstraints);
+
+// set the grid operator
+SET_PROP(Stokes2cniSubProblem, GridOperator)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, ConstraintsTrafo) ConstraintsTrafo;
+ typedef typename GET_PROP_TYPE(TypeTag, GridFunctionSpace) GridFunctionSpace;
+ typedef Dumux::PDELab::MultiDomainLocalOperator<TypeTag> LocalOperator;
+
+ enum{numEq = GET_PROP_VALUE(TypeTag, NumEq)};
+public:
+ typedef Dune::PDELab::GridOperator<GridFunctionSpace,
+ GridFunctionSpace, LocalOperator,
+ Dune::PDELab::ISTLBCRSMatrixBackend<numEq, numEq>,
+ Scalar, Scalar, Scalar,
+ ConstraintsTrafo, ConstraintsTrafo,
+ true> type;
+};
+
+SET_PROP(Stokes2cniSubProblem, FluidSystem)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, MultiDomainTypeTag) CoupledTypeTag;
+ typedef typename GET_PROP_TYPE(CoupledTypeTag, FluidSystem) FluidSystem;
+public:
+ typedef FluidSystem type;
+};
+
+// use formulation based on mass fractions
+SET_BOOL_PROP(Stokes2cniSubProblem, UseMoles, false);
+
+// Disable gravity in the Stokes domain
+SET_BOOL_PROP(Stokes2cniSubProblem, ProblemEnableGravity, false);
+
+// switch inertia term on or off
+SET_BOOL_PROP(Stokes2cniSubProblem, EnableNavierStokes, false);
+}
+
+/*!
+ * \ingroup Stokes2cniBoxProblems
+ * \brief Stokes2cni problem with air flowing
+ * from the left to the right.
+ *
+ * The stokes subdomain is sized 1m times 1m. The boundary conditions for the momentum balances
+ * are all set to Dirichlet. The mass balance receives
+ * outflow bcs, which are replaced in the localresidual by the sum
+ * of the two momentum balances. In the middle of the right boundary,
+ * one vertex receives Dirichlet bcs, to set the pressure level.
+ *
+ * This sub problem uses the \ref Stokes2cniModel. It is part of the 2cnistokes2p2cni model and
+ * is combined with the 2p2cnisubproblem for the Darcy domain.
+ *
+ */
+template <class TypeTag>
+class Stokes2cniSubProblem : public StokesProblem<TypeTag>
+{
+ typedef Stokes2cniSubProblem<TypeTag> ThisType;
+ typedef StokesProblem<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+ // soil parameters for beavers & joseph
+ typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+
+ enum {
+ // Number of equations and grid dimension
+ numEq = GET_PROP_VALUE(TypeTag, NumEq),
+ dim = GridView::dimension
+ };
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ enum {
+ // equation indices
+ massBalanceIdx = Indices::massBalanceIdx,
+ momentumXIdx = Indices::momentumXIdx, //!< Index of the x-component of the momentum balance
+ momentumYIdx = Indices::momentumYIdx, //!< Index of the y-component of the momentum balance
+ momentumZIdx = Indices::momentumZIdx, //!< Index of the z-component of the momentum balance
+ transportEqIdx = Indices::transportEqIdx, //!< Index of the transport equation (massfraction)
+ energyEqIdx = Indices::energyEqIdx //!< Index of the energy equation (temperature)
+ };
+ enum { // primary variable indices
+ pressureIdx = Indices::pressureIdx,
+ velocityXIdx = Indices::velocityXIdx,
+ velocityYIdx = Indices::velocityYIdx,
+ velocityZIdx = Indices::velocityZIdx,
+ massOrMoleFracIdx = Indices::massOrMoleFracIdx,
+ temperatureIdx = Indices::temperatureIdx
+ };
+ enum { phaseIdx = Indices::phaseIdx };
+ enum { numComponents = Indices::numComponents };
+ enum {
+ transportCompIdx = Indices::transportCompIdx, //!< water component index
+ phaseCompIdx = Indices::phaseCompIdx //!< air component index
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::ctype CoordScalar;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef Dune::FieldVector<CoordScalar, dim> GlobalPosition;
+
+ /*!
+ * \brief docme
+ */
+public:
+ Stokes2cniSubProblem(TimeManager &timeManager, const GridView &gridView)
+ : ParentType(timeManager, gridView),
+ spatialParams_(gridView)
+ {
+ try
+ {
+ bboxMin_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMin);
+ bboxMax_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMax);
+ bboxMin_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePos);
+ bboxMax_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, YMax);
+
+ refTemperature_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefTemperature);
+ refPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefPressure);
+ refMassfrac_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefMassfrac);
+ vxMax_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, VxMax);
+ bjSlipVel_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, BeaversJosephSlipVel);
+
+ sinusVelVar_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusVelVar);
+ sinusPVar_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusPVar);
+ sinusTVar_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusTVar);
+ sinusXVar_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusXVar);
+
+ xMaterialInterface_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, MaterialInterfaceX);
+ runUpDistanceX_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, RunUpDistanceX); // first part of the interface without coupling
+ initializationTime_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, InitTime);
+
+ alphaBJ_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, AlphaBJ);
+ }
+ catch (Dumux::ParameterException &e) {
+ std::cerr << e << ". Abort!\n";
+ exit(1) ;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ exit(1);
+ }
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ //! \copydoc Dumux::ImplicitProblem::name()
+ const std::string &name() const
+ { return GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Vtk, NameFF); }
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ //! \copydoc Dumux::ImplicitProblem::boundaryTypes()
+ void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
+ {
+ const GlobalPosition &globalPos = vertex.geometry().center();
+ const Scalar time = this->timeManager().time();
+
+ values.setAllDirichlet();
+
+ if (onUpperBoundary_(globalPos))
+ {
+ values.setNeumann(transportEqIdx);
+ values.setNeumann(energyEqIdx);
+ }
+
+ // Left inflow boundaries should be Neumann, otherwise the
+ // evaporative fluxes are much more grid dependent
+ if (onLeftBoundary_(globalPos))
+ {
+ values.setNeumann(transportEqIdx);
+ values.setNeumann(energyEqIdx);
+
+ if (onUpperBoundary_(globalPos)) // corner point
+ values.setAllDirichlet();
+ }
+
+ if (onRightBoundary_(globalPos))
+ {
+ values.setAllOutflow();
+
+ if (onUpperBoundary_(globalPos)) // corner point
+ values.setAllDirichlet();
+ }
+
+ if (onLowerBoundary_(globalPos))
+ {
+ values.setAllDirichlet();
+ if (!onLeftBoundary_(globalPos)) // is this required?
+ {
+ values.setNeumann(transportEqIdx);
+ values.setNeumann(energyEqIdx);
+ }
+
+ if (globalPos[0] > runUpDistanceX_-eps_ && time > initializationTime_)
+ {
+ values.setAllCouplingOutflow();
+// values.setCouplingInflow(energyEqIdx);
+ }
+ }
+
+ // the mass balance has to be of type outflow
+ // it does not get a coupling condition, since pn is a condition for stokes
+ values.setOutflow(massBalanceIdx);
+
+ // set pressure at one point, do NOT specify this
+ // if the Darcy domain has a Dirichlet condition for pressure
+ if (onRightBoundary_(globalPos))
+ {
+ if (time > initializationTime_)
+ values.setDirichlet(pressureIdx, massBalanceIdx);
+ else
+ if (!onLowerBoundary_(globalPos) && !onUpperBoundary_(globalPos))
+ values.setDirichlet(pressureIdx, massBalanceIdx);
+ }
+ }
+
+ //! \copydoc Dumux::ImplicitProblem::dirichlet()
+ void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
+ {
+ const GlobalPosition globalPos = vertex.geometry().center();
+
+// initial_(values, globalPos);
+ values = 0.0;
+
+ FluidState fluidState;
+ updateFluidStateForBC_(fluidState);
+
+ const Scalar density =
+ FluidSystem::density(fluidState, phaseIdx);
+
+ values[velocityXIdx] = xVelocity_(globalPos);
+ values[velocityYIdx] = 0.0;
+ values[pressureIdx] = refPressure() +
+ density*this->gravity()[1]*(globalPos[1] - bboxMin_[1]);
+ values[massOrMoleFracIdx] = refMassfrac();
+ values[temperatureIdx] = refTemperature();
+ }
+
+ //! \copydoc Dumux::ImplicitProblem::neumann()
+ void neumann(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const Intersection &is,
+ const int scvIdx,
+ const int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos =
+ fvGeometry.boundaryFace[boundaryFaceIdx].ipGlobal;
+
+ values = 0.;
+
+ FluidState fluidState;
+ updateFluidStateForBC_(fluidState);
+
+ const Scalar density =
+ FluidSystem::density(fluidState, phaseIdx);
+ const Scalar enthalpy =
+ FluidSystem::enthalpy(fluidState, phaseIdx);
+ const Scalar xVelocity = xVelocity_(globalPos);
+
+ if (onLeftBoundary_(globalPos)
+ && globalPos[1] > bboxMin_[1] && globalPos[1] < bboxMax_[1])
+ {
+ values[transportEqIdx] = -xVelocity*density*refMassfrac();
+ values[energyEqIdx] = -xVelocity*density*enthalpy;
+ }
+ }
+
+ /*!
+ * \brief Evaluate the Beavers-Joseph coefficient
+ * at the center of a given intersection
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry
+ * \param is The intersection between element and boundary
+ * \param scvIdx The local subcontrolvolume index
+ * \param boundaryFaceIdx The index of the boundary face
+ *
+ * \return Beavers-Joseph coefficient
+ */
+ Scalar beaversJosephCoeff(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const Intersection &is,
+ const int scvIdx,
+ const int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos =
+ fvGeometry.boundaryFace[boundaryFaceIdx].ipGlobal;
+
+ if (onLowerBoundary_(globalPos))
+ return alphaBJ_;
+ else
+ return 0.0;
+ }
+
+ /*!
+ * \brief Evaluate the intrinsic permeability
+ * at the corner of a given element
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry
+ * \param scvIdx The local subcontrolvolume index
+ *
+ * \return permeability in x-direction
+ */
+ Scalar permeability(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ return spatialParams_.intrinsicPermeability(element,
+ fvGeometry,
+ scvIdx);
+ }
+
+ // \}
+
+ //! \copydoc Dumux::ImplicitProblem::source()
+ void source(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ // ATTENTION: The source term of the mass balance has to be chosen as
+ // div (q_momentum) in the problem file
+ values = Scalar(0);
+ }
+
+ //! \copydoc Dumux::ImplicitProblem::initial()
+ void initial(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ const GlobalPosition &globalPos
+ = element.geometry().corner(scvIdx);
+
+ initial_(values, globalPos);
+ }
+ // \}
+
+ /*!
+ * \brief Determine if we are on a corner of the grid
+ *
+ * \param global Pos The global position
+ *
+ */
+ bool isCornerPoint(const GlobalPosition &globalPos)
+ {
+ if ((onLeftBoundary_(globalPos) && onLowerBoundary_(globalPos)) ||
+ (onLeftBoundary_(globalPos) && onUpperBoundary_(globalPos)) ||
+ (onRightBoundary_(globalPos) && onLowerBoundary_(globalPos)) ||
+ (onRightBoundary_(globalPos) && onUpperBoundary_(globalPos)))
+ return true;
+ else
+ return false;
+ }
+
+ // required in case of mortar coupling
+ // otherwise it should return false
+ /*!
+ * \brief docme
+ *
+ * \param global Pos The global position
+ *
+ */
+ bool isInterfaceCornerPoint(const GlobalPosition &globalPos) const
+ { return false; }
+
+ /*!
+ * \brief Returns the spatial parameters object.
+ */
+ SpatialParams &spatialParams()
+ { return spatialParams_; }
+ const SpatialParams &spatialParams() const
+ { return spatialParams_; }
+
+ const Scalar refTemperature() const
+ { return refTemperature_+ diurnalVariation_(sinusTVar_); }
+
+ const Scalar refMassfrac() const
+ { return refMassfrac_ + diurnalVariation_(sinusXVar_); }
+
+ const Scalar refPressure() const
+ { return refPressure_ + diurnalVariation_(sinusPVar_); }
+
+private:
+ // internal method for the initial condition (reused for the
+ // dirichlet conditions!)
+
+
+ void initial_(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ FluidState fluidState;
+ updateFluidStateForBC_(fluidState);
+
+ const Scalar density =
+ FluidSystem::density(fluidState, phaseIdx);
+
+ values[velocityXIdx] = xVelocity_(globalPos);
+ values[velocityYIdx] = 0.;
+
+ values[pressureIdx] = refPressure()
+ + density*this->gravity()[1]*(globalPos[1] - bboxMin_[1]);
+ values[massOrMoleFracIdx] = refMassfrac();
+ values[temperatureIdx] = refTemperature();
+ }
+
+
+ const Scalar xVelocity_(const GlobalPosition &globalPos) const
+ {
+ const Scalar vmax = vxMax_ + hourlyVariation_(sinusVelVar_);
+// const Scalar relativeHeight = (globalPos[1]-bboxMin_[1])/height_();
+ // linear profile
+// return vmax*relativeHeight + bjSlipVel_; // BJ slip velocity is added as sqrt(Kxx)
+ // parabolic profile
+ return 4*vmax*(globalPos[1] - bboxMin_[1])*(bboxMax_[1] - globalPos[1])
+ / (height_()*height_()) + bjSlipVel_;
+// Scalar tempVel = (globalPos[1] - bboxMin_[1])/height_();
+// return vmax*tempVel*tempVel + bjSlipVel_;
+ // logarithmic profile
+// return 0.1*vmax*log((relativeHeight+1e-3)/1e-3) + bjSlipVel_;
+ }
+
+
+ void updateFluidStateForBC_(FluidState& fluidState) const
+ {
+ fluidState.setTemperature(refTemperature());
+ fluidState.setPressure(phaseIdx, refPressure());
+
+ Scalar massFraction[numComponents];
+ massFraction[transportCompIdx] = refMassfrac();
+ massFraction[phaseCompIdx] = 1 - massFraction[transportCompIdx];
+
+ // calculate average molar mass of the gas phase
+ Scalar M1 = FluidSystem::molarMass(transportCompIdx);
+ Scalar M2 = FluidSystem::molarMass(phaseCompIdx);
+ Scalar X2 = massFraction[phaseCompIdx];
+ Scalar massToMoleDenominator = M2 + X2*(M1 - M2);
+
+ fluidState.setMoleFraction(phaseIdx, transportCompIdx, massFraction[transportCompIdx]*M2/massToMoleDenominator);
+ fluidState.setMoleFraction(phaseIdx, phaseCompIdx, massFraction[phaseCompIdx]*M1/massToMoleDenominator);
+ }
+
+
+ const Scalar diurnalVariation_(const Scalar value) const
+ {
+ const Scalar time = this->timeManager().time() + this->timeManager().timeStepSize();
+ return sin(2*M_PI*time/86400) * value;
+ }
+
+
+ const Scalar hourlyVariation_(const Scalar value) const
+ {
+ const Scalar time = this->timeManager().time();
+ return sin(2*M_PI*time/3600) * value;
+ }
+
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < bboxMin_[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > bboxMax_[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < bboxMin_[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > bboxMax_[1] - eps_; }
+
+ bool onBoundary_(const GlobalPosition &globalPos) const
+ {
+ return (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)
+ || onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos));
+ }
+
+ const Scalar height_() const
+ { return bboxMax_[1] - bboxMin_[1]; }
+
+ // spatial parameters
+ SpatialParams spatialParams_;
+
+ static constexpr Scalar eps_ = 1e-8;
+
+ GlobalPosition bboxMin_;
+ GlobalPosition bboxMax_;
+
+ Scalar refPressure_;
+ Scalar refTemperature_;
+ Scalar refMassfrac_;
+
+ Scalar vxMax_;
+ Scalar bjSlipVel_;
+ Scalar alphaBJ_;
+
+ Scalar sinusVelVar_;
+ Scalar sinusPVar_;
+ Scalar sinusTVar_;
+ Scalar sinusXVar_;
+
+ Scalar xMaterialInterface_;
+ Scalar runUpDistanceX_;
+ Scalar initializationTime_;
+};
+} //end namespace
+
+#endif
diff --git a/test/multidomain/2cnistokes2p2cni/test_2cnistokes2p2cni.cc b/test/multidomain/2cnistokes2p2cni/test_2cnistokes2p2cni.cc
new file mode 100644
index 0000000000..5230301847
--- /dev/null
+++ b/test/multidomain/2cnistokes2p2cni/test_2cnistokes2p2cni.cc
@@ -0,0 +1,279 @@
+// -*- 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 coupled 2cni-Stokes and 2p2cni-Darcy model
+ */
+#include "config.h"
+#include <iostream>
+#include <boost/format.hpp>
+
+#include <dune/common/version.hh>
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3)
+#include <dune/common/parallel/mpihelper.hh>
+#else
+#include <dune/common/mpihelper.hh>
+#endif
+
+#include <dune/common/parametertreeparser.hh>
+#include <dumux/io/interfacemeshcreator.hh>
+
+#include "2cnistokes2p2cniproblem.hh"
+
+
+/*!
+ * \brief Print a usage string for simulations.
+ *
+ * \param progname The name of the executable
+ */
+void printUsage(const char *progName)
+{
+ std::cout << "usage: " << progName
+ << " [--restart restartTime] -parameterFile test_2cnistokes2p2cni.input\n";
+ exit(1);
+}
+
+template <class TypeTag>
+int start_(int argc,
+ char **argv)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, MultiDomainGrid) MDGrid;
+ typedef typename GET_PROP_TYPE(TypeTag, GridCreator) GridCreator;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+// typedef Dune::GridPtr<Grid> GridPointer;
+
+ ////////////////////////////////////////////////////////////
+ // Load the input parameters
+ ////////////////////////////////////////////////////////////
+
+ typedef typename GET_PROP(TypeTag, ParameterTree) ParameterTree;
+ Dune::ParameterTreeParser::readOptions(argc, argv, ParameterTree::tree());
+
+ if (ParameterTree::tree().hasKey("parameterFile") or argc==1)
+ {
+ // read input file, but do not overwrite options specified
+ // on the command line, since the latter have precedence.
+ std::string inputFileName ;
+ if(argc==1) // if there are no arguments given (and there is a file ./<programname>.input) we use it as input file
+ {
+ std::cout<< "\nNo parameter file given. \n"
+ << "Defaulting to '"
+ << argv[0]
+ << ".input' for input file.\n";
+ inputFileName = argv[0];
+ inputFileName += ".input";
+ }
+ else
+ inputFileName = GET_RUNTIME_PARAM(TypeTag, std::string, parameterFile); // otherwise we read from the command line
+
+ std::ifstream parameterFile;
+
+ // check whether the parameter file exists.
+ parameterFile.open(inputFileName.c_str());
+ if (not parameterFile.is_open()){
+ std::cout<< "\n\t -> Could not open file"
+ << inputFileName
+ << ". <- \n\n\n\n";
+ printUsage(argv[0]);
+ return 1;
+ }
+ parameterFile.close();
+
+ Dune::ParameterTreeParser::readINITree(inputFileName,
+ ParameterTree::tree(),
+ /*overwrite=*/false);
+ }
+
+ // initialize MPI, finalize is done automatically on exit
+ static Dune::MPIHelper& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // define the problem dimensions
+ const int dim=2;
+
+ // deal with the restart stuff
+ int argIdx = 1;
+ bool restart = false;
+ double tStart = 0.0;
+ if (argc > 1 && std::string("--restart") == argv[argIdx])
+ {
+ restart = true;
+ ++argIdx;
+
+ std::istringstream(argv[argIdx++]) >> tStart;
+ }
+
+ std::string dgfFileName;
+ Scalar dt, tEnd;
+ Dune::FieldVector<int, dim> nElements;
+ Scalar interfacePos, gradingFactor;
+ int gridRefinement;
+ bool useInterfaceMeshCreator;
+
+ try
+ {
+ dgfFileName = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Grid, File);
+ dt = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, DtInitial);
+ tEnd = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, TEnd);
+ nElements[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsX);
+ if (dim>1) nElements[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsY);
+ if (dim==3) nElements[2] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsZ);
+ interfacePos = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePos);
+ gradingFactor = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, Grading);
+ gridRefinement = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, Refinement);
+ useInterfaceMeshCreator = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, Grid, UseInterfaceMeshCreator);
+ }
+ catch (Dumux::ParameterException &e) {
+ std::cerr << e << ". Abort!\n";
+ exit(1) ;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ exit(1);
+ }
+ std::cout << "Starting with timestep size = " << dt << "s, simulation end = " << tEnd << "s\n";
+
+ if (useInterfaceMeshCreator)
+ {
+ Dumux::InterfaceMeshCreator<Grid> interfaceMeshCreator;
+ GridCreator::gridPtr() = interfaceMeshCreator.create(dgfFileName, nElements, interfacePos, gradingFactor);
+ }
+ else
+ {
+ // try to create a grid (from the given grid file)
+ try { GridCreator::makeGrid(); }
+ catch (...) {
+ std::string usageMessage = "\n\t -> Creation of the grid failed! <- \n\n\n\n";
+// usageMessage += usageTextBlock();
+// usage(argv[0], usageMessage);
+ throw;
+ }
+ }
+
+ if (gridRefinement)
+ GridCreator::grid().globalRefine(gridRefinement);
+
+ if (mpiHelper.size() > 1) {
+ if (!Dune::Capabilities::isParallel<Grid>::v) {
+ std::cerr << "WARNING: THE PROGRAM IS STARTED USING MPI, BUT THE GRID IMPLEMENTATION\n"
+ << " YOU HAVE CHOSEN IS NOT PARALLEL!\n";
+ }
+ GridCreator::loadBalance();
+ }
+
+ // Instantiate the time manager
+ TimeManager timeManager;
+
+ // instantiate coupled problem
+ Dune::shared_ptr<MDGrid> mdGrid_ = Dune::make_shared<MDGrid> (GridCreator::grid());
+
+ Problem problem(*mdGrid_,
+ timeManager);
+
+ Dumux::Parameters::print<TypeTag>();
+
+ // run the simulation
+ timeManager.init(problem,
+ tStart, // initial time
+ dt, // initial time step
+ tEnd, // final time
+ restart);
+
+ // print all properties
+ Dumux::Properties::print<TypeTag>();
+
+ timeManager.run();
+
+ return 0;
+}
+
+/*!
+ * \ingroup Start
+ *
+ * \brief Provides a main function which reads in parameters from the
+ * command line and a parameter file.
+ *
+ * In this function only the differentiation between debugger
+ * or not is made.
+ *
+ * \tparam TypeTag The type tag of the problem which needs to be solved
+ *
+ * \param argc The number of command line arguments of the program
+ * \param argv The contents of the command line arguments of the program
+ * \param usage Callback function for printing the usage message
+ */
+template <class TypeTag>
+int start(int argc,
+ char **argv)
+{
+ try {
+ return start_<TypeTag>(argc, argv);
+ }
+ catch (Dumux::ParameterException &e)
+ {
+ std::cerr << e << ". Abort!\n";
+ printUsage(argv[0]);
+ return 1;
+ }
+ catch (Dune::Exception &e) {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ return 2;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ return 3;
+ }
+}
+
+//
+///*!
+// * \brief Print a usage string for simulations using
+// * Dumux::startFromDGF() as their main() function.
+// *
+// * \param progname The name of the executable
+// */
+//void printUsage(const char *progname)
+//{
+// std::cout << boost::format(
+// "usage: %s [options]\n"
+// "\n"
+// "Options starting with a lowercase letter must be specified somewhere,\n"
+// "options beginning with an uppercase letter do have fallback values.\n"
+// "Some important options include:\n"
+// " -Help Print this usage message and quit\n"
+// " -tEnd time Time [s] of the simulation's end\n"
+// " -dtInitial time Initial time step size [s]\n"
+// " -parameterFile parameter file A INI file with options\n"
+// " -Restart time Restart a previous run from a DRS file\n"
+// " -PrintParams (true|false) Print the parameters used at the end of the simulation\n"
+// " -CellsX n Number of cells in x direction\n"
+// " -CellsY n Number of cells in y direction\n"
+// " -Newton.RelTolerance value Tolerated relative error for the Newton solver\n"
+// )%progname;
+//}
+
+int main(int argc, char** argv)
+{
+ typedef TTAG(TwoCNIStokesTwoPTwoCNIProblem) ProblemTypeTag;
+ return start<ProblemTypeTag>(argc, argv);//, usage);
+}
diff --git a/test/multidomain/2cnistokes2p2cni/test_2cnistokes2p2cni.input b/test/multidomain/2cnistokes2p2cni/test_2cnistokes2p2cni.input
new file mode 100644
index 0000000000..7b736ca114
--- /dev/null
+++ b/test/multidomain/2cnistokes2p2cni/test_2cnistokes2p2cni.input
@@ -0,0 +1,193 @@
+#############################################################
+#Configuration file for test_2cnistokes2p2cni
+#############################################################
+
+#############################################################
+#General conditions: Problem, TimeManager, Vtk, Grid
+#############################################################
+
+#############################################################
+[Problem]
+#############################################################
+#First part of the interface, where it is not coupled;
+#set to zero or negative if not desired
+RunUpDistanceX = 0.0 # 0.06, has to be larger than NoDarcyX !?
+NoDarcyX = 0.0 # 0.05
+
+#############################################################
+[TimeManager]
+#############################################################
+#Initial and maximum time-step size
+DtInitial = 5e-1
+MaxTimeStepSize = 360
+
+#Initialization time without coupling
+#set to zero or negative if not desired
+InitTime = 864
+
+#Simulation end
+TEnd = 0.864e6
+
+#Define the length of an episode (for output)
+EpisodeLength = 43200
+
+#############################################################
+[Vtk]
+#############################################################
+#Names for VTK output
+NameFF = stokes2cni
+NamePM = 2p2cni
+AddVelocity = 1
+
+#############################################################
+[Grid]
+#############################################################
+UseInterfaceMeshCreator = true
+
+#Name of the dgf file (grid)
+File = grids/interfacedomain.dgf
+#File = grids/test_2cnistokes2p2cni.dgf
+Refinement = 0
+
+#Number of elements in x-, y-, z-direction
+CellsX = 30
+CellsY = 62
+CellsZ = 1
+#Grading and refinement of the mesh
+Grading = 1.13
+
+#Extend of entire domain
+XMin = 0.0
+XMax = 0.25
+YMin = 0.0
+YMax = 0.5
+
+#Vertical position of coupling interface
+InterfacePos = 0.25
+
+#############################################################
+[Output]
+#############################################################
+#Frequency of restart file, flux and VTK output
+FreqRestart = 1000 # how often restart files are written out
+FreqOutput = 50 # frequency of VTK output
+FreqMassOutput = 2 # frequency of mass and evaporation rate output (Darcy)
+FreqFluxOutput = 1000 # frequency of detailed flux output
+FreqVaporFluxOutput = 2 # frequency of summarized flux output
+
+#############################################################
+[Stokes]
+#############################################################
+StabilizationAlpha = -1.0
+
+#############################################################
+[FreeFlow]
+#############################################################
+RefPressure = 1e5
+RefTemperature = 298.15
+RefMassfrac = 0.008
+VxMax = 3.5
+BeaversJosephSlipVel = 0.00134
+ExponentMTC = 0.0 # 1./6., Mass transfer coefficient for S^MTC
+UseBoundaryLayerModel = 0 # integer, 0 for no boundary layer model, 1 for Blasius, 2 for turbulent BL, 9 for constant thickness
+BoundaryLayerOffset = 0.00 # For BL model like Blasius, determines a virtual run-up distance for the flow
+ConstThickness = 0.0016 # for a constant BL thickness, use BL model 3
+MassTransferModel = 0 # 0 for none, 1 for power law, 2 for Schluender model
+
+SinusVelVar = 0.0 # hourly velocity variation
+SinusPVar = 0.0 # diurnal pressure variation
+SinusTVar = 0.0 # diurnal temperature variation
+SinusXVar = 0.0 # diurnal variation of the vapor concentration (massfraction)
+
+#############################################################
+[PorousMedium]
+#############################################################
+RefPressurePM = 1e5
+RefTemperaturePM = 298.15
+InitialSw1 = 0.98
+InitialSw2 = 0.98
+CharPoreRadius = 1e-4
+PlotFluidProperties = false
+
+#############################################################
+[SpatialParams]
+#############################################################
+MaterialInterfaceX = -100.0
+AlphaBJ = 1.0
+
+# for homogeneous setups (0 for heterogeneous):
+SoilType = 2
+# threshold of the linearization of the capillary pressure curve
+RegularizationThreshold = 1e-2
+
+# GStat stuff, only required for SoilType=0
+GenerateNewPermeability = true
+GStatControlFileName = gstatControl_2D.txt
+GStatInputFileName = gstatInput.txt
+PermeabilityInputFileName = permeab.dat
+
+### SoilType = 1 ### (Zurich coarse, MUSIS sand)
+[SpatialParams.Coarse]
+Permeability1 = 7.0e-10 #Kozeny-Carman; 1.387e-11 measured?
+Porosity1 = 0.44
+Swr1 = 0.005
+Snr1 = 0.01
+VgAlpha1 = 8.74e-4
+VgN1 = 3.35
+PlotMaterialLaw1 = false
+LambdaSolid1 = 5.26
+
+Pentry1 = 1012
+BCLambda1 = 3.277
+
+### SoilType = 2 ### (Zurich fine)
+[SpatialParams.Medium]
+Permeability2 = 2.65e-10
+Porosity2 = 0.41
+Swr2 = 0.005
+Snr2 = 0.01
+VgAlpha2 = 6.371e-4
+VgN2 = 8.0 # before: 6.9
+PlotMaterialLaw2 = false
+LambdaSolid2 = 5.26
+
+Pentry2 = 1357
+BCLambda2 = 6.960
+
+### SoilType = 4 ### (Zurich fine)
+[SpatialParams.LeverettJ]
+PermeabilityJ = 1.0e-09
+PorosityJ = 0.41
+PlotMaterialLawJ = false
+
+### SoilType = 3 ### (Colorado)
+[SpatialParams.Fine]
+Permeability3 = 1.06e-10
+Porosity3 = 0.334
+Swr3 = 0.028
+Snr3 = 0.01
+VgAlpha3 = 5.81e-4
+VgN3 = 17.8
+PlotMaterialLaw3 = false
+LambdaSolid3 = 5.26
+
+Pentry3 = 1012
+BCLambda3 = 3.277
+
+#############################################################
+[Newton]
+#############################################################
+RelTolerance = 1e-5
+TargetSteps = 8
+MaxSteps = 12
+WriteConvergence = false
+MaxTimeStepDivisions = 20 # should be in group Implicit
+
+#############################################################
+[LinearSolver]
+#############################################################
+ResidualReduction = 1e-9
+Verbosity = 0
+MaxIterations = 200
+#GMResRestart = 100
+#PreconditionerIterations = 2
diff --git a/test/multidomain/Makefile.am b/test/multidomain/Makefile.am
index 08f06370a2..2db4fddc7f 100644
--- a/test/multidomain/Makefile.am
+++ b/test/multidomain/Makefile.am
@@ -1,5 +1,6 @@
SUBDIRS = \
- 2cstokes2p2c
+ 2cstokes2p2c \
+ 2cnistokes2p2cni
EXTRA_DIST = CMakeLists.txt
--
GitLab