From 0e4cd9efbe5000f6cebc47aa0f8f4d1742fc8c03 Mon Sep 17 00:00:00 2001
From: Klaus Mosthaf <klmos@env.dtu.dk>
Date: Fri, 24 Jan 2014 14:41:43 +0000
Subject: [PATCH] added 2cstokes2p2c test (first coupled ff-pm example);
adapted makefiles and configure, added stokesnccouplinglocalresidual
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@12349 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
configure.ac | 2 +
.../stokesnccouplinglocalresidual.hh | 429 ++++++++++
test/Makefile.am | 1 +
.../2cstokes2p2c/2cstokes2p2cproblem.hh | 738 ++++++++++++++++++
.../2cstokes2p2c/2cstokes2p2cspatialparams.hh | 489 ++++++++++++
.../2cstokes2p2c/2p2csubproblem.hh | 459 +++++++++++
test/multidomain/2cstokes2p2c/Makefile.am | 5 +
.../2cstokes2p2c/grids/interfacedomain.dgf | 15 +
.../2cstokes2p2c/grids/test_2cstokes2p2c.dgf | 17 +
.../2cstokes2p2c/stokes2csubproblem.hh | 598 ++++++++++++++
.../2cstokes2p2c/test_2cstokes2p2c.cc | 252 ++++++
.../2cstokes2p2c/test_2cstokes2p2c.input | 170 ++++
test/multidomain/Makefile.am | 7 +
13 files changed, 3182 insertions(+)
create mode 100644 dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh
create mode 100644 test/multidomain/2cstokes2p2c/2cstokes2p2cproblem.hh
create mode 100644 test/multidomain/2cstokes2p2c/2cstokes2p2cspatialparams.hh
create mode 100644 test/multidomain/2cstokes2p2c/2p2csubproblem.hh
create mode 100644 test/multidomain/2cstokes2p2c/Makefile.am
create mode 100644 test/multidomain/2cstokes2p2c/grids/interfacedomain.dgf
create mode 100644 test/multidomain/2cstokes2p2c/grids/test_2cstokes2p2c.dgf
create mode 100644 test/multidomain/2cstokes2p2c/stokes2csubproblem.hh
create mode 100644 test/multidomain/2cstokes2p2c/test_2cstokes2p2c.cc
create mode 100644 test/multidomain/2cstokes2p2c/test_2cstokes2p2c.input
create mode 100644 test/multidomain/Makefile.am
diff --git a/configure.ac b/configure.ac
index 85f211277d..39272660ee 100644
--- a/configure.ac
+++ b/configure.ac
@@ -128,6 +128,8 @@ AC_CONFIG_FILES([dumux.pc
test/material/ncpflash/Makefile
test/material/pengrobinson/Makefile
test/material/tabulation/Makefile
+ test/multidomain/Makefile
+ test/multidomain/2cstokes2p2c/Makefile
test/references/Makefile
tutorial/Makefile
])
diff --git a/dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh b/dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh
new file mode 100644
index 0000000000..27de52eef2
--- /dev/null
+++ b/dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh
@@ -0,0 +1,429 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2012 by Katherina Baber, Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Bernd Flemisch, Andreas Lauser *
+ * Institute for Modelling Hydraulic and Environmental Systems *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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 Element-wise calculation of the Jacobian matrix for problems
+ * using the stokes box model.
+ */
+
+#ifndef DUMUX_STOKESNC_COUPLING_LOCAL_RESIDUAL_HH
+#define DUMUX_STOKESNC_COUPLING_LOCAL_RESIDUAL_HH
+
+#include <dumux/freeflow/stokesnc/stokesncmodel.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup BoxStokesModel
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the Stokes box model.
+ *
+ * This class is also used for the stokes transport
+ * model, which means that it uses static polymorphism.
+ */
+template<class TypeTag>
+class StokesncCouplingLocalResidual : public StokesncLocalResidual<TypeTag>
+{
+protected:
+ typedef StokesncCouplingLocalResidual<TypeTag> ThisType;
+ typedef StokesncLocalResidual<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+ numEq = GET_PROP_VALUE(TypeTag, NumEq),
+ numComponents = FluidSystem::numComponents
+ };
+ enum {
+ //indices of the equations
+ massBalanceIdx = Indices::massBalanceIdx, //!< Index of the mass balance
+
+ 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
+ lastMomentumIdx = Indices::lastMomentumIdx, //!< Index of the last component of the momentum balance
+ transportEqIdx = Indices::transportEqIdx//!< Index of the transport equation
+ };
+ enum {
+ //indices of phase and transported component
+ phaseIdx = Indices::phaseIdx,
+ transportCompIdx = Indices::transportCompIdx
+ };
+ enum {
+ dimXIdx = Indices::dimXIdx, //!< Index for the first component of a vector
+ dimYIdx = Indices::dimYIdx, //!< Index for the second component of a vector
+ dimZIdx = Indices::dimZIdx //!< Index for the third component of a vector
+ };
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::template Codim<dim>::EntityPointer VertexPointer;
+ typedef typename GridView::ctype CoordScalar;
+ typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElements;
+ typedef Dune::GenericReferenceElement<Scalar, dim> ReferenceElement;
+
+ typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<Scalar, dim> DimVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, DofMapper) DofMapper;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementSolutionVector) ElementSolutionVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+
+ typedef typename GridView::Intersection Intersection;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+
+ static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+
+public:
+ // the stokes model needs a modified treatment of the BCs
+ void evalBoundary_()
+ {
+ assert(this->residual_.size() == this->fvGeometry_().numScv);
+ const ReferenceElement &refElement = ReferenceElements::general(this->element_().geometry().type());
+
+ // loop over vertices of the element
+ for (int idx = 0; idx < this->fvGeometry_().numScv; idx++)
+ {
+ // consider only SCVs on the boundary
+ if (this->fvGeometry_().subContVol[idx].inner)
+ continue;
+
+ // important at corners of the grid
+ DimVector momentumResidual(0.0);
+ DimVector averagedNormal(0.0);
+ int numberOfOuterFaces = 0;
+ // evaluate boundary conditions for the intersections of
+ // the current element
+ const BoundaryTypes &bcTypes = this->bcTypes_(idx);
+ IntersectionIterator isIt = this->gridView_().ibegin(this->element_());
+ const IntersectionIterator &endIt = this->gridView_().iend(this->element_());
+ for (; isIt != endIt; ++isIt)
+ {
+ // handle only intersections on the boundary
+ if (!isIt->boundary())
+ continue;
+
+ // assemble the boundary for all vertices of the current face
+ const int faceIdx = isIt->indexInInside();
+ const int numFaceVertices = refElement.size(faceIdx, 1, dim);
+
+ // loop over the single vertices on the current face
+ for (int faceVertIdx = 0; faceVertIdx < numFaceVertices; ++faceVertIdx)
+ {
+ // only evaluate, if we consider the same face vertex as in the outer
+ // loop over the element vertices
+ if (refElement.subEntity(faceIdx, 1, faceVertIdx, dim)
+ != idx)
+ continue;
+
+ const int boundaryFaceIdx = this->fvGeometry_().boundaryFaceIndex(faceIdx, faceVertIdx);
+ const FluxVariables boundaryVars(this->problem_(),
+ this->element_(),
+ this->fvGeometry_(),
+ boundaryFaceIdx,
+ this->curVolVars_(),
+ true);
+
+ // the computed residual of the momentum equations is stored
+ // into momentumResidual for the replacement of the mass balance
+ // in case of Dirichlet conditions for the momentum balance;
+ // the fluxes at the boundary are added in the second step
+ if (this->momentumBalanceDirichlet_(bcTypes))
+ {
+ DimVector muGradVelNormal(0.);
+ const DimVector &boundaryFaceNormal =
+ boundaryVars.face().normal;
+
+ boundaryVars.velocityGrad().umv(boundaryFaceNormal, muGradVelNormal);
+ muGradVelNormal *= boundaryVars.dynamicViscosity();
+
+ for (int i=0; i < this->residual_.size(); i++)
+ Valgrind::CheckDefined(this->residual_[i]);
+ for (int dimIdx=0; dimIdx < dim; ++dimIdx)
+ momentumResidual[dimIdx] = this->residual_[idx][momentumXIdx+dimIdx];
+
+ //Sign is right!!!: boundary flux: -mu grad v n
+ //but to compensate outernormal -> residual - (-mu grad v n)
+ momentumResidual += muGradVelNormal;
+ averagedNormal += boundaryFaceNormal;
+ }
+
+ // evaluate fluxes at a single boundary segment
+ asImp_()->evalNeumannSegment_(isIt, idx, boundaryFaceIdx, boundaryVars);
+ asImp_()->evalOutflowSegment_(isIt, idx, boundaryFaceIdx, boundaryVars);
+
+ //for the corner points, the boundary flux across the vertical non-coupling boundary faces
+ //has to be calculated to fulfill the mass balance
+ //convert suddomain intersection into multidomain intersection and check whether it is an outer boundary
+ if(!GridView::Grid::multiDomainIntersection(*isIt).neighbor()
+ && this->boundaryHasMortarCoupling_(this->bcTypes_(idx)))
+ {
+ const GlobalPosition& globalPos = this->fvGeometry_().subContVol[idx].global;
+ //problem specific function, in problem orientation of interface is known
+ if(this->problem_().isInterfaceCornerPoint(globalPos))
+ {
+ PrimaryVariables priVars(0.0);
+// DimVector faceCoord = this->fvGeometry_().boundaryFace[boundaryFaceIdx].ipGlobal;
+// std::cout<<faceCoord<<std::endl;
+
+ const int numVertices = refElement.size(dim);
+ bool evalBoundaryFlux = false;
+ for(int equationIdx = 0; equationIdx < numEq; ++equationIdx)
+ {
+ for(int i= 0; i < numVertices; i++)
+ {
+ //if vertex is on boundary and not the coupling vertex: check whether an outflow condition is set
+ if(this->model_().onBoundary(this->element_(), i) && i!=idx)
+ if (!this->bcTypes_(i).isOutflow(equationIdx))
+ evalBoundaryFlux = true;
+ }
+
+ //calculate the actual boundary fluxes and add to residual (only for momentum and transport equation, mass balance already has outflow)
+ asImp_()->computeFlux(priVars, boundaryFaceIdx, true/*on boundary*/);
+ if(evalBoundaryFlux)
+ this->residual_[idx][equationIdx] += priVars[equationIdx];
+ }
+ }
+ }
+ // Beavers-Joseph condition at the coupling boundary/interface
+ if(boundaryHasCoupling_(bcTypes) || boundaryHasMortarCoupling_(bcTypes))
+ {
+ evalBeaversJoseph_(isIt, idx, boundaryFaceIdx, boundaryVars);
+ asImp_()->evalCouplingVertex_(isIt, idx, boundaryFaceIdx, boundaryVars);
+ }
+
+ // count the number of outer faces to determine, if we are on
+ // a corner point and if an interpolation should be done
+ numberOfOuterFaces++;
+ } // end loop over face vertices
+ } // end loop over intersections
+
+ // replace defect at the corner points of the grid
+ // by the interpolation of the primary variables
+ if(!bcTypes.isDirichlet(massBalanceIdx))
+ {
+ if (this->momentumBalanceDirichlet_(bcTypes))
+ this->replaceMassbalanceResidual_(momentumResidual, averagedNormal, idx);
+ else // de-stabilize (remove alpha*grad p - alpha div f
+ // from computeFlux on the boundary)
+ this->removeStabilizationAtBoundary_(idx);
+ }
+ if (numberOfOuterFaces == 2)
+ this->interpolateCornerPoints_(bcTypes, idx);
+ } // end loop over element vertices
+
+ // evaluate the dirichlet conditions of the element
+ if (this->bcTypes_().hasDirichlet())
+ asImp_()->evalDirichlet_();
+ }
+
+ void evalBoundaryPDELab_()
+ {
+ // loop over vertices of the element
+ for (int idx = 0; idx < this->fvGeometry_().numScv; idx++)
+ {
+ // consider only SCVs on the boundary
+ if (this->fvGeometry_().subContVol[idx].inner)
+ continue;
+
+ this->removeStabilizationAtBoundary_(idx);
+ } // end loop over vertices
+ }
+
+protected:
+ /*!
+ * \brief Evaluate one part of the Dirichlet-like coupling conditions for a single
+ * sub-control volume face; rest is done in the local coupling operator
+ */
+ void evalCouplingVertex_(const IntersectionIterator &isIt,
+ const int scvIdx,
+ const int boundaryFaceIdx,
+ const FluxVariables &boundaryVars)
+ {
+ const VolumeVariables &volVars = this->curVolVars_()[scvIdx];
+ const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
+
+ // TODO: beaversJosephCoeff is used to determine, if we are on a coupling face
+ // this is important at the corners. However, a better way should be found.
+
+ // check if BJ-coeff is not zero to decide, if coupling condition
+ // for the momentum balance (Dirichlet vor v.n) has to be set;
+ // may be important at corners
+ const Scalar beaversJosephCoeff = this->problem_().beaversJosephCoeff(this->element_(),
+ this->fvGeometry_(),
+ *isIt,
+ scvIdx,
+ boundaryFaceIdx);
+ // set velocity normal to the interface
+ if (bcTypes.isCouplingInflow(momentumYIdx) && beaversJosephCoeff)
+ this->residual_[scvIdx][momentumYIdx] =
+ volVars.velocity() *
+ boundaryVars.face().normal /
+ boundaryVars.face().normal.two_norm();
+ Valgrind::CheckDefined(this->residual_[scvIdx][momentumYIdx]);
+
+ // add pressure correction - required for pressure coupling,
+ // if p.n comes from the pm
+ if ((bcTypes.isCouplingOutflow(momentumYIdx) && beaversJosephCoeff) || bcTypes.isMortarCoupling(momentumYIdx))
+ {
+ DimVector pressureCorrection(isIt->centerUnitOuterNormal());
+ pressureCorrection *= volVars.pressure(); // TODO: 3D
+ this->residual_[scvIdx][momentumYIdx] += pressureCorrection[momentumYIdx]*
+ boundaryVars.face().area;
+ }
+
+ // set mole or mass fraction for the transported components
+ for (int compIdx = 0; compIdx < numComponents; compIdx++)
+ {
+ int eqIdx = dim + compIdx;
+ if (eqIdx != massBalanceIdx) {
+ if (bcTypes.isCouplingOutflow(eqIdx))
+ {
+ if(useMoles)
+ this->residual_[scvIdx][eqIdx] = volVars.fluidState().moleFraction(phaseIdx, compIdx);
+ else
+ this->residual_[scvIdx][eqIdx] = volVars.fluidState().massFraction(phaseIdx, compIdx);
+ Valgrind::CheckDefined(this->residual_[scvIdx][eqIdx]);
+ }
+ }
+ }
+ }
+
+ void evalBeaversJoseph_(const IntersectionIterator &isIt,
+ const int scvIdx,
+ const int boundaryFaceIdx,
+ const FluxVariables &boundaryVars) //TODO: required
+ {
+ const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
+
+ Scalar beaversJosephCoeff = this->problem_().beaversJosephCoeff(this->element_(),
+ this->fvGeometry_(),
+ *isIt,
+ scvIdx,
+ boundaryFaceIdx);
+
+ // only enter here, if we are on a coupling boundary (see top)
+ // and the BJ coefficient is not zero
+ if (beaversJosephCoeff)
+ {
+ const Scalar Kxx = this->problem_().permeability(this->element_(),
+ this->fvGeometry_(),
+ scvIdx);
+ beaversJosephCoeff /= sqrt(Kxx);
+ const DimVector& elementUnitNormal = isIt->centerUnitOuterNormal();
+
+ // implementation as NEUMANN condition /////////////////////////////////////////////
+ // (v.n)n
+ if (bcTypes.isCouplingOutflow(momentumXIdx))
+ {
+ const Scalar normalComp = boundaryVars.velocity()*elementUnitNormal;
+ DimVector normalV = elementUnitNormal;
+ normalV *= normalComp; // v*n*n
+
+ // v_tau = v - (v.n)n
+ const DimVector tangentialV = boundaryVars.velocity() - normalV;
+ const Scalar boundaryFaceArea = boundaryVars.face().area;
+
+ for (int dimIdx=0; dimIdx < dim; ++dimIdx)
+ this->residual_[scvIdx][dimIdx] += beaversJosephCoeff *
+ boundaryFaceArea *
+ tangentialV[dimIdx] *
+ boundaryVars.dynamicViscosity();
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //normal component has only to be set if no coupling conditions are defined
+ //set NEUMANN flux (set equal to pressure in problem)
+// PrimaryVariables priVars(0.0);
+// this->problem_().neumann(priVars, this->element_(), this->fvGeometry_(),
+// *isIt, scvIdx, boundaryFaceIdx);
+// for (int dimIdx=0; dimIdx < dim; ++dimIdx)
+// this->residual_[scvIdx][dimIdx] += priVars[dimIdx]*
+// boundaryFaceArea;
+ }
+ if (bcTypes.isCouplingInflow(momentumXIdx)) //TODO: 3D
+ {
+ ///////////////////////////////////////////////////////////////////////////////////////////
+ //IMPLEMENTATION AS DIRICHLET CONDITION
+ //tangential componment: vx = sqrt K /alpha * (grad v n(unity))t
+ DimVector tangentialVelGrad(0);
+ boundaryVars.velocityGrad().umv(elementUnitNormal, tangentialVelGrad);
+ tangentialVelGrad /= -beaversJosephCoeff; // was - before
+ this->residual_[scvIdx][momentumXIdx] =
+ tangentialVelGrad[momentumXIdx] - this->curPriVars_(scvIdx)[momentumXIdx];
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //for testing Beavers and Joseph without adjacent porous medium set vy to 0
+// Scalar normalVel(0);
+// this->residual_[scvIdx][momentumYIdx] = this->curPrimaryVars_(scvIdx)[momentumYIdx] - normalVel;
+ ////////////////////////////////////////////////////////////////////////////////////
+ }
+ }
+ }
+
+ // return true, if at least one equation on the boundary has a coupling condition
+ bool boundaryHasCoupling_(const BoundaryTypes& bcTypes) const
+ {
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ if (bcTypes.isCouplingInflow(eqIdx) || bcTypes.isCouplingOutflow(eqIdx))
+ return true;
+ return false;
+ }
+
+ // return true, if at least one equation on the boundary has a mortar coupling condition
+ bool boundaryHasMortarCoupling_(const BoundaryTypes& bcTypes) const
+ {
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ if (bcTypes.isMortarCoupling(eqIdx))
+ return true;
+ return false;
+ }
+
+
+private:
+ Implementation *asImp_()
+ { return static_cast<Implementation *>(this); }
+ const Implementation *asImp_() const
+ { return static_cast<const Implementation *>(this); }
+
+};
+
+}
+
+#endif
diff --git a/test/Makefile.am b/test/Makefile.am
index 06b1fd9513..58568807e9 100644
--- a/test/Makefile.am
+++ b/test/Makefile.am
@@ -4,6 +4,7 @@ SUBDIRS = common \
geomechanics \
implicit \
material \
+ multidomain \
references
EXTRA_DIST = CMakeLists.txt
diff --git a/test/multidomain/2cstokes2p2c/2cstokes2p2cproblem.hh b/test/multidomain/2cstokes2p2c/2cstokes2p2cproblem.hh
new file mode 100644
index 0000000000..728995cd9b
--- /dev/null
+++ b/test/multidomain/2cstokes2p2c/2cstokes2p2cproblem.hh
@@ -0,0 +1,738 @@
+// -*- 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_2CSTOKES2P2CPROBLEM_HH
+#define DUMUX_2CSTOKES2P2CPROBLEM_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/2cstokes2p2clocaloperator.hh>
+#include <dumux/multidomain/2cstokes2p2c/2cstokes2p2cnewtoncontroller.hh>
+#include <dumux/material/fluidsystems/h2oairfluidsystem.hh>
+#include <dumux/linear/pardisobackend.hh>
+
+#include "2cstokes2p2cspatialparams.hh"
+#include "stokes2csubproblem.hh"
+#include "2p2csubproblem.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class TwoCStokesTwoPTwoCProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(TwoCStokesTwoPTwoCProblem, INHERITS_FROM(MultiDomain));
+
+// Set the grid type
+SET_PROP(TwoCStokesTwoPTwoCProblem, 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(TwoCStokesTwoPTwoCProblem, GridView)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, MultiDomainGrid) MDGrid;
+ public:
+ typedef typename MDGrid::LeafGridView type;
+};
+
+// Set the global problem
+SET_TYPE_PROP(TwoCStokesTwoPTwoCProblem, Problem, Dumux::TwoCStokesTwoPTwoCProblem<TypeTag>);
+
+// Set the two sub-problems of the global problem
+SET_TYPE_PROP(TwoCStokesTwoPTwoCProblem, SubDomain1TypeTag, TTAG(Stokes2cSubProblem));
+SET_TYPE_PROP(TwoCStokesTwoPTwoCProblem, SubDomain2TypeTag, TTAG(TwoPTwoCSubProblem));
+
+// Set the global problem in the context of the two sub-problems
+SET_TYPE_PROP(Stokes2cSubProblem, MultiDomainTypeTag, TTAG(TwoCStokesTwoPTwoCProblem));
+SET_TYPE_PROP(TwoPTwoCSubProblem, MultiDomainTypeTag, TTAG(TwoCStokesTwoPTwoCProblem));
+
+// Set the other sub-problem for each of the two sub-problems
+SET_TYPE_PROP(Stokes2cSubProblem, OtherSubDomainTypeTag, TTAG(TwoPTwoCSubProblem));
+SET_TYPE_PROP(TwoPTwoCSubProblem, OtherSubDomainTypeTag, TTAG(Stokes2cSubProblem));
+
+SET_PROP(Stokes2cSubProblem, SpatialParams)
+{ typedef Dumux::TwoCStokesTwoPTwoCSpatialParams<TypeTag> type; };
+SET_PROP(TwoPTwoCSubProblem, SpatialParams)
+{ typedef Dumux::TwoCStokesTwoPTwoCSpatialParams<TypeTag> type; };
+
+// Set the fluid system
+SET_PROP(TwoCStokesTwoPTwoCProblem, FluidSystem)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef Dumux::FluidSystems::H2OAir<Scalar,
+ Dumux::H2O<Scalar>,
+ /*useComplexrelations=*/false> type;
+};
+
+SET_TYPE_PROP(TwoCStokesTwoPTwoCProblem, JacobianAssembler, Dumux::MultiDomainAssembler<TypeTag>);
+
+// Set the local coupling operator
+SET_TYPE_PROP(TwoCStokesTwoPTwoCProblem, MultiDomainCouplingLocalOperator, Dumux::TwoCStokesTwoPTwoCLocalOperator<TypeTag>);
+
+SET_PROP(TwoCStokesTwoPTwoCProblem, SolutionVector)
+{
+ private:
+ typedef typename GET_PROP_TYPE(TypeTag, MultiDomainGridOperator) MDGridOperator;
+ public:
+ typedef typename MDGridOperator::Traits::Domain type;
+};
+
+// newton controller
+SET_TYPE_PROP(TwoCStokesTwoPTwoCProblem, NewtonController, Dumux::TwoCStokesTwoPTwoCNewtonController<TypeTag>);
+
+// if you do not have PARDISO, the SuperLU solver is used:
+#ifdef HAVE_PARDISO
+SET_TYPE_PROP(TwoCStokesTwoPTwoCProblem, LinearSolver, PardisoBackend<TypeTag>);
+#else
+SET_TYPE_PROP(TwoCStokesTwoPTwoCProblem, 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(TwoCStokesTwoPTwoCProblem, NumEq, 1);
+
+// Write the solutions of individual newton iterations?
+SET_BOOL_PROP(TwoCStokesTwoPTwoCProblem, NewtonWriteConvergence, false);
+}
+
+/*!
+ * \brief The problem class for the coupling of an isothermal two-component Stokes
+ * and an isothermal two-phase two-component Darcy model.
+ *
+ * The problem class for the coupling of an isothermal two-component Stokes (stokes2c)
+ * and an isothermal two-phase two-component Darcy model (2p2c).
+ * It uses the 2p2cCoupling model and the Stokes2ccoupling 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,
+ * 2p2csubproblem.hh and stokes2csubproblem.hh, which are accessible via the coupled problem.
+ */
+template <class TypeTag = TTAG(TwoCStokesTwoPTwoCProblem) >
+class TwoCStokesTwoPTwoCProblem : public MultiDomainProblem<TypeTag>
+{
+ template<int dim>
+ struct VertexLayout
+ {
+ bool contains(Dune::GeometryType geomtype)
+ { return geomtype.dim() == 0; }
+ };
+
+ typedef TwoCStokesTwoPTwoCProblem<TypeTag> ThisType;
+ typedef MultiDomainProblem<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, NewtonController) NewtonController;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+ typedef typename GET_PROP_TYPE(TypeTag, SubDomain1TypeTag) Stokes2cTypeTag;
+ typedef typename GET_PROP_TYPE(TypeTag, SubDomain2TypeTag) TwoPTwoCTypeTag;
+
+ typedef typename GET_PROP_TYPE(Stokes2cTypeTag, Problem) Stokes2cSubProblem;
+ typedef typename GET_PROP_TYPE(TwoPTwoCTypeTag, Problem) TwoPTwoCSubProblem;
+
+ typedef typename GET_PROP_TYPE(Stokes2cTypeTag, ElementSolutionVector) ElementSolutionVector1;
+ typedef typename GET_PROP_TYPE(TwoPTwoCTypeTag, ElementSolutionVector) ElementSolutionVector2;
+
+ typedef typename GET_PROP_TYPE(Stokes2cTypeTag, ElementVolumeVariables) ElementVolumeVariables1;
+ typedef typename GET_PROP_TYPE(TwoPTwoCTypeTag, ElementVolumeVariables) ElementVolumeVariables2;
+
+ typedef typename GET_PROP_TYPE(Stokes2cTypeTag, FluxVariables) BoundaryVariables1;
+
+ typedef typename GET_PROP_TYPE(Stokes2cTypeTag, FVElementGeometry) FVElementGeometry1;
+ typedef typename GET_PROP_TYPE(TwoPTwoCTypeTag, FVElementGeometry) FVElementGeometry2;
+
+ typedef typename GET_PROP_TYPE(TypeTag, MultiDomainGrid) MDGrid;
+ typedef typename MDGrid::LeafGridView MDGridView;
+ enum { dim = MDGridView::dimension };
+
+ typedef typename MDGrid::Traits::template Codim<0>::EntityPointer MDElementPointer;
+ typedef typename MDGrid::Traits::template Codim<0>::Entity MDElement;
+ typedef typename MDGrid::SubDomainGrid SDGrid;
+ typedef typename SDGrid::template Codim<0>::EntityPointer SDElementPointer;
+
+ typedef typename GET_PROP_TYPE(Stokes2cTypeTag, GridView) Stokes2cGridView;
+ typedef typename GET_PROP_TYPE(TwoPTwoCTypeTag, GridView) TwoPTwoCGridView;
+ typedef typename Stokes2cGridView::template Codim<0>::Entity SDElement1;
+ typedef typename TwoPTwoCGridView::template Codim<0>::Entity SDElement2;
+
+ typedef typename MDGrid::template Codim<0>::LeafIterator ElementIterator;
+ typedef typename MDGrid::LeafSubDomainInterfaceIterator SDInterfaceIterator;
+
+ typedef Dune::FieldVector<Scalar, dim> GlobalPosition;
+ typedef Dune::FieldVector<Scalar, dim> FieldVector;
+
+ typedef typename GET_PROP_TYPE(Stokes2cTypeTag, Indices) Stokes2cIndices;
+ typedef typename GET_PROP_TYPE(TwoPTwoCTypeTag, Indices) TwoPTwoCIndices;
+
+ typedef typename GET_PROP_TYPE(Stokes2cTypeTag, FluidSystem) FluidSystem;
+
+ enum { numEq1 = GET_PROP_VALUE(Stokes2cTypeTag, NumEq) };
+ enum { // indices in the Stokes domain
+ momentumXIdx1 = Stokes2cIndices::momentumXIdx, //!< Index of the x-component of the momentum balance
+ momentumYIdx1 = Stokes2cIndices::momentumYIdx, //!< Index of the y-component of the momentum balance
+ momentumZIdx1 = Stokes2cIndices::momentumZIdx, //!< Index of the z-component of the momentum balance
+ lastMomentumIdx1 = Stokes2cIndices::lastMomentumIdx, //!< Index of the last component of the momentum balance
+ massBalanceIdx1 = Stokes2cIndices::massBalanceIdx, //!< Index of the mass balance
+ transportEqIdx1 = Stokes2cIndices::transportEqIdx, //!< Index of the transport equation
+ };
+
+ enum { numEq2 = GET_PROP_VALUE(TwoPTwoCTypeTag, NumEq) };
+ enum { // indices in the Darcy domain
+ massBalanceIdx2 = TwoPTwoCIndices::pressureIdx,
+ switchIdx2 = TwoPTwoCIndices::switchIdx,
+ contiTotalMassIdx2 = TwoPTwoCIndices::contiNEqIdx,
+ contiWEqIdx2 = TwoPTwoCIndices::contiWEqIdx,
+
+ wCompIdx2 = TwoPTwoCIndices::wCompIdx,
+ nCompIdx2 = TwoPTwoCIndices::nCompIdx,
+
+ wPhaseIdx2 = TwoPTwoCIndices::wPhaseIdx,
+ nPhaseIdx2 = TwoPTwoCIndices::nPhaseIdx
+
+ };
+ enum { phaseIdx = GET_PROP_VALUE(Stokes2cTypeTag, PhaseIdx) };
+ enum { transportCompIdx1 = Stokes2cIndices::transportCompIdx };
+
+
+public:
+ /*!
+ * \brief docme
+ *
+ * \param hostGrid docme
+ * \param timeManager The time manager
+ *
+ */
+ TwoCStokesTwoPTwoCProblem(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);
+
+ // 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);
+ }
+
+ stokes2c_ = this->sdID1();
+ twoPtwoC_ = this->sdID2();
+
+ initializeGrid();
+
+ // initialize the tables of the fluid system
+ FluidSystem::init(/*tempMin=*/273.15, /*tempMax=*/373.15, /*numTemp=*/200,
+ /*pMin=*/1e4, /*pMax=*/2e5, /*numP=*/200);
+
+ if (initializationTime_ > 0.0)
+ this->timeManager().startNextEpisode(initializationTime_);
+ else
+ this->timeManager().startNextEpisode(episodeLength_);
+ }
+
+ ~TwoCStokesTwoPTwoCProblem()
+ {
+ 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; flux2; wPhaseFlux2; nPhaseFlux2\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(stokes2c_,*elementIt);
+ else
+ if(globalPos[0] > noDarcyX_)
+ mdGrid.addToSubDomain(twoPtwoC_,*elementIt);
+ }
+ mdGrid.preUpdateSubDomains();
+ mdGrid.updateSubDomains();
+ mdGrid.postUpdateSubDomains();
+
+ gridinfo(this->sdGrid1());
+ gridinfo(this->sdGrid2());
+
+ this->sdProblem1().spatialParams().loadIntrinsicPermeability(this->sdGridView1());
+ this->sdProblem2().spatialParams().loadIntrinsicPermeability(this->sdGridView2());
+ }
+
+ //! \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_); }
+
+ /*
+ * \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.;
+
+ // count number of elements to determine number of interface nodes
+ int numElements = 0;
+ const SDInterfaceIterator endIfIt = mdGrid.leafSubDomainInterfaceEnd(stokes2c_, twoPtwoC_);
+ for (SDInterfaceIterator ifIt =
+ mdGrid.leafSubDomainInterfaceBegin(stokes2c_, twoPtwoC_); 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 faces on the interface
+ for (SDInterfaceIterator ifIt =
+ mdGrid.leafSubDomainInterfaceBegin(stokes2c_, twoPtwoC_); 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(stokes2c_, 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];
+ }
+ }
+ } // 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 "
+ << "count1 "
+ << 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].count << " "
+ << std::endl;
+ }
+ outfile.close();
+ }
+ if (shouldWriteVaporFlux())
+ fluxFile_ << this->timeManager().time() + this->timeManager().timeStepSize() << "; "
+ << sumVaporFluxes << "; " << advectiveVaporFlux << "; " << diffusiveVaporFlux << "; "
+ << advectiveVaporFlux-diffusiveVaporFlux << "; ";
+ }
+
+ /*
+ * \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 sumWaterFluxInGasPhase = 0.;
+
+ // count number of elements to determine number of interface nodes
+ int numElements = 0;
+ const SDInterfaceIterator endIfIt = mdGrid.leafSubDomainInterfaceEnd(stokes2c_, twoPtwoC_);
+ for (SDInterfaceIterator ifIt =
+ mdGrid.leafSubDomainInterfaceBegin(stokes2c_, twoPtwoC_); 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(stokes2c_, twoPtwoC_); 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(twoPtwoC_,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);
+ }
+ }
+ }
+
+ 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 "
+ << 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
+ << std::endl;
+ }
+ outfile.close();
+ }
+ if (shouldWriteVaporFlux()){
+ Scalar sumWaterFluxInLiquidPhase = sumVaporFluxes - sumWaterFluxInGasPhase;
+ fluxFile_ << sumVaporFluxes << "; "
+ << sumWaterFluxInLiquidPhase << "; "
+ << sumWaterFluxInGasPhase << std::endl;
+ }
+ }
+
+ /*!
+ * \brief docme
+ *
+ * \param elemVolVars1 docme
+ * \param boundaryVars1 docme
+ * \param vertInElem1 docme
+ *
+ */
+ 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;
+ }
+
+ /*!
+ * \brief docme
+ *
+ * \param elemVolVars1 docme
+ * \param boundaryVars1 docme
+ * \param vertInElem1 docme
+ *
+ */
+ 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;
+
+ }
+
+ Stokes2cSubProblem& stokes2cProblem()
+ { return this->sdProblem1(); }
+ const Stokes2cSubProblem& stokes2cProblem() const
+ { return this->sdProblem1(); }
+
+ TwoPTwoCSubProblem& twoPtwoCProblem()
+ { return this->sdProblem2(); }
+ const TwoPTwoCSubProblem& twoPtwoCProblem() const
+ { return this->sdProblem2(); }
+
+private:
+ typename MDGrid::SubDomainType stokes2c_;
+ typename MDGrid::SubDomainType twoPtwoC_;
+
+ unsigned counter_;
+ unsigned freqRestart_;
+ unsigned freqOutput_;
+ unsigned freqFluxOutput_;
+ unsigned freqVaporFluxOutput_;
+
+ Scalar interface_;
+ Scalar noDarcyX_;
+ Scalar episodeLength_;
+ Scalar initializationTime_;
+
+ 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/2cstokes2p2c/2cstokes2p2cspatialparams.hh b/test/multidomain/2cstokes2p2c/2cstokes2p2cspatialparams.hh
new file mode 100644
index 0000000000..2707b898e9
--- /dev/null
+++ b/test/multidomain/2cstokes2p2c/2cstokes2p2cspatialparams.hh
@@ -0,0 +1,489 @@
+// -*- 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_TWOCSTOKES_2P2C_SPATIALPARAMS_HH
+#define DUMUX_TWOCSTOKES_2P2C_SPATIALPARAMS_HH
+
+#include <dune/grid/io/file/vtk/common.hh>
+
+//#include <dumux/material/spatialparameters/gstatrandompermeability.hh>
+#include <dumux/material/spatialparams/implicitspatialparams.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
+//#include <dumux/material/fluidmatrixinteractions/2p/linearizedregvangenuchtenparams.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 TwoCStokesTwoPTwoCSpatialParams;
+
+namespace Properties
+{
+// The spatial parameters TypeTag
+NEW_TYPE_TAG(TwoCStokesTwoPTwoCSpatialParams);
+
+// Set the spatial parameters
+SET_TYPE_PROP(TwoCStokesTwoPTwoCSpatialParams, SpatialParams,
+ TwoCStokesTwoPTwoCSpatialParams<TypeTag>);
+
+// Set the material Law
+SET_PROP(TwoCStokesTwoPTwoCSpatialParams, MaterialLaw)
+{
+private:
+ // define the material law which is parameterized by effective
+ // saturations
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+// typedef RegularizedBrooksCorey<Scalar> EffMaterialLaw;
+// typedef RegularizedVanGenuchten<Scalar, LinearizedRegVanGenuchtenParams<Scalar, TypeTag> > EffMaterialLaw;
+ typedef RegularizedVanGenuchten<Scalar> EffMaterialLaw;
+public:
+ // define the material law parameterized by absolute saturations
+ typedef EffToAbsLaw<EffMaterialLaw> type;
+};
+}
+
+
+/** \todo Please doc me! */
+template<class TypeTag>
+class TwoCStokesTwoPTwoCSpatialParams : 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 typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ enum {
+ wPhaseIdx = FluidSystem::wPhaseIdx,
+ nPhaseIdx = FluidSystem::nPhaseIdx
+ };
+
+ 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 Spatial parameters for the coupled 2cstokes and 2p2c test
+ *
+ * \param gridView The GridView which is used by the problem
+ */
+ TwoCStokesTwoPTwoCSpatialParams(const GridView& gridView)
+ : ParentType(gridView),
+// permeability_(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);
+ }
+ }
+
+ ~TwoCStokesTwoPTwoCSpatialParams()
+ {}
+
+ /*!
+ * \brief The intrinsic permeability of the porous medium.
+ *
+ * \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
+ {
+// // heterogeneous parameter field computed with GSTAT
+// if (soilType_ == 0)
+// return permeability_[indexSet_.index(
+// *(element.template subEntity<dim> (scvIdx)))];
+// // soil can be chosen by setting the parameter soiltype accordingly
+// else
+// {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ if (checkSoilType(globalPos) == 1)
+ return coarsePermeability_;
+ if (checkSoilType(globalPos) == 3)
+ 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) == 3)
+ return finePorosity_;
+ else
+ return mediumPorosity_;
+ }
+
+
+ //! \copydoc Dumux::ImplicitSpatialParams::materialLawParams()
+ const MaterialLawParams& materialLawParams(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+// if (soilType_ == 0)
+// return materialLawParams_[indexSet_.index(
+// *(element.template subEntity<dim> (scvIdx)))];
+// else
+// {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+ if (checkSoilType(globalPos)==1)
+ return coarseParams_;
+ if (checkSoilType(globalPos)==3)
+ return fineParams_;
+ 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
+ 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) == 3)
+ return fineLambdaSolid_;
+ else
+ return mediumLambdaSolid_;
+ }
+
+ /*!
+ * \brief docme
+ *
+ * \param global Pos The global position
+ *
+ */
+ const unsigned checkSoilType(const GlobalPosition &globalPos) const
+ {
+ // one soil, which can be chosen as runtime parameter:
+ // 1: coarse
+ // 2: medium
+ // 3: fine
+ return soilType_;
+ }
+
+ /*!
+ * \brief This method allows the generation of a statistical field for the soil parameters
+ * using GStat
+ *
+ * \param gridView The GridView which is used by the problem
+ *
+ */
+ void loadIntrinsicPermeability(const GridView& gridView)
+ {
+// // only load random field, if soilType is set to 0
+// if (soilType_ != 0)
+// return;
+//
+// const unsigned size = gridView.size(dim);
+// permeability_.resize(size, 0.0);
+// materialLawParams_.resize(size);
+//
+// bool create = true;
+// if (ParameterTree::tree().hasKey("SpatialParams.GenerateNewPermeability"))
+// create = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool , SpatialParams, GenerateNewPermeability);
+//
+// std::string gStatControlFileName("gstatControl_2D.txt");
+// if (ParameterTree::tree().hasKey("SpatialParams.GStatControlFileName"))
+// gStatControlFileName = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, SpatialParams,
+// GStatControlFileName);
+//
+//
+// std::string gStatInputFileName("gstatInput.txt");
+// if (ParameterTree::tree().hasKey("SpatialParams.GStatInputFileName"))
+// gStatInputFileName = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, SpatialParams,
+// GStatInputFileName);
+//
+// std::string permeabilityFileName("permeab.dat");
+// if (ParameterTree::tree().hasKey("SpatialParams.PermeabilityInputFileName"))
+// permeabilityFileName = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, SpatialParams,
+// PermeabilityInputFileName);
+//
+// // create random permeability object
+// GstatRandomPermeability<GridView, Scalar> randomPermeability(gridView,
+// create,
+// permeabilityFileName.c_str(),
+// gStatControlFileName.c_str(),
+// gStatInputFileName.c_str());
+//
+// Scalar totalVolume = 0;
+// Scalar meanPermeability = 0;
+//
+//// for (unsigned n=0; n<numVertices; ++n)
+//// {
+//// //TODO the following can be done in a loop
+//// }
+// // iterate over elements
+// ElementIterator eItEnd = gridView.template end<0> ();
+// for (ElementIterator eIt = gridView.template begin<0> (); eIt
+// != eItEnd; ++eIt)
+// {
+// Scalar perm = randomPermeability.K(*eIt)[0][0];
+// permeability_[indexSet_.index(*(eIt->template subEntity<dim> (3)))]
+// = perm;
+//
+// Scalar volume = eIt->geometry().volume();
+// totalVolume += volume;
+//
+// meanPermeability += volume / perm;
+// }
+//
+// for (ElementIterator eIt = gridView.template begin<0> (); eIt
+// != eItEnd; ++eIt)
+// {
+// Scalar perm = randomPermeability.K(*eIt)[0][0];
+// permeability_[indexSet_.index(*(eIt->template subEntity<dim> (1)))]
+// = perm;
+//
+// Scalar volume = eIt->geometry().volume();
+// totalVolume += volume;
+//
+// meanPermeability += volume / perm;
+// }
+//
+// for (ElementIterator eIt = gridView.template begin<0> (); eIt
+// != eItEnd; ++eIt)
+// {
+// Scalar perm = randomPermeability.K(*eIt)[0][0];
+// permeability_[indexSet_.index(*(eIt->template subEntity<dim> (2)))]
+// = perm;
+//
+// Scalar volume = eIt->geometry().volume();
+// totalVolume += volume;
+//
+// meanPermeability += volume / perm;
+// }
+//
+// for (ElementIterator eIt = gridView.template begin<0> (); eIt
+// != eItEnd; ++eIt)
+// {
+// Scalar perm = randomPermeability.K(*eIt)[0][0];
+// permeability_[indexSet_.index(*(eIt->template subEntity<dim> (0)))]
+// = perm;
+//
+// Scalar volume = eIt->geometry().volume();
+// totalVolume += volume;
+//
+// meanPermeability += volume / perm;
+// }
+//
+// meanPermeability /= totalVolume;
+// meanPermeability = 1.0 / meanPermeability;
+//
+// //Iterate over elements
+// VertexIterator vItEnd = gridView.template end<dim> ();
+// for (VertexIterator vIt = gridView.template begin<dim> (); vIt
+// != vItEnd; ++vIt)
+// {
+// int globalIdx = indexSet_.index(*vIt);
+//
+// Scalar elementEntryPressure = 1./GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, VgAlpha2);
+// elementEntryPressure *= sqrt(meanPermeability / permeability_[globalIdx]);
+//
+// materialLawParams_[globalIdx].setVgAlpha(1./elementEntryPressure);
+// materialLawParams_[globalIdx].setVgn(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, VgN2));
+// materialLawParams_[globalIdx].setSwr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Swr2));
+// materialLawParams_[globalIdx].setSnr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Snr2));
+// }
+//
+// // if (create)
+// // {
+// Dune::VTKWriter<GridView> vtkwriter(gridView);
+// vtkwriter.addVertexData(permeability_, "absolute permeability");
+// PermeabilityType logPerm(size);
+// for (unsigned i = 0; i < size; i++)
+// logPerm[i] = log10(permeability_[i]);
+// vtkwriter.addVertexData(logPerm, "logarithm of permeability");
+// vtkwriter.write("permeability", Dune::VTK::OutputType::ascii);
+ }
+
+ //
+ /*!
+ * \brief This is called from the coupled problem
+ * and creates a gnuplot output of the Pc-Sw curve -- currently commented
+ */
+ 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));
+// }
+// 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));
+// }
+// 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));
+// }
+ }
+
+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_;
+ const IndexSet& indexSet_;
+
+ MaterialLawParams coarseParams_;
+ MaterialLawParams mediumParams_;
+ MaterialLawParams fineParams_;
+};
+
+} // end namespace
+
+#endif // DUMUX_TWOCSTOKES_2P2C_SPATIALPARAMS_HH
diff --git a/test/multidomain/2cstokes2p2c/2p2csubproblem.hh b/test/multidomain/2cstokes2p2c/2p2csubproblem.hh
new file mode 100644
index 0000000000..de6edfcaff
--- /dev/null
+++ b/test/multidomain/2cstokes2p2c/2p2csubproblem.hh
@@ -0,0 +1,459 @@
+// -*- 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_2P2C_SUBPROBLEM_HH
+#define DUMUX_2P2C_SUBPROBLEM_HH
+
+#include <dumux/implicit/common/implicitporousmediaproblem.hh>
+#include <dumux/multidomain/couplinglocalresiduals/2p2ccouplinglocalresidual.hh>
+#include <dumux/multidomain/common/subdomainpropertydefaults.hh>
+#include <dumux/multidomain/common/multidomainlocaloperator.hh>
+
+#include "2cstokes2p2cproblem.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class TwoPTwoCSubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(TwoPTwoCSubProblem,
+ INHERITS_FROM(BoxTwoPTwoC, SubDomain, TwoCStokesTwoPTwoCSpatialParams));
+
+// Set the problem property
+SET_TYPE_PROP(TwoPTwoCSubProblem, Problem, TwoPTwoCSubProblem<TTAG(TwoPTwoCSubProblem)>);
+
+// Set the local residual extended for the coupling
+SET_TYPE_PROP(TwoPTwoCSubProblem, LocalResidual, TwoPTwoCCouplingLocalResidual<TypeTag>);
+
+// choose pn and Sw as primary variables
+SET_INT_PROP(TwoPTwoCSubProblem, Formulation, TwoPTwoCFormulation::pnsw);
+
+
+// the gas component balance (air) is replaced by the total mass balance
+SET_PROP(TwoPTwoCSubProblem, ReplaceCompEqIdx)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ static const int value = Indices::contiNEqIdx;
+};
+
+// set the constraints for multidomain / pdelab // not required here
+SET_TYPE_PROP(TwoPTwoCSubProblem, Constraints,
+ Dune::PDELab::NonoverlappingConformingDirichletConstraints);
+
+// set the grid operator
+SET_PROP(TwoPTwoCSubProblem, 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(TwoPTwoCSubProblem, 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(TwoPTwoCSubProblem, UseMoles, false);
+
+// enable/disable velocity output
+SET_BOOL_PROP(TwoPTwoCSubProblem, VtkAddVelocity, true);
+
+// Enable gravity
+SET_BOOL_PROP(TwoPTwoCSubProblem, ProblemEnableGravity, true);
+}
+
+template <class TypeTag = TTAG(TwoPTwoCSubProblem) >
+class TwoPTwoCSubProblem : 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 TwoPTwoCSubProblem<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 { // the equation indices
+ contiTotalMassIdx = Indices::contiNEqIdx,
+ contiWEqIdx = Indices::contiWEqIdx,
+ };
+ enum { // the indices of the primary variables
+ pressureIdx = Indices::pressureIdx,
+ switchIdx = Indices::switchIdx,
+ };
+ enum { // the indices for the phase presence
+ wPhaseOnly = Indices::wPhaseOnly,
+ nPhaseOnly = Indices::nPhaseOnly,
+ bothPhases = Indices::bothPhases
+ };
+ 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 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:
+ TwoPTwoCSubProblem(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);
+ refTemperatureFF_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefTemperature);
+ 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);
+
+ freqMassOutput_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqMassOutput);
+
+ storageLastTimestep_ = Scalar(0);
+ lastMassOutputTime_ = Scalar(0);
+
+ outfile.open("evaporation.out");
+ outfile << "time; evaporationRate; massH2O" << std::endl;
+ }
+ catch (Dumux::ParameterException &e) {
+ std::cerr << e << ". Abort!\n";
+ exit(1) ;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ exit(1);
+ }
+ }
+
+ ~TwoPTwoCSubProblem()
+ {
+ 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;
+ }
+
+ //! \copydoc Dumux::ImplicitPorousMediaProblem::temperatureAtPos()
+ Scalar temperatureAtPos(const GlobalPosition &globalPos) const
+ {
+ return refTemperature_;
+ };
+
+ // \}
+
+ //! \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 (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 = Scalar(0);
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ //! \copydoc Dumux::ImplicitProblem::source()
+ void source(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ 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);
+
+ values = Scalar(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;
+ }
+
+ //! \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]
+ << " EvaporationRate: " << evaporationRate[contiWEqIdx]
+ << " Time: " << time <<std::endl;
+ if (this->timeManager().time() != 0.)
+ outfile << time <<
+ "; " << evaporationRate[contiTotalMassIdx] <<
+ "; " << evaporationRate[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 globalPos 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
+ {
+ values[pressureIdx] = refPressure_
+ + 1000.*this->gravity()[1]*(globalPos[1]-bboxMax_[1]);
+ if (globalPos[0] < xMaterialInterface_)
+ values[switchIdx] = initialSw1_;
+ else
+ values[switchIdx] = initialSw2_;
+
+// if (onUpperBoundary_(globalPos) && onLeftBoundary_(globalPos))
+// {
+ // values[switchIdx] = 0.99e-4;
+// }
+
+ // values[pressureIdx] = 1e5*(2.0 - globalPos[0]);
+ // if (onLeftBoundary_(globalPos))
+ // values[switchIdx] = 0.9e-4;
+ // else
+ // values[switchIdx] = 1e-4;
+ }
+
+ 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 refTemperatureFF_;
+ Scalar refPressure_;
+ Scalar initialSw1_;
+ Scalar initialSw2_;
+
+ Scalar runUpDistanceX_;
+ Scalar initializationTime_;
+ std::ofstream outfile;
+};
+} //end namespace
+
+#endif // DUMUX_2P2C_SUBPROBLEM_HH
diff --git a/test/multidomain/2cstokes2p2c/Makefile.am b/test/multidomain/2cstokes2p2c/Makefile.am
new file mode 100644
index 0000000000..ac442bfc93
--- /dev/null
+++ b/test/multidomain/2cstokes2p2c/Makefile.am
@@ -0,0 +1,5 @@
+# tests where program to build and program to run are equal
+check_PROGRAMS = test_2cstokes2p2c
+test_2cstokes2p2c_SOURCES = test_2cstokes2p2c.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/multidomain/2cstokes2p2c/grids/interfacedomain.dgf b/test/multidomain/2cstokes2p2c/grids/interfacedomain.dgf
new file mode 100644
index 0000000000..b0290172cb
--- /dev/null
+++ b/test/multidomain/2cstokes2p2c/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/2cstokes2p2c/grids/test_2cstokes2p2c.dgf b/test/multidomain/2cstokes2p2c/grids/test_2cstokes2p2c.dgf
new file mode 100644
index 0000000000..3085c34981
--- /dev/null
+++ b/test/multidomain/2cstokes2p2c/grids/test_2cstokes2p2c.dgf
@@ -0,0 +1,17 @@
+DGF
+Interval
+0 0 % first corner
+0.25 0.5 % second corner
+32 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/2cstokes2p2c/stokes2csubproblem.hh b/test/multidomain/2cstokes2p2c/stokes2csubproblem.hh
new file mode 100644
index 0000000000..f523d7537f
--- /dev/null
+++ b/test/multidomain/2cstokes2p2c/stokes2csubproblem.hh
@@ -0,0 +1,598 @@
+// -*- 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 Stokes2cProblems
+ * \brief Definition of an isothermal compositional Stokes problem
+ */
+#ifndef DUMUX_STOKES2C_SUBPROBLEM_HH
+#define DUMUX_STOKES2C_SUBPROBLEM_HH
+
+#include <dumux/freeflow/stokesnc/stokesncmodel.hh>
+#include <dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh>
+#include <dumux/multidomain/common/subdomainpropertydefaults.hh>
+
+namespace Dumux
+{
+
+template <class TypeTag>
+class Stokes2cSubProblem;
+
+//////////
+// Specify the properties for the Stokes problem
+//////////
+namespace Properties
+{
+NEW_TYPE_TAG(Stokes2cSubProblem,
+ INHERITS_FROM(BoxStokesnc, SubDomain, TwoCStokesTwoPTwoCSpatialParams));
+
+// Set the problem property
+SET_TYPE_PROP(Stokes2cSubProblem, Problem, Dumux::Stokes2cSubProblem<TypeTag>);
+
+/*!
+ * \brief Set the property for the material parameters by extracting
+ * it from the material law.
+ */
+SET_PROP(Stokes2cSubProblem, MaterialLawParams)
+{
+ private:
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+ public:
+ typedef typename MaterialLaw::Params type;
+};
+
+//! Use the Stokes2cCouplingLocalResidual for the computation of the local residual in the Stokes domain
+SET_TYPE_PROP(Stokes2cSubProblem,
+ LocalResidual,
+ StokesncCouplingLocalResidual<TypeTag>);
+
+SET_TYPE_PROP(Stokes2cSubProblem, Constraints,
+ Dune::PDELab::NonoverlappingConformingDirichletConstraints);
+
+SET_PROP(Stokes2cSubProblem, 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(Stokes2cSubProblem, FluidSystem)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, MultiDomainTypeTag) CoupledTypeTag;
+ typedef typename GET_PROP_TYPE(CoupledTypeTag, FluidSystem) FluidSystem;
+public:
+ typedef FluidSystem type;
+};
+
+//! Set Scalar to type long double for higher accuracy
+SET_TYPE_PROP(BoxStokes, Scalar, double);
+//SET_TYPE_PROP(BoxStokes, Scalar, long double);
+
+// use formulation based on mass fractions
+SET_BOOL_PROP(Stokes2cSubProblem, UseMoles, false);
+
+// Disable gravity
+SET_BOOL_PROP(Stokes2cSubProblem, ProblemEnableGravity, false);
+
+// switch inertia term on or off
+SET_BOOL_PROP(Stokes2cSubProblem, EnableNavierStokes, false);
+}
+
+/*!
+ * \ingroup Stokes2cBoxProblems
+ * \brief Stokes2c 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 Stokes2cModel. It is part of the 2cstokes2p2c model and
+ * is combined with the 2p2csubproblem for the Darcy domain.
+ *
+ */
+template <class TypeTag>
+class Stokes2cSubProblem : public StokesProblem<TypeTag>
+{
+ typedef Stokes2cSubProblem<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;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ enum {
+
+ // Number of equations and grid dimension
+ numEq = GET_PROP_VALUE(TypeTag, NumEq),
+ dim = GridView::dimension,
+
+ // copy some indices for convenience
+ 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)
+ };
+ enum { // primary variable indices
+ pressureIdx = Indices::pressureIdx,
+ velocityXIdx = Indices::velocityXIdx,
+ velocityYIdx = Indices::velocityYIdx,
+ velocityZIdx = Indices::velocityZIdx,
+ massOrMoleFracIdx = Indices::massOrMoleFracIdx
+ };
+ enum { phaseIdx = Indices::phaseIdx };
+ enum { numComponents = Indices::numComponents };
+// enum { numComponents = GET_PROP_VALUE(TypeTag, 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;
+
+
+public:
+ /*!
+ * \brief docme
+ *
+ * \param timeManager The time manager
+ * \param gridView The grid view
+ *
+ */
+ Stokes2cSubProblem(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);
+
+ 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); }
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature() const
+ {
+ return refTemperature_;
+ };
+
+ // \}
+
+ /*!
+ * \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);
+
+ // Left inflow boundaries should be Neumann, otherwise the
+ // evaporative fluxes are much more grid dependent
+ if (onLeftBoundary_(globalPos))
+ {
+ values.setNeumann(transportEqIdx);
+
+ if (onUpperBoundary_(globalPos)) // corner point
+ values.setAllDirichlet();
+ }
+
+
+ if (onRightBoundary_(globalPos))
+ {
+ values.setAllOutflow();
+
+ if (onUpperBoundary_(globalPos)) // corner point
+ values.setAllDirichlet();
+ }
+
+ if (onLowerBoundary_(globalPos))
+ {
+ values.setAllDirichlet();
+ values.setNeumann(transportEqIdx);
+
+ 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 all corner points to dirichlet
+ // if ((onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)) &&
+ // (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos)))
+ // {
+ // values.setAllDirichlet();
+ //// values.setCouplingOutflow(momentumXIdx);
+ // }
+
+ // 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);
+
+ 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_;
+ }
+
+ //! \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 xVelocity = xVelocity_(globalPos);
+
+ if (onLeftBoundary_(globalPos)
+ && globalPos[1] > bboxMin_[1] && globalPos[1] < bboxMax_[1])
+ {
+ values[transportEqIdx] = -xVelocity*density*refMassfrac_;
+
+// (globalPos[1] - bboxMin_[1])*(bboxMax_[1] - globalPos[1])
+// / (0.25*height_()*height_()));
+ }
+ }
+
+ /*!
+ * \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);
+ }
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ //! \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 globalPos 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_; }
+
+ /*!
+ * \brief Returns the reference pressure.
+ */
+ const Scalar refPressure() const
+ { return refPressure_; }
+
+ const Scalar refTemperature() const
+ { return refTemperature_; }
+
+ const Scalar refMassfrac() const
+ { return refMassfrac_; }
+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_;
+ }
+
+ 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_;
+ // 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();
+ 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 sinusVelVar_;
+ Scalar alphaBJ_;
+
+ Scalar xMaterialInterface_;
+ Scalar runUpDistanceX_;
+ Scalar initializationTime_;
+};
+} //end namespace
+
+#endif // DUMUX_STOKES2C_SUBPROBLEM_HH
diff --git a/test/multidomain/2cstokes2p2c/test_2cstokes2p2c.cc b/test/multidomain/2cstokes2p2c/test_2cstokes2p2c.cc
new file mode 100644
index 0000000000..03bb706839
--- /dev/null
+++ b/test/multidomain/2cstokes2p2c/test_2cstokes2p2c.cc
@@ -0,0 +1,252 @@
+// -*- 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 2c-Stokes and 2p2c-Darcy model
+ */
+
+#include "config.h"
+#include <iostream>
+
+#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 "2cstokes2p2cproblem.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_2cstokes2p2c.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") || 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());
+
+ // instantiate coupled problem
+ 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;
+ }
+}
+
+int main(int argc, char** argv)
+{
+ typedef TTAG(TwoCStokesTwoPTwoCProblem) ProblemTypeTag;
+ return start<ProblemTypeTag>(argc, argv);
+}
diff --git a/test/multidomain/2cstokes2p2c/test_2cstokes2p2c.input b/test/multidomain/2cstokes2p2c/test_2cstokes2p2c.input
new file mode 100644
index 0000000000..819ae289ff
--- /dev/null
+++ b/test/multidomain/2cstokes2p2c/test_2cstokes2p2c.input
@@ -0,0 +1,170 @@
+#############################################################
+#Configuration file for test_2cstokes2p2c
+#############################################################
+
+#############################################################
+#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 = 180
+
+#Initialization time without coupling
+#set to zero or negative if not desired
+InitTime = 864
+
+#Simulation end
+TEnd= 1.3e6
+
+#Define the length of an episode (for output)
+EpisodeLength = 43200
+
+#############################################################
+[Vtk]
+#############################################################
+#Names for VTK output
+NameFF = stokes2c
+NamePM = 2p2c
+AddVelocity = 1
+
+#############################################################
+[Grid]
+#############################################################
+UseInterfaceMeshCreator = true
+
+#Name of the dgf file (grid)
+File = grids/interfacedomain.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
+SinusVelVar = 0.0
+ExponentMTC = 0.0 # 1./6., Mass transfer coefficient for S^MTC
+UseBoundaryLayerModel = 9 # integer, 0 for no boundary layer model, 1 for Blasius, 2 for turbulent BL, 3 for constant thickness
+BoundaryLayerOffset = 0.0 # For BL model like Blasius, determines a virtual run-up distance for the flow
+ConstThickness = 0.0018 # for a constant BL thickness, use BL model 3
+MassTransferModel = 2 # 0 for none, 1 for power law, 2 for Schluender model
+
+#############################################################
+[PorousMedium]
+#############################################################
+RefPressurePM = 1e5
+RefTemperaturePM = 298.15
+InitialSw1 = 0.98
+InitialSw2 = 0.98
+CharPoreDiameter = 1e-4
+
+#############################################################
+[SpatialParams]
+#############################################################
+MaterialInterfaceX = 100.0
+AlphaBJ = 1.0
+
+# for homogeneous setups (0 for heterogeneous):
+SoilType = 2
+RegularizationThreshold = 1e-2
+
+# GStat stuff, only required for SoilType=0
+GenerateNewPermeability = true
+GStatControlFileName = gstatControl_2D.txt
+GStatInputFileName = gstatInput.txt
+PermeabilityInputFileName = permeab.dat
+
+### SoilType = 1 ###
+[SpatialParams.Coarse]
+Permeability1 = 7.0e-10
+Porosity1 = 0.44
+Swr1 = 0.005
+Snr1 = 0.01
+VgAlpha1 = 8.74e-4
+VgN1 = 3.35
+PlotMaterialLaw1 = false
+LambdaSolid1 = 5.3
+
+### SoilType = 2 ###
+[SpatialParams.Medium]
+Permeability2 = 2.65e-10
+Porosity2 = 0.41
+Swr2 = 0.005
+Snr2 = 0.01
+VgAlpha2 = 6.371e-4
+VgN2 = 6.9
+PlotMaterialLaw2 = false
+LambdaSolid2 = 5.3
+
+### SoilType = 3 ###
+[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.3
+
+#############################################################
+[Newton]
+#############################################################
+RelTolerance = 1e-5
+TargetSteps = 8
+MaxSteps = 12
+WriteConvergence = false
+MaxTimeStepDivisions = 20
+
+#############################################################
+[LinearSolver]
+#############################################################
+ResidualReduction = 1e-10
+Verbosity = 0
+MaxIterations = 200
+
diff --git a/test/multidomain/Makefile.am b/test/multidomain/Makefile.am
new file mode 100644
index 0000000000..08f06370a2
--- /dev/null
+++ b/test/multidomain/Makefile.am
@@ -0,0 +1,7 @@
+SUBDIRS = \
+ 2cstokes2p2c
+
+EXTRA_DIST = CMakeLists.txt
+
+include $(top_srcdir)/am/global-rules
+
--
GitLab