GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_freeflow.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 81 94 86.2%
Functions: 11 15 73.3%
Branches: 94 219 42.9%
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 free-flow sub-problem of coupled FreeFlow/Darcy convergence test
11 */
12
13 #ifndef DUMUX_FREEFLOW_SUBPROBLEM_HH
14 #define DUMUX_FREEFLOW_SUBPROBLEM_HH
15
16 #include <dumux/common/properties.hh>
17 #include <dumux/freeflow/navierstokes/momentum/fluxhelper.hh>
18 #include <dumux/freeflow/navierstokes/scalarfluxhelper.hh>
19 #include <dumux/freeflow/navierstokes/mass/1p/advectiveflux.hh>
20 #include <dune/common/fmatrix.hh>
21
22 #include "testcase.hh"
23 #include "analyticalsolutions.hh"
24
25 namespace Dumux {
26
27 /*!
28 * \ingroup BoundaryTests
29 * \brief The Stokes sub-problem of coupled Stokes-Darcy convergence test
30 */
31 template <class TypeTag, class BaseProblem>
32 class FreeFlowSubProblem : public BaseProblem
33 {
34 using ParentType = BaseProblem;
35 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
36 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
37 using GridView = typename GridGeometry::GridView;
38 using FVElementGeometry = typename GridGeometry::LocalView;
39 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
40 using Element = typename GridView::template Codim<0>::Entity;
41 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
42 using InitialValues = typename ParentType::InitialValues;
43 using Sources = typename ParentType::Sources;
44 using DirichletValues = typename ParentType::DirichletValues;
45 using BoundaryFluxes = typename ParentType::BoundaryFluxes;
46 using BoundaryTypes = typename ParentType::BoundaryTypes;
47 using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
48
49 using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
50
51 enum class BC {
52 dirichlet, stress, mixed
53 };
54
55 public:
56 //! export the Indices
57 using Indices = typename ModelTraits::Indices;
58
59 48 FreeFlowSubProblem(std::shared_ptr<const GridGeometry> gridGeometry,
60 std::shared_ptr<CouplingManager> couplingManager,
61 const DarcyStokesTestCase testCase,
62 const std::string& name)
63 : ParentType(gridGeometry, couplingManager, "FreeFlow")
64 , couplingManager_(couplingManager)
65
6/20
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 24 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 24 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 24 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 24 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
 192 , testCase_(testCase)
66 {
67
5/16
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 24 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 24 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 24 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
 48 problemName_ = name + "_"
68
2/4
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
 96 + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
69
70
1/4
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 48 auto bc = getParamFromGroup<std::string>(this->paramGroup(), "Problem.BoundaryConditions", "Dirichlet");
71
2/2
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 18 times.
 48 if (bc == "Dirichlet")
72 12 boundaryConditions_ = BC::dirichlet;
73
2/2
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 12 times.
 36 else if (bc == "Stress")
74 12 boundaryConditions_ = BC::stress;
75
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 24 else if (bc == "Mixed")
76 24 boundaryConditions_ = BC::mixed;
77 else
78 ✗ DUNE_THROW(Dune::Exception, "Wrong BC type choose: Dirichlet, Stress or Mixed");
79
80
3/6
✓ Branch 2 taken 24 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 24 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 24 times.
 144 std::cout << "Free flow domain: Using " << bc << " boundary conditions" << std::endl;
81 48 }
82
83 /*!
84 * \name Problem parameters
85 */
86 // \{
87
88 const std::string& name() const
89
2/4
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
 24 { return problemName_; }
90
91 // \}
92
93 /*!
94 * \name Boundary conditions
95 */
96 // \{
97
98 /*!
99 * \brief Specifies which kind of boundary condition should be
100 * used for which equation on a given boundary segment.
101 *
102 * \param element The finite element
103 * \param scvf The sub control volume face
104 */
105 141516 BoundaryTypes boundaryTypes(const Element& element,
106 const SubControlVolumeFace& scvf) const
107 {
108
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 141516 BoundaryTypes values;
109
110 if constexpr (ParentType::isMomentumProblem())
111 {
112
2/2
✓ Branch 0 taken 17685 times.
✓ Branch 1 taken 53073 times.
 141516 if (couplingManager().isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf))
113 {
114 35370 values.setCouplingNeumann(Indices::momentumYBalanceIdx);
115 35370 values.setCouplingNeumann(Indices::momentumXBalanceIdx);
116 }
117 else
118 {
119
2/2
✓ Branch 0 taken 10098 times.
✓ Branch 1 taken 42975 times.
 106146 if (boundaryConditions_ == BC::dirichlet)
120 values.setAllDirichlet();
121
2/2
✓ Branch 0 taken 10134 times.
✓ Branch 1 taken 32841 times.
 85950 else if (boundaryConditions_ == BC::stress)
122 values.setAllNeumann();
123 else
124 {
125
2/2
✓ Branch 0 taken 21894 times.
✓ Branch 1 taken 10947 times.
 65682 if (onLeftBoundary_(scvf.center()))
126 values.setAllNeumann();
127 else
128 values.setAllDirichlet();
129 }
130 }
131 }
132 else
133 {
134 ✗ if (couplingManager().isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, scvf))
135 values.setAllCouplingNeumann();
136 else
137 values.setAllNeumann();
138 }
139
140 141516 return values;
141 }
142
143 /*!
144 * \brief Evaluates the boundary conditions for a Dirichlet control volume.
145 */
146 DirichletValues dirichletAtPos(const GlobalPosition& globalPos) const
147 76222 { return analyticalSolution(globalPos); }
148
149 /*!
150 * \brief Evaluates the boundary conditions for a Neumann control volume.
151 *
152 * \param element The element for which the Neumann boundary condition is set
153 * \param fvGeometry The fvGeometry
154 * \param elemVolVars The element volume variables
155 * \param elemFaceVars The element face variables
156 * \param scvf The boundary sub control volume face
157 */
158 template<class ElementVolumeVariables, class ElementFluxVariablesCache>
159 948784 BoundaryFluxes neumann(const Element& element,
160 const FVElementGeometry& fvGeometry,
161 const ElementVolumeVariables& elemVolVars,
162 const ElementFluxVariablesCache& elemFluxVarsCache,
163 const SubControlVolumeFace& scvf) const
164 {
165 948784 BoundaryFluxes values(0.0);
166 948784 const auto& globalPos = scvf.ipGlobal();
167 using SlipVelocityPolicy = NavierStokesSlipVelocity<typename GridGeometry::DiscretizationMethod, NavierStokes::SlipConditions::BJS>;
168 using FluxHelper = NavierStokesMomentumBoundaryFlux<typename GridGeometry::DiscretizationMethod, SlipVelocityPolicy>;
169
170 // momentum boundary conditions
171 if constexpr (ParentType::isMomentumProblem())
172 {
173 // We need to take care here: If the integration point of an scvf lies both on the coupling interface and the
174 // domain boundary (i.e., the leftmost and rightmost points of the interface), do not evaluate the coupling or slip condition
175 // because they would require data from the boundary not available in this case (velocities for evaluating gradients,
176 // those would only be available for Dirichlet BCs). Instead, directly use the given Neumann boundary stress.
177 // TODO: Maybe couplingManager().isCoupled(...) could return false for these scvfs.
178 //
179 // | <-- left Neumann boundary
180 // |
181 // integration point --> o##### <-- scvf lying on coupling interface with integration point touching left Neumann boundary
182 //
183 1379600 if (couplingManager().isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf)
184
12/12
✓ Branch 0 taken 173970 times.
✓ Branch 1 taken 270 times.
✓ Branch 2 taken 173970 times.
✓ Branch 3 taken 270 times.
✓ Branch 4 taken 173970 times.
✓ Branch 5 taken 270 times.
✓ Branch 6 taken 173970 times.
✓ Branch 7 taken 270 times.
✓ Branch 8 taken 173970 times.
✓ Branch 9 taken 270 times.
✓ Branch 10 taken 173970 times.
✓ Branch 11 taken 270 times.
 2090880 && scvf.ipGlobal()[0] > this->gridGeometry().bBoxMin()[0] + eps_
185
14/14
✓ Branch 0 taken 174240 times.
✓ Branch 1 taken 170660 times.
✓ Branch 2 taken 173910 times.
✓ Branch 3 taken 60 times.
✓ Branch 4 taken 173910 times.
✓ Branch 5 taken 60 times.
✓ Branch 6 taken 173910 times.
✓ Branch 7 taken 60 times.
✓ Branch 8 taken 173910 times.
✓ Branch 9 taken 60 times.
✓ Branch 10 taken 173910 times.
✓ Branch 11 taken 60 times.
✓ Branch 12 taken 173910 times.
✓ Branch 13 taken 60 times.
 2777440 && scvf.ipGlobal()[0] < this->gridGeometry().bBoxMax()[0] - eps_)
186 {
187 1043460 values += couplingManager().momentumCouplingCondition(
188 CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex,
189 fvGeometry, scvf, elemVolVars
190 );
191
192 695640 values += FluxHelper::slipVelocityMomentumFlux(
193 *this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache
194 );
195 }
196 else
197 {
198 341980 const auto stress = stressTensorAtPos(globalPos);
199 341980 const auto n = scvf.unitOuterNormal();
200 683960 stress.mv(n, values);
201 }
202 }
203
204 // mass boundary conditions
205 else
206 {
207
2/2
✓ Branch 0 taken 36960 times.
✓ Branch 1 taken 92532 times.
 258984 if (couplingManager().isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, scvf))
