GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/analyticalsolutions.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 175 177 98.9%
Functions: 18 20 90.0%
Branches: 0 0 -%
Line Branch Exec Source
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 // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup BoundaryTests
10 * \brief Analytical solutions and rhs for the different test cases
11 */
12
13 #ifndef DUMUX_CONVERGENCE_TEST_ANALYTICAL_SOLUTIONS_HH
14 #define DUMUX_CONVERGENCE_TEST_ANALYTICAL_SOLUTIONS_HH
15
16 #include <cmath>
17 #include <dune/common/fvector.hh>
18 #include <dune/common/fmatrix.hh>
19
20 namespace Dumux::Solution::DarcyStokes {
21
22 /////////////////////////////////////////////////////////////////////////////////////////////////
23 // see Rybak et al., 2015: "Multirate time integration for coupled
24 // saturated/unsaturated porous medium and free flow systems"
25 /////////////////////////////////////////////////////////////////////////////////////////////////
26 namespace Rybak {
27
28 // 0: velocity-x | 1: velocity-y | 2: pressure
29 72288 Dune::FieldVector<double, 3> darcy(const Dune::FieldVector<double, 2>& globalPos)
30 {
31 72288 Dune::FieldVector<double, 3> sol(0.0);
32 144576 const double x = globalPos[0];
33 144576 const double y = globalPos[1];
34
35 72288 using std::exp; using std::sin; using std::cos;
36 144576 sol[0] = -0.5*M_PI*y*y*cos(M_PI*x);
37 144576 sol[1] = -1.0*y*sin(M_PI*x);
38 144576 sol[2] = 0.5*y*y*sin(M_PI*x);
39 72288 return sol;
40 }
41
42 // 0: mass
43 136640 Dune::FieldVector<double, 1> darcyRHS(const Dune::FieldVector<double, 2>& globalPos)
44 {
45 273280 const double x = globalPos[0];
46 273280 const double y = globalPos[1];
47 136640 using std::sin;
48 273280 return { (0.5*M_PI*y*M_PI*y - 1)*sin(M_PI*x) };
49 }
50
51 // 0: velocity-x | 1: velocity-y | 2: pressure
52 571200 Dune::FieldVector<double, 3> stokes(const Dune::FieldVector<double, 2>& globalPos)
53 {
54 571200 Dune::FieldVector<double, 3> sol(0.0);
55 1142400 const double x = globalPos[0];
56 1142400 const double y = globalPos[1];
57
58 571200 using std::sin; using std::cos;
59 1142400 sol[0] = -cos(M_PI*x)*sin(M_PI*y);
60 1142400 sol[1] = sin(M_PI*x)*cos(M_PI*y);
61 1142400 sol[2] = 0.5*y*sin(M_PI*x);
62 571200 return sol;
63 }
64
65 // 0: mom-x | 1: mom-y | 2: mass
66 3625440 Dune::FieldVector<double, 3> stokesRHS(const Dune::FieldVector<double, 2>& globalPos)
67 {
68 3625440 const double x = globalPos[0];
69 3625440 const double y = globalPos[1];
70 3625440 Dune::FieldVector<double, 3> source(0.0);
71 3625440 using std::sin; using std::cos;
72 7250880 source[0] = 0.5*M_PI*y*cos(M_PI*x) - 2*M_PI*M_PI*cos(M_PI*x)*sin(M_PI*y);
73 7250880 source[1] = 0.5*sin(M_PI*x) + 2*M_PI*M_PI*sin(M_PI*x)*cos(M_PI*y);
74 3625440 return source;
75 }
76
77 // sigma_ff = -sym(grad(v)) + pI
78 78708 Dune::FieldMatrix<double, 2, 2> stokesStress(const Dune::FieldVector<double, 2>& globalPos)
79 {
80 78708 const double x = globalPos[0];
81 78708 const double y = globalPos[1];
82
83 78708 using std::exp; using std::sin; using std::cos;
84
85 78708 Dune::FieldMatrix<double, 2, 2> stress(0.0);
86 236124 stress[0][0] = (0.5*y - 2*M_PI*sin(M_PI*y))*sin(M_PI*x);
87 236124 stress[1][0] = 0.0;
88 393540 stress[0][1] = stress[1][0]; // symmetric
89 236124 stress[1][1] = (0.5*y - 2*M_PI*sin(M_PI*y))*sin(M_PI*x);
90 78708 return stress;
91 }
92
93 } // end namespace Rybak
94
95 /////////////////////////////////////////////////////////////////////////////////////////////////
96 // see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
97 // Example 1
98 /////////////////////////////////////////////////////////////////////////////////////////////////
99 namespace ShiueOne {
100
101 // 0: velocity-x | 1: velocity-y | 2: pressure
102 70560 Dune::FieldVector<double, 3> darcy(const Dune::FieldVector<double, 2>& globalPos)
103 {
104 70560 Dune::FieldVector<double, 3> sol(0.0);
105 141120 const double x = globalPos[0];
106 141120 const double y = globalPos[1];
107
108 70560 using std::exp; using std::sin; using std::cos;
109 141120 sol[0] = M_PI*(-exp(1)*y + exp(y))*sin(M_PI*x);
110 141120 sol[1] = (exp(1) - exp(y))*cos(M_PI*x);
111 141120 sol[2] = (exp(y) - y*exp(1)) * cos(M_PI*x);
112
113 70560 return sol;
114 }
115
116 // 0: mass
117 136640 Dune::FieldVector<double, 1> darcyRHS(const Dune::FieldVector<double, 2>& globalPos)
118 {
119 273280 const double x = globalPos[0];
120 273280 const double y = globalPos[1];
121 136640 using std::exp; using std::sin; using std::cos;
122 273280 return { cos(M_PI*x) * (M_PI*M_PI*(exp(y) - y*exp(1)) - exp(y)) };
123 }
124
125 // 0: velocity-x | 1: velocity-y | 2: pressure
126 577938 Dune::FieldVector<double, 3> stokes(const Dune::FieldVector<double, 2>& globalPos)
127 {
128 577938 Dune::FieldVector<double, 3> sol(0.0);
129 1155876 const double x = globalPos[0];
130 1155876 const double y = globalPos[1];
131
132 577938 using std::exp; using std::sin; using std::cos;
133 1155876 sol[0] = -1/M_PI * exp(y) * sin(M_PI*x);
134 1155876 sol[1] = (exp(y) - exp(1)) * cos(M_PI*x);
135 1155876 sol[2] = 2*exp(y) * cos(M_PI*x);
136 577938 return sol;
137 }
138
139 // 0: mom-x | 1: mom-y | 2: mass
140 3625440 Dune::FieldVector<double, 3> stokesRHS(const Dune::FieldVector<double, 2>& globalPos)
141 {
142 3625440 const double x = globalPos[0];
143 3625440 const double y = globalPos[1];
144 3625440 using std::exp; using std::sin; using std::cos;
145 3625440 Dune::FieldVector<double, 3> source(0.0);
146 7250880 source[0] = exp(y)*sin(M_PI*x) * (1/M_PI -3*M_PI);
147 7250880 source[1] = cos(M_PI*x) * (M_PI*M_PI*(exp(y)- exp(1)) + exp(y));
148 3625440 return source;
149 }
150
151 // sigma_ff = -sym(grad(v)) + pI
152 120 Dune::FieldMatrix<double, 2, 2> stokesStress(const Dune::FieldVector<double, 2>& globalPos)
153 {
154 120 const double x = globalPos[0];
155 120 const double y = globalPos[1];
156
157 120 using std::exp; using std::sin; using std::cos;
158
159 120 Dune::FieldMatrix<double, 2, 2> stress(0.0);
160 360 stress[0][0] = 4.0*exp(y)*cos(M_PI*x);
161 360 stress[1][0] = (M_PI*M_PI*(exp(y) - exp(1)) + 1.0*exp(y))*sin(M_PI*x)/M_PI;
162 600 stress[0][1] = stress[1][0]; // symmetric
163 360 stress[1][1] = 0.0;
164 120 return stress;
165 }
166
167 } // end namespace ShiueExOne
168
169 /////////////////////////////////////////////////////////////////////////////////////////////////
170 // see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
171 // Example 2
172 /////////////////////////////////////////////////////////////////////////////////////////////////
173 namespace ShiueTwo {
174
175 // 0: velocity-x | 1: velocity-y | 2: pressure
176 71144 Dune::FieldVector<double, 3> darcy(const Dune::FieldVector<double, 2>& globalPos)
177 {
178 71144 Dune::FieldVector<double, 3> sol(0.0);
179 142288 const double x = globalPos[0];
180 142288 const double y = globalPos[1];
181
182 142288 sol[0] = x*(y - 1.0) + (x - 1.0)*(y - 1.0) - 2.0;
183 142288 sol[1] = x*(x - 1.0) - 1.0*(y - 1.0)*(y - 1.0) - 2.0;
184 142288 sol[2] = x*(1.0-x)*(y-1.0) + (y-1.0)*(y-1.0)*(y-1.0)/3.0 + 2.0*x + 2.0*y + 4.0;
185
186 71144 return sol;
187 }
188
189 // 0: mass
190 ✗ Dune::FieldVector<double, 1> darcyRHS(const Dune::FieldVector<double, 2>& globalPos)
191 273280 { return { 0.0 }; }
192
193 // 0: velocity-x | 1: velocity-y | 2: pressure
194 575698 Dune::FieldVector<double, 3> stokes(const Dune::FieldVector<double, 2>& globalPos)
195 {
196 575698 Dune::FieldVector<double, 3> sol(0.0);
197 1151396 const double x = globalPos[0];
198 1151396 const double y = globalPos[1];
199
200 1151396 sol[0] = (y-1.0)*(y-1.0) + x*(y-1.0) + 3.0*x - 1.0;
201 1151396 sol[1] = x*(x-1.0) - 0.5*(y-1.0)*(y-1.0) - 3.0*y + 1.0;
202 1151396 sol[2] = 2.0*x + y - 1.0;
203 575698 return sol;
204 }
205
206 // 0: mom-x | 1: mom-y | 2: mass
207 ✗ Dune::FieldVector<double, 3> stokesRHS(const Dune::FieldVector<double, 2>& globalPos)
208 3625440 { return { 0.0, 0.0, 0.0 }; }
209
210 // sigma_ff = -sym(grad(v)) + pI
211 26332 Dune::FieldMatrix<double, 2, 2> stokesStress(const Dune::FieldVector<double, 2>& globalPos)
212 {
213 26332 const double x = globalPos[0];
214 26332 const double y = globalPos[1];
215
216 26332 Dune::FieldMatrix<double, 2, 2> stress(0.0);
217 78996 stress[0][0] = 2.0*x - y - 5.0;
218 78996 stress[1][0] = -3*x - 2*y + 3.0;
219 131660 stress[0][1] = stress[1][0]; // symmetric
220 78996 stress[1][1] = 2.0*x + 3.0*y + 3.0;
221 26332 return stress;
222 }
223
224 } // end namespace ShiueTwo
225
226 /////////////////////////////////////////////////////////////////////////////////////////////////
227 // see Schneider et al., 2019: "Coupling staggered-grid and MPFA finite volume methods for
228 // free flow/porous-medium flow problems (with c = 0)"
229 /////////////////////////////////////////////////////////////////////////////////////////////////
230 namespace Schneider {
231
232 // 0: velocity-x | 1: velocity-y | 2: pressure
233 74540 Dune::FieldVector<double, 3> darcy(const Dune::FieldVector<double, 2>& globalPos)
234 {
235 74540 Dune::FieldVector<double, 3> sol(0.0);
236 149080 const double x = globalPos[0];
237 149080 const double y = globalPos[1];
238 74540 static constexpr double omega = M_PI;
239 74540 static constexpr double c = 0.0;
240 74540 using std::exp; using std::sin; using std::cos;
241 74540 const double sinOmegaX = sin(omega*x);
242 74540 const double cosOmegaX = cos(omega*x);
243 74540 static const double expTwo = exp(2);
244 74540 const double expYPlusOne = exp(y+1);
245
246 223620 sol[0] = c/(2*omega)*expYPlusOne*sinOmegaX*sinOmegaX
247 74540 -omega*(expYPlusOne + 2 - expTwo)*cosOmegaX;
248 223620 sol[1] = (0.5*c*(expYPlusOne + 2 - expTwo)*cosOmegaX
249 74540 -(c*cosOmegaX + 1)*exp(y-1))*sinOmegaX;
250 149080 sol[2] = (expYPlusOne + 2 - expTwo)*sinOmegaX;
251
252 74540 return sol;
253 }
254
255 // 0: mass
256 341600 Dune::FieldVector<double, 1> darcyRHS(const Dune::FieldVector<double, 2>& globalPos)
257 {
258 683200 const double x = globalPos[0];
259 683200 const double y = globalPos[1];
260 341600 using std::exp; using std::sin; using std::cos;
261 341600 static constexpr double omega = M_PI;
262 341600 static constexpr double c = 0.0;
263 341600 const double cosOmegaX = cos(omega*x);
264 341600 static const double expTwo = exp(2);
265 341600 const double expYPlusOne = exp(y+1);
266
267 683200 return { sin(omega*x)*(-(c*cosOmegaX + 1)*exp(y - 1)
268 341600 + 1.5*c*expYPlusOne*cosOmegaX
269 1024800 + omega*omega*(expYPlusOne - expTwo + 2)) };
270 }
271
272 // 0: velocity-x | 1: velocity-y | 2: pressure
273 636186 Dune::FieldVector<double, 3> stokes(const Dune::FieldVector<double, 2>& globalPos)
274 {
275 636186 Dune::FieldVector<double, 3> sol(0.0);
276 1272372 const double x = globalPos[0];
277 1272372 const double y = globalPos[1];
278 636186 using std::sin;
279 636186 static constexpr double omega = M_PI;
280 636186 const double sinOmegaX = sin(omega*x);
281
282 1272372 sol[0] = y;
283 1272372 sol[1] = -y*sinOmegaX;
284 1272372 sol[2] = -y*y*sinOmegaX*sinOmegaX;
285 636186 return sol;
286 }
287
288 // 0: mom-x | 1: mom-y | 2: mass
289 9063600 Dune::FieldVector<double, 3> stokesRHS(const Dune::FieldVector<double, 2>& globalPos)
290 {
291 9063600 const double x = globalPos[0];
292 9063600 const double y = globalPos[1];
293 9063600 using std::exp; using std::sin; using std::cos;
294 9063600 static constexpr double omega = M_PI;
295 9063600 const double sinOmegaX = sin(omega*x);
296 9063600 const double cosOmegaX = cos(omega*x);
297
298 9063600 Dune::FieldVector<double, 3> source(0.0);
299 27190800 source[0] = -2*omega*y*y*sinOmegaX*cosOmegaX
300 9063600 -2*y*sinOmegaX + omega*cosOmegaX;
301 18127200 source[1] = -omega*y*y*cosOmegaX - omega*omega*y*sinOmegaX;
302 18127200 source[2] = -sinOmegaX;
303 9063600 return source;
304 }
305
306 // sigma_ff = vvT + sym(grad(v)) + pI
307 65830 Dune::FieldMatrix<double, 2, 2> stokesStress(const Dune::FieldVector<double, 2>& globalPos)
308 {
309 65830 using std::sin; using std::cos;
310 65830 const double x = globalPos[0];
311 65830 const double y = globalPos[1];
312 65830 static constexpr double omega = M_PI;
313 65830 const double sinOmegaX = sin(omega*x);
314 65830 const double cosOmegaX = cos(omega*x);
315
316 65830 Dune::FieldMatrix<double, 2, 2> stress(0.0);
317 197490 stress[0][0] = y*y * cosOmegaX*cosOmegaX;
318 197490 stress[0][1] = -y*y*sinOmegaX + omega*y*cosOmegaX - 1.0;
319 329150 stress[1][0] = stress[0][1]; // symmetric
320 197490 stress[1][1] = 2*sinOmegaX;
321 65830 return stress;
322 }
323
324 } // end namespace Schneider
325
326 } // end namespace Dumux::Solution::DarcyStokes
327
328 #endif
329