3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
multidomain/boundary/darcydarcy/couplingmanager.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
25#ifndef DUMUX_MULTIDOMAIN_DARCYDARCY_BOUNDARY_COUPLINGMANAGER_HH
26#define DUMUX_MULTIDOMAIN_DARCYDARCY_BOUNDARY_COUPLINGMANAGER_HH
27
28#include <iostream>
29#include <vector>
30
31#include <dune/common/indices.hh>
32#include <dune/common/exceptions.hh>
33
35#include <dumux/common/math.hh>
41
42namespace Dumux {
43
48template<class MDTraits>
50: public CouplingManager<MDTraits>
51{
53
54 using Scalar = typename MDTraits::Scalar;
55 using SolutionVector = typename MDTraits::SolutionVector;
56
57 template<std::size_t i> using SubDomainTypeTag = typename MDTraits::template SubDomain<i>::TypeTag;
58 template<std::size_t i> using Problem = GetPropType<SubDomainTypeTag<i>, Properties::Problem>;
59 template<std::size_t i> using PrimaryVariables = GetPropType<SubDomainTypeTag<i>, Properties::PrimaryVariables>;
60 template<std::size_t i> using NumEqVector = GetPropType<SubDomainTypeTag<i>, Properties::NumEqVector>;
61 template<std::size_t i> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<i>, Properties::GridVolumeVariables>::LocalView;
62 template<std::size_t i> using VolumeVariables = typename GetPropType<SubDomainTypeTag<i>, Properties::GridVolumeVariables>::VolumeVariables;
63 template<std::size_t i> using GridGeometry = typename MDTraits::template SubDomain<i>::GridGeometry;
64 template<std::size_t i> using FVElementGeometry = typename GridGeometry<i>::LocalView;
65 template<std::size_t i> using SubControlVolumeFace = typename GridGeometry<i>::SubControlVolumeFace;
66 template<std::size_t i> using SubControlVolume = typename GridGeometry<i>::SubControlVolume;
67 template<std::size_t i> using GridView = typename GridGeometry<i>::GridView;
68 template<std::size_t i> using Element = typename GridView<i>::template Codim<0>::Entity;
69
70 template<std::size_t i>
71 static constexpr auto domainIdx()
72 { return typename MDTraits::template SubDomain<i>::Index{}; }
73
74 template<std::size_t i>
75 static constexpr bool isCCTpfa()
76 { return GridGeometry<i>::discMethod == DiscretizationMethod::cctpfa; }
77
78 using CouplingStencil = std::vector<std::size_t>;
79public:
81 using MultiDomainTraits = MDTraits;
83 using CouplingStencils = std::unordered_map<std::size_t, CouplingStencil>;
84
88 // \{
89
90 void init(std::shared_ptr<Problem<domainIdx<0>()>> problem0,
91 std::shared_ptr<Problem<domainIdx<1>()>> problem1,
92 const SolutionVector& curSol)
93 {
94 static_assert(isCCTpfa<0>() && isCCTpfa<1>(), "Only cctpfa implemented!");
95
96 this->setSubProblems(std::make_tuple(problem0, problem1));
97 this->updateSolution(curSol);
98 couplingMapper_.update(*this);
99 }
100
101 // \}
102
106 // \{
107
116 template<std::size_t i, std::size_t j>
117 const CouplingStencil& couplingStencil(Dune::index_constant<i> domainI,
118 const Element<i>& element,
119 Dune::index_constant<j> domainJ) const
120 {
121 static_assert(i != j, "A domain cannot be coupled to itself!");
122 const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element);
123 return couplingMapper_.couplingStencil(domainI, eIdx, domainJ);
124 }
125
126 // \}
127
133 template<std::size_t i>
134 bool isCoupled(Dune::index_constant<i> domainI,
135 const SubControlVolumeFace<i>& scvf) const
136 {
137 return couplingMapper_.isCoupled(domainI, scvf);
138 }
139
156 template<std::size_t i>
157 Scalar advectiveFluxCoupling(Dune::index_constant<i> domainI,
158 const Element<i>& element,
159 const FVElementGeometry<i>& fvGeometry,
160 const ElementVolumeVariables<i>& elemVolVars,
161 const SubControlVolumeFace<i>& scvf,
162 int phaseIdx = 0) const
163 {
164 Scalar flux = 0.0;
165 static const bool enableGravity = getParamFromGroup<bool>(this->problem(domainI).paramGroup(), "Problem.EnableGravity");
166 constexpr auto otherDomainIdx = domainIdx<1-i>();
167
168 const auto& outsideElement = this->problem(otherDomainIdx).gridGeometry().element(couplingMapper_.outsideElementIndex(domainI, scvf));
169 auto fvGeometryOutside = localView(this->problem(otherDomainIdx).gridGeometry());
170 fvGeometryOutside.bindElement(outsideElement);
171
172 const auto& flipScvf = fvGeometryOutside.scvf(couplingMapper_.flipScvfIndex(domainI, scvf));
173 const auto& outsideScv = fvGeometryOutside.scv(flipScvf.insideScvIdx());
174 const auto outsideVolVars = volVars(otherDomainIdx, outsideElement, outsideScv);
175
176 const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
177 const auto& insideVolVars = elemVolVars[insideScv];
178
179 const auto ti = computeTpfaTransmissibility(scvf, insideScv, insideVolVars.permeability(), insideVolVars.extrusionFactor());
180 const auto tj = computeTpfaTransmissibility(flipScvf, outsideScv, outsideVolVars.permeability(), outsideVolVars.extrusionFactor());
181 Scalar tij = 0.0;
182 if (ti*tj > 0.0)
183 tij = scvf.area()*(ti * tj)/(ti + tj);
184
185 if (enableGravity)
186 {
187 // do averaging for the density over all neighboring elements
188 const auto rho = (insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx))*0.5;
189 const auto& g = this->problem(domainI).spatialParams().gravity(scvf.ipGlobal());
190
192 const auto alpha_inside = vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g)*insideVolVars.extrusionFactor();
193 const auto alpha_outside = vtmv(scvf.unitOuterNormal(), outsideVolVars.permeability(), g)*outsideVolVars.extrusionFactor();
194
195 // Obtain inside and outside pressures
196 const auto pInside = insideVolVars.pressure(0);
197 const auto pOutside = outsideVolVars.pressure(0);
198
199 flux = tij*(pInside - pOutside) + rho*scvf.area()*alpha_inside - rho*tij/tj*(alpha_inside - alpha_outside);
200
201 }
202 else
203 {
204 // Obtain inside and outside pressures
205 const auto pInside = insideVolVars.pressure(phaseIdx);
206 const auto pOutside = outsideVolVars.pressure(phaseIdx);
207
208 // return flux
209 flux = tij*(pInside - pOutside);
210 }
211
212 // upwind scheme
213 static const Scalar upwindWeight = getParam<Scalar>("Flux.UpwindWeight");
214 auto upwindTerm = [phaseIdx](const auto& volVars){ return volVars.density(phaseIdx)*volVars.mobility(phaseIdx); };
215 if (std::signbit(flux)) // if sign of flux is negative
216 flux *= (upwindWeight*upwindTerm(outsideVolVars)
217 + (1.0 - upwindWeight)*upwindTerm(insideVolVars));
218 else
219 flux *= (upwindWeight*upwindTerm(insideVolVars)
220 + (1.0 - upwindWeight)*upwindTerm(outsideVolVars));
221
222 // make this a area-specific flux
223 flux /= scvf.area()*insideVolVars.extrusionFactor();
224 return flux;
225 }
226
228 template<std::size_t i>
229 VolumeVariables<i> volVars(Dune::index_constant<i> domainI,
230 const Element<i>& element,
231 const SubControlVolume<i>& scv) const
232 {
233 VolumeVariables<i> volVars;
234 const auto elemSol = elementSolution(element, this->curSol()[domainI], this->problem(domainI).gridGeometry());
235 volVars.update(elemSol, this->problem(domainI), element, scv);
236 return volVars;
237 }
238
239private:
241};
242
243} // end namespace Dumux
244
245#endif
Define some often used mathematical functions.
Element solution classes and factory functions.
The available discretization methods in Dumux.
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:38
auto elementSolution(const Element &element, const SolutionVector &sol, const GridGeometry &gg) -> std::enable_if_t< GridGeometry::discMethod==DiscretizationMethod::box, BoxElementSolution< typename GridGeometry::LocalView, std::decay_t< decltype(std::declval< SolutionVector >()[0])> > >
Make an element solution for box schemes.
Definition: box/elementsolution.hh:115
Tensor::field_type computeTpfaTransmissibility(const SubControlVolumeFace &scvf, const SubControlVolume &scv, const Tensor &T, typename SubControlVolume::Traits::Scalar extrusionFactor)
Free function to evaluate the Tpfa transmissibility associated with the flux (in the form of flux = T...
Definition: tpfa/computetransmissibility.hh:47
Dune::DenseMatrix< MAT >::value_type vtmv(const Dune::DenseVector< V1 > &v1, const Dune::DenseMatrix< MAT > &M, const Dune::DenseVector< V2 > &v2)
Evaluates the scalar product of a vector v2, projected by a matrix M, with a vector v1.
Definition: math.hh:840
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
A vector of primary variables.
Definition: common/properties.hh:59
A vector of size number equations that can be used for Neumann fluxes, sources, residuals,...
Definition: common/properties.hh:61
Property to specify the type of a problem which has to be solved.
Definition: common/properties.hh:69
The type for a global container for the volume variables.
Definition: common/properties.hh:176
Coupling manager for equal-dimension boundary coupling of darcy models.
Definition: multidomain/boundary/darcydarcy/couplingmanager.hh:51
std::unordered_map< std::size_t, CouplingStencil > CouplingStencils
export stencil types
Definition: multidomain/boundary/darcydarcy/couplingmanager.hh:83
Scalar advectiveFluxCoupling(Dune::index_constant< i > domainI, const Element< i > &element, const FVElementGeometry< i > &fvGeometry, const ElementVolumeVariables< i > &elemVolVars, const SubControlVolumeFace< i > &scvf, int phaseIdx=0) const
Evaluate the boundary conditions for a coupled Neumann boundary segment.
Definition: multidomain/boundary/darcydarcy/couplingmanager.hh:157
const CouplingStencil & couplingStencil(Dune::index_constant< i > domainI, const Element< i > &element, Dune::index_constant< j > domainJ) const
Methods to be accessed by the assembly.
Definition: multidomain/boundary/darcydarcy/couplingmanager.hh:117
VolumeVariables< i > volVars(Dune::index_constant< i > domainI, const Element< i > &element, const SubControlVolume< i > &scv) const
Return the volume variables of domain i for a given element and scv.
Definition: multidomain/boundary/darcydarcy/couplingmanager.hh:229
void init(std::shared_ptr< Problem< domainIdx< 0 >()> > problem0, std::shared_ptr< Problem< domainIdx< 1 >()> > problem1, const SolutionVector &curSol)
Methods to be accessed by main.
Definition: multidomain/boundary/darcydarcy/couplingmanager.hh:90
bool isCoupled(Dune::index_constant< i > domainI, const SubControlVolumeFace< i > &scvf) const
If the boundary entity is on a coupling boundary.
Definition: multidomain/boundary/darcydarcy/couplingmanager.hh:134
MDTraits MultiDomainTraits
export traits
Definition: multidomain/boundary/darcydarcy/couplingmanager.hh:81
the default mapper for conforming equal dimension boundary coupling between two domains (box or cc)
Definition: boundary/darcydarcy/couplingmapper.hh:49
Definition: multidomain/couplingmanager.hh:46
void setSubProblems(const std::tuple< std::shared_ptr< SubProblems >... > &problems)
set the pointers to the sub problems
Definition: multidomain/couplingmanager.hh:247
const Problem< i > & problem(Dune::index_constant< i > domainIdx) const
Return a reference to the sub problem.
Definition: multidomain/couplingmanager.hh:264
SolutionVector & curSol()
the solution vector of the coupled problem
Definition: multidomain/couplingmanager.hh:278
void updateSolution(const SolutionVector &curSol)
Updates the entire solution vector, e.g. before assembly or after grid adaption.
Definition: multidomain/couplingmanager.hh:186
Declares all properties used in Dumux.
Free functions to evaluate the transmissibilities associated with flux evaluations across sub-control...
The interface of the coupling manager for multi domain problems.