208 {
209 147840 values = couplingManager().massCouplingCondition(
210 CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex,
211 fvGeometry, scvf, elemVolVars
212 );
213 }
214 else
215 {
216 using FluxHelper = NavierStokesScalarBoundaryFluxHelper<AdvectiveFlux<ModelTraits>>;
217 185064 values = FluxHelper::scalarOutflowFlux(
218 *this, element, fvGeometry, scvf, elemVolVars
219 );
220 }
221 }
222
223 948784 return values;
224 }
225
226 /*!
227 * \brief Returns true if the scvf lies on a porous slip boundary
228 */
229 ✗ bool onSlipBoundary(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
230 {
231 return scvf.boundary()
232
4/12
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 73128 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 73128 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 66968 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 66968 times.
✗ Branch 11 not taken.
 280192 && couplingManager().isCoupled(
233 CouplingManager::freeFlowMomentumIndex,
234 CouplingManager::porousMediumIndex,
235 scvf
236 ✗ );
237 }
238
239 // \}
240
241 //! Get the coupling manager
242 const CouplingManager& couplingManager() const
243
0/4
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
 1792232 { return *couplingManager_; }
244
245 /*!
246 * \name Volume terms
247 */
248 // \{
249
250 /*!
251 * \brief Returns the sources within the domain.
252 * \param globalPos The global position
253 */
254 39879840 Sources sourceAtPos(const GlobalPosition &globalPos) const
255 {
256 ✗ const auto select = [&](const auto& rhs) -> Sources {
257 if constexpr (ParentType::isMomentumProblem())
258 64951040 return { rhs[0], rhs[1] };
259 else
260 7404320 return { rhs[2] };
261 };
262
263 using namespace Solution::DarcyStokes;
264
4/5
✓ Branch 0 taken 3625440 times.
✓ Branch 1 taken 3625440 times.
✓ Branch 2 taken 3625440 times.
✓ Branch 3 taken 9063600 times.
✗ Branch 4 not taken.
 39879840 switch (testCase_)
265 {
266 7250880 case DarcyStokesTestCase::ShiueExampleOne:
267 7250880 return select(ShiueOne::stokesRHS(globalPos));
268 7250880 case DarcyStokesTestCase::ShiueExampleTwo:
269 7250880 return select(ShiueTwo::stokesRHS(globalPos));
270 7250880 case DarcyStokesTestCase::Rybak:
271 7250880 return select(Rybak::stokesRHS(globalPos));
272 18127200 case DarcyStokesTestCase::Schneider:
273 18127200 return select(Schneider::stokesRHS(globalPos));
274 ✗ default:
275 ✗ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
276 }
277 }
278
279 /*!
280 * \brief Evaluates the initial value for a control volume.
281 * \param globalPos The global position
282 */
283 InitialValues initialAtPos(const GlobalPosition& globalPos) const
284 { return InitialValues(0.0); }
285
286 /*!
287 * \brief Returns the beta value which is the alpha value divided by the square root of the (scalar-valued) interface permeability.
288 */
289 140096 Scalar betaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf, const GlobalPosition& tangentialVector) const
290 {
291 280192 const auto& K = couplingManager().darcyPermeability(fvGeometry, scvf);
292
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 140096 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 140096 times.
 280192 const auto alphaBJ = couplingManager().problem(CouplingManager::porousMediumIndex).spatialParams().beaversJosephCoeffAtPos(scvf.ipGlobal());
293 140096 const auto interfacePermeability = vtmv(tangentialVector, K, tangentialVector);
294 using std::sqrt;
295 140096 return alphaBJ / sqrt(interfacePermeability);
296 }
297
298 /*!
299 * \brief Returns the analytical solution of the problem at a given position.
300 * \param globalPos The global position
301 * Returns vector with entries: (velocity-x | velocity-y | pressure)
302 */
303 4722044 Dune::FieldVector<Scalar, 3> fullAnalyticalSolution(const GlobalPosition& globalPos) const
304 {
305 using namespace Solution::DarcyStokes;
306
4/5
✓ Branch 0 taken 577938 times.
✓ Branch 1 taken 575698 times.
✓ Branch 2 taken 571200 times.
✓ Branch 3 taken 636186 times.
✗ Branch 4 not taken.
 4722044 switch (testCase_)
307 {
308 1155876 case DarcyStokesTestCase::ShiueExampleOne:
309 1155876 return ShiueOne::stokes(globalPos);
310 1151396 case DarcyStokesTestCase::ShiueExampleTwo:
311 1151396 return ShiueTwo::stokes(globalPos);
312 1142400 case DarcyStokesTestCase::Rybak:
313 1142400 return Rybak::stokes(globalPos);
314 1272372 case DarcyStokesTestCase::Schneider:
315 1272372 return Schneider::stokes(globalPos);
316 ✗ default:
317 ✗ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
318 }
319 }
320
321 /*!
322 * \brief Returns the analytical solution of the problem at a given position.
323 * \param globalPos The global position
324 * Returns analytical solution depending on the type of sub-problem
325 */
326 ✗ DirichletValues analyticalSolution(const GlobalPosition& globalPos, Scalar time = 0.0) const
327 {
328 using namespace Solution::DarcyStokes;
329
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 403200 const auto sol = fullAnalyticalSolution(globalPos);
330 if constexpr (ParentType::isMomentumProblem())
331 ✗ return { sol[0], sol[1] };
332 else
333
6/6
✓ Branch 0 taken 3354 times.
✓ Branch 1 taken 265446 times.
✓ Branch 2 taken 3354 times.
✓ Branch 3 taken 265446 times.
✓ Branch 4 taken 3354 times.
✓ Branch 5 taken 265446 times.
 1209600 return { sol[2] };
334 }
335
336 // \}
337
338 private:
339 /*!
340 * \brief Returns the stress tensor within the domain.
341 * \param globalPos The global position
342 */
343 170990 Dune::FieldMatrix<double, 2, 2> stressTensorAtPos(const GlobalPosition &globalPos) const
344 {
345 using namespace Solution::DarcyStokes;
346
4/5
✓ Branch 0 taken 120 times.
✓ Branch 1 taken 26332 times.
✓ Branch 2 taken 78708 times.
✓ Branch 3 taken 65830 times.
✗ Branch 4 not taken.
 170990 switch (testCase_)
347 {
348 120 case DarcyStokesTestCase::ShiueExampleOne:
349 120 return ShiueOne::stokesStress(globalPos);
350 26332 case DarcyStokesTestCase::ShiueExampleTwo:
351 26332 return ShiueTwo::stokesStress(globalPos);
352 78708 case DarcyStokesTestCase::Rybak:
353 78708 return Rybak::stokesStress(globalPos);
354 65830 case DarcyStokesTestCase::Schneider:
355 65830 return Schneider::stokesStress(globalPos);
356 ✗ default:
357 ✗ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
358 }
359 }
360
361 bool onLeftBoundary_(const GlobalPosition& globalPos) const
362
10/10
✓ Branch 0 taken 21894 times.
✓ Branch 1 taken 10947 times.
✓ Branch 2 taken 21894 times.
✓ Branch 3 taken 10947 times.
✓ Branch 4 taken 21894 times.
✓ Branch 5 taken 10947 times.
✓ Branch 6 taken 21894 times.
✓ Branch 7 taken 10947 times.
✓ Branch 8 taken 21894 times.
✓ Branch 9 taken 10947 times.
 164205 { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
363
364 static constexpr Scalar eps_ = 1e-7;
365 std::string problemName_;
366 std::shared_ptr<CouplingManager> couplingManager_;
367 DarcyStokesTestCase testCase_;
368 BC boundaryConditions_;
369 };
370 } // end namespace Dumux
371
372 #endif
373