GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/problem_darcy.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 53 96 55.2%
Functions: 6 13 46.2%
Branches: 52 138 37.7%
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 The Darcy sub-problem of coupled Stokes-Darcy convergence test
11 */
12
13 #ifndef DUMUX_DARCY_SUBPROBLEM_HH
14 #define DUMUX_DARCY_SUBPROBLEM_HH
15
16 #include <dune/common/fvector.hh>
17 #include <dune/grid/yaspgrid.hh>
18
19 #include <dumux/discretization/cctpfa.hh>
20
21 #include <dumux/common/boundarytypes.hh>
22 #include <dumux/common/numeqvector.hh>
23
24 #include <dumux/material/components/constant.hh>
25 #include <dumux/material/fluidsystems/1pliquid.hh>
26
27 #include <dumux/porousmediumflow/1p/model.hh>
28 #include <dumux/porousmediumflow/problem.hh>
29
30 #include "spatialparams.hh"
31 #include "testcase.hh"
32
33 namespace Dumux {
34 template <class TypeTag>
35 class DarcySubProblem;
36
37 /*!
38 * \ingroup BoundaryTests
39 * \brief The Darcy sub-problem of coupled Stokes-Darcy convergence test
40 */
41 template <class TypeTag>
42 class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
43 {
44 using ParentType = PorousMediumFlowProblem<TypeTag>;
45 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
46 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
47 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
48 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
49 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
50 using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
51 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
52 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
53 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
54 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
55 using Element = typename GridView::template Codim<0>::Entity;
56 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
57
58 using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
59
60 static constexpr auto velocityXIdx = 0;
61 static constexpr auto velocityYIdx = 1;
62 static constexpr auto pressureIdx = 2;
63
64 public:
65 //! export the Indices
66 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
67
68 6 DarcySubProblem(std::shared_ptr<const GridGeometry> gridGeometry,
69 std::shared_ptr<CouplingManager> couplingManager,
70 std::shared_ptr<typename ParentType::SpatialParams> spatialParams,
71 const TestCase testCase,
72 const std::string& name)
73 : ParentType(gridGeometry, spatialParams, "Darcy")
74 , couplingManager_(couplingManager)
75
7/22
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 6 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 6 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 6 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 6 times.
✓ Branch 20 taken 6 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
 24 , testCase_(testCase)
76 {
77
7/20
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 6 times.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 6 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 6 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 6 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
 12 problemName_ = name + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
78 6 }
79
80 /*!
81 * \brief The problem name.
82 */
83 const std::string& name() const
84 {
85
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 12 return problemName_;
86 }
87
88 /*!
89 * \name Boundary conditions
90 */
91 // \{
92
93 /*!
94 * \brief Specifies which kind of boundary condition should be
95 * used for which equation on a given boundary control volume.
96 *
97 * \param element The element
98 * \param scvf The boundary sub control volume face
99 */
100 BoundaryTypes boundaryTypes(const Element &element, const SubControlVolumeFace &scvf) const
101 {
102
8/8
✓ Branch 0 taken 1120 times.
✓ Branch 1 taken 3360 times.
✓ Branch 2 taken 1120 times.
✓ Branch 3 taken 3360 times.
✓ Branch 4 taken 11760 times.
✓ Branch 5 taken 11928 times.
✓ Branch 6 taken 3080 times.
✓ Branch 7 taken 9240 times.
 44968 BoundaryTypes values;
103
104
8/8
✓ Branch 0 taken 1120 times.
✓ Branch 1 taken 3360 times.
✓ Branch 2 taken 1120 times.
✓ Branch 3 taken 3360 times.
✓ Branch 4 taken 11760 times.
✓ Branch 5 taken 11928 times.
✓ Branch 6 taken 3080 times.
✓ Branch 7 taken 9240 times.
 44968 if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
105 values.setAllCouplingNeumann();
106 else
107 values.setAllDirichlet();
108
109 return values;
110 }
111
112 /*!
113 * \brief Evaluates the boundary conditions for a Dirichlet control volume.
114 *
115 * \param element The element for which the Dirichlet boundary condition is set
116 * \param scvf The boundary subcontrolvolumeface
117 *
118 * For this method, the \a values parameter stores primary variables.
119 */
120 ✗ PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
121 {
122
0/4
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 18480 const auto p = analyticalSolution(scvf.center())[pressureIdx];
123 9240 return PrimaryVariables(p);
124 }
125
126 /*!
127 * \brief Evaluates the boundary conditions for a Neumann control volume.
128 *
129 * \param element The element for which the Neumann boundary condition is set
130 * \param fvGeometry The fvGeometry
131 * \param elemVolVars The element volume variables
132 * \param elemFluxVarsCache Flux variables caches for all faces in stencil
133 * \param scvf The boundary sub control volume face
134 *
135 * For this method, the \a values variable stores primary variables.
136 */
137 template<class ElementVolumeVariables, class ElementFluxVarsCache>
138 ✗ NumEqVector neumann(const Element& element,
139 const FVElementGeometry& fvGeometry,
140 const ElementVolumeVariables& elemVolVars,
141 const ElementFluxVarsCache& elemFluxVarsCache,
142 const SubControlVolumeFace& scvf) const
143 {
144 ✗ NumEqVector values(0.0);
145
146 ✗ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
147 ✗ values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf);
148
149 ✗ return values;
150 }
151
152 // \}
153
154 /*!
155 * \name Volume terms
156 */
157 // \{
158 /*!
159 * \brief Evaluates the source term for all phases within a given
160 * sub control volume.
161 *
162 * \param globalPos The global position
163 */
164 478240 NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
165 {
166
2/5
✗ Branch 0 not taken.
✓ Branch 1 taken 136640 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 341600 times.
✗ Branch 4 not taken.
 478240 switch (testCase_)
167 {
168 ✗ case TestCase::ShiueExampleOne:
169 ✗ return rhsShiueEtAlExampleOne_(globalPos);
170 case TestCase::ShiueExampleTwo:
171 136640 return rhsShiueEtAlExampleTwo_(globalPos);
172 ✗ case TestCase::Rybak:
173 ✗ return rhsRybak_(globalPos);
174 341600 case TestCase::Schneider:
175 341600 return rhsSchneiderEtAl_(globalPos);
176 ✗ default:
177 ✗ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
178 }
179 }
180
181 // \}
182
183 /*!
184 * \brief Evaluates the initial value for a control volume.
185 *
186 * \param element The element
187 *
188 * For this method, the \a priVars parameter stores primary
189 * variables.
190 */
191 PrimaryVariables initial(const Element &element) const
192 {
193 return PrimaryVariables(0.0);
194 }
195
196 /*!
197 * \brief Returns the analytical solution of the problem at a given position.
198 *
199 * \param globalPos The global position
200 */
201 143640 auto analyticalSolution(const GlobalPosition& globalPos) const
202 {
203
2/5
✗ Branch 0 not taken.
✓ Branch 1 taken 70560 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 73080 times.
✗ Branch 4 not taken.
 143640 switch (testCase_)
204 {
205 ✗ case TestCase::ShiueExampleOne:
206 ✗ return analyticalSolutionShiueEtAlExampleOne_(globalPos);
207 70560 case TestCase::ShiueExampleTwo:
208 70560 return analyticalSolutionShiueEtAlExampleTwo_(globalPos);
209 ✗ case TestCase::Rybak:
210 ✗ return analyticalSolutionRybak_(globalPos);
211 73080 case TestCase::Schneider:
212 73080 return analyticalSolutionSchneiderEtAl_(globalPos);
213 ✗ default:
214 ✗ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
215 }
216 }
217
218 // \}
219
220 //! Get the coupling manager
221 const CouplingManager& couplingManager() const
222
16/24
✓ Branch 0 taken 1120 times.
✓ Branch 1 taken 3360 times.
✓ Branch 2 taken 1120 times.
✓ Branch 3 taken 3360 times.
✓ Branch 4 taken 1120 times.
✓ Branch 5 taken 3360 times.
✓ Branch 6 taken 1120 times.
✓ Branch 7 taken 3360 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✓ Branch 18 taken 11760 times.
✓ Branch 19 taken 11928 times.
✓ Branch 20 taken 11760 times.
✓ Branch 21 taken 11928 times.
✓ Branch 22 taken 3080 times.
✓ Branch 23 taken 9240 times.
✓ Branch 24 taken 3080 times.
✓ Branch 25 taken 9240 times.
 89936 { return *couplingManager_; }
223
224 private:
225
226 // see Rybak et al., 2015: "Multirate time integration for coupled saturated/unsaturated porous medium and free flow systems"
227 ✗ Dune::FieldVector<Scalar, 3> analyticalSolutionRybak_(const GlobalPosition& globalPos) const
228 {
229 ✗ Dune::FieldVector<Scalar, 3> sol(0.0);
230 ✗ const Scalar x = globalPos[0];
231 ✗ const Scalar y = globalPos[1];
232
233 using std::exp; using std::sin; using std::cos;
234 ✗ sol[velocityXIdx] = -0.5*M_PI*y*y*cos(M_PI*x);
235 ✗ sol[velocityYIdx] = -1.0*y*sin(M_PI*x);
236 ✗ sol[pressureIdx] = 0.5*y*y*sin(M_PI*x);
237 ✗ return sol;
238 }
239
240 // see Rybak et al., 2015: "Multirate time integration for coupled saturated/unsaturated porous medium and free flow systems"
241 ✗ NumEqVector rhsRybak_(const GlobalPosition& globalPos) const
242 {
243 ✗ const Scalar x = globalPos[0];
244 ✗ const Scalar y = globalPos[1];
245 using std::sin;
246 ✗ return NumEqVector((0.5*M_PI*y*M_PI*y - 1)*sin(M_PI*x));
247 }
248
249 // see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
250 ✗ Dune::FieldVector<Scalar, 3> analyticalSolutionShiueEtAlExampleOne_(const GlobalPosition& globalPos) const
251 {
252 ✗ Dune::FieldVector<Scalar, 3> sol(0.0);
253 ✗ const Scalar x = globalPos[0];
254 ✗ const Scalar y = globalPos[1];
255
256 using std::exp; using std::sin; using std::cos;
257 ✗ sol[pressureIdx] = (exp(y) - y*exp(1)) * cos(M_PI*x);
258 ✗ sol[velocityXIdx] = M_PI*(-exp(1)*y + exp(y))*sin(M_PI*x);
259 ✗ sol[velocityYIdx] = (exp(1) - exp(y))*cos(M_PI*x);
260
261 ✗ return sol;
262 }
263
264 // see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
265 ✗ NumEqVector rhsShiueEtAlExampleOne_(const GlobalPosition& globalPos) const
266 {
267 ✗ const Scalar x = globalPos[0];
268 ✗ const Scalar y = globalPos[1];
269 using std::exp; using std::sin; using std::cos;
270 ✗ return NumEqVector(cos(M_PI*x) * (M_PI*M_PI*(exp(y) - y*exp(1)) - exp(y)));
271 }
272
273 // see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
274 70560 Dune::FieldVector<Scalar, 3> analyticalSolutionShiueEtAlExampleTwo_(const GlobalPosition& globalPos) const
275 {
276 70560 Dune::FieldVector<Scalar, 3> sol(0.0);
277 141120 const Scalar x = globalPos[0];
278 141120 const Scalar y = globalPos[1];
279
280 141120 sol[pressureIdx] = 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;
281 141120 sol[velocityXIdx] = x*(y - 1.0) + (x - 1.0)*(y - 1.0) - 2.0;
282 141120 sol[velocityYIdx] = x*(x - 1.0) - 1.0*(y - 1.0)*(y - 1.0) - 2.0;
283
284 70560 return sol;
285 }
286
287 // see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
288 ✗ NumEqVector rhsShiueEtAlExampleTwo_(const GlobalPosition& globalPos) const
289 136640 { return NumEqVector(0.0); }
290
291 // see Schneider et al., 2019: "Coupling staggered-grid and MPFA finite volume methods for
292 // free flow/porous-medium flow problems"
293 73080 Dune::FieldVector<Scalar, 3> analyticalSolutionSchneiderEtAl_(const GlobalPosition& globalPos) const
294 {
295 73080 Dune::FieldVector<Scalar, 3> sol(0.0);
296 146160 const Scalar x = globalPos[0];
297 146160 const Scalar y = globalPos[1];
298 static constexpr Scalar omega = M_PI;
299 static constexpr Scalar c = 0.0;
300 using std::exp; using std::sin; using std::cos;
301 73080 const Scalar sinOmegaX = sin(omega*x);
302 73080 const Scalar cosOmegaX = cos(omega*x);
303 static const Scalar expTwo = exp(2);
304 73080 const Scalar expYPlusOne = exp(y+1);
305
306 146160 sol[pressureIdx] = (expYPlusOne + 2 - expTwo)*sinOmegaX;
307 219240 sol[velocityXIdx] = c/(2*omega)*expYPlusOne*sinOmegaX*sinOmegaX
308 73080 -omega*(expYPlusOne + 2 - expTwo)*cosOmegaX;
309 219240 sol[velocityYIdx] = (0.5*c*(expYPlusOne + 2 - expTwo)*cosOmegaX
310 73080 -(c*cosOmegaX + 1)*exp(y-1))*sinOmegaX;
311
312 73080 return sol;
313 }
314
315 // see Schneider et al., 2019: "Coupling staggered-grid and MPFA finite volume methods for
316 // free flow/porous-medium flow problems (with c = 0)"
317 341600 NumEqVector rhsSchneiderEtAl_(const GlobalPosition& globalPos) const
318 {
319 683200 const Scalar x = globalPos[0];
320 683200 const Scalar y = globalPos[1];
321 using std::exp; using std::sin; using std::cos;
322 static constexpr Scalar omega = M_PI;
323 static constexpr Scalar c = 0.0;
324 341600 const Scalar cosOmegaX = cos(omega*x);
325 static const Scalar expTwo = exp(2);
326 341600 const Scalar expYPlusOne = exp(y+1);
327
328 683200 const Scalar result = (-(c*cosOmegaX + 1)*exp(y - 1)
329 341600 + 1.5*c*expYPlusOne*cosOmegaX
330 341600 + omega*omega*(expYPlusOne - expTwo + 2))
331 341600 *sin(omega*x);
332 341600 return NumEqVector(result);
333 }
334
335 static constexpr Scalar eps_ = 1e-7;
336 std::shared_ptr<CouplingManager> couplingManager_;
337 std::string problemName_;
338 TestCase testCase_;
339 };
340 } // end namespace Dumux
341
342 #endif
343