"README.md" did not exist on "d57c55199b7be6fffd04c202f4c6f3ac07a6f520"
Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
---
title: The DuMux Multidomain Framework
author: The DuMux Development Team, IWS-LH2, University of Stuttgart
date: April 03, 2023
---
# Coupling Multiple Domains
## Coupling Multiple Domains
{width=900px}
## Design Goals
* Reuse existing DuMux models in the coupled setting
* Extensible to more than two subdomains
* Common structure for different coupling modes
# General Structure
## General Structure
## The Sub Problems
* any DuMux problem with
- initial conditions
- boundary conditions
- associated DuMux model
* get a pointer to the coupling manager
## The Coupling Manager
* Transfer data from one subproblem to another
- e.g. give me the soil pressure (from the well domain)
- e.g. give me the radius of the embedded well (from the soil domain)
* Compute the coupling stencil
* Compute the coupling residual (numerical differentiation)
# Example: Embedded Fracture Model
## Coupling Residual
$$
\begin{aligned}
\frac{\partial \rho \phi}{\partial t} - \nabla\cdot \left[ \frac{\rho}{\mu}\mathbf{K_m} (\nabla p_m - \rho \mathbf{g}) \right] - q \delta_\Gamma = 0 \quad \text{in} \, \Omega \\
\frac{\partial \rho \phi}{\partial t} - \nabla_\tau\cdot \left[ \frac{\rho}{\mu}\mathbf{K_f} (\nabla p_f - \rho \mathbf{g}) \right] + q = 0 \quad \text{in} \, \Lambda \\
q = \zeta (p_f - p_m) \quad \text{and} \quad \Gamma = \Omega \cap \Lambda
\end{aligned}
$$
## Coupling Residual / Data Transfer (1/3)
```c++
/*!*
Evaluates the point sources (added by addPointSources)
* for all phases within a given sub-control volume.!*/
template<class ElementVolumeVariables>
void pointSource(PointSource& source...) const
{
// compute source at every integration point
const Scalar pressure3D = this->couplingManager()
.bulkPriVars(source.id())[Indices::pressureIdx];
const Scalar pressure1D = this->couplingManager()
.lowDimPriVars(source.id())[Indices::pressureIdx];
...
```
## Coupling Residual / Data Transfer (2/3)
```c++
template<class ElementVolumeVariables>
void pointSource(PointSource& source...) const
{
...
// calculate the source
const Scalar meanDistance = this->couplingManager()
.averageDistance(source.id());
static const Scalar matrixPerm = getParamFromGroup<Scalar>("...");
static const Scalar rho = getParam<Scalar>("...");
static const Scalar mu = getParam<Scalar>("...")*rho;
...
}
```
## Coupling Residual / Data Transfer (3/3)
```c++
template<class ElementVolumeVariables>
void pointSource(PointSource& source...) const
{
...
const Scalar sourceValue = rho*(pressure3D - pressure1D)
/meanDistance*matrixPerm/mu;
source = sourceValue*source.quadratureWeight()
*source.integrationElement();
}
```
## Coupling Residual (1/2)
```c++
/*! evaluates the element residual of a coupled element of
domain i which depends on the variables at the degree of
freedom with index dofIdxGlobalJ of domain j */
template<std::size_t i, std::size_t j, class LocalAssemblerI>
decltype(auto) evalCouplingResidual( domainI ,localAssemblerI
,domainJ ,dofIdxGlobalJ)
{
static_assert(i != j, "A domain cannot be coupled to itself!");
typename LocalAssemblerI::LocalResidual::ElementResidualVector residual;
const auto& element = localAssemblerI.element();
const auto& fvGeometry = localAssemblerI.fvGeometry();
const auto& curElemVolVars = localAssemblerI.curElemVolVars();
residual.resize(fvGeometry.numScv()); ...}
```
## Coupling Residual (2/2)
```c++
decltype(auto) evalCouplingResidual( domainI ,localAssemblerI
,domainJ ,dofIdxGlobalJ)
{...
for (const auto& scv : scvs(fvGeometry))
{
auto couplingSource = this->problem(domainI)
.scvPointSources(...);
couplingSource += this->problem(domainI)
.source(...);
couplingSource *= -GridGeometry<i>::Extrusion::volume(...)
*curElemVolVars[scv].extrusionFactor();
residual[scv.indexInElement()] = couplingSource;
}
return residual;}
```
## Coupling Stencil
## Coupling Stencil (1/2)
```c++
/*!
* returns an iterable container of all indices of degrees of
*freedom of domain j that couple with / influence the element
*residual of the given element of domain i */
template<std::size_t i, std::size_t j>
const CouplingStencil<j>& couplingStencil(domainI,
const Element<i>& element,
domainJ) const
{
static_assert(i != j, "A domain cannot be coupled to itself!");
const auto eIdx = this->problem(domainI).gridGeometry()
.elementMapper().index(element);
...}
```
## Coupling Stencil (2/2)
```c++
template<std::size_t i, std::size_t j>
const CouplingStencil<j>& couplingStencil(domainI,
const Element<i>& element,
domainJ) const
{...
if (couplingStencils(domainI).count(eIdx))
return couplingStencils(domainI).at(eIdx);
else
return emptyStencil(domainI);
}
```
# Examples
## Embedded Fracture Model
## Root Water Uptake
