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 OnePTests | ||
10 | * \brief The problem setup for the convergence test with analytic solution | ||
11 | */ | ||
12 | #ifndef DUMUX_CONVERGENCE_TEST_ONEP_PROBLEM_HH | ||
13 | #define DUMUX_CONVERGENCE_TEST_ONEP_PROBLEM_HH | ||
14 | |||
15 | #include <cmath> | ||
16 | #include <dune/common/fvector.hh> | ||
17 | |||
18 | #include <dumux/common/properties.hh> | ||
19 | #include <dumux/common/boundarytypes.hh> | ||
20 | #include <dumux/common/numeqvector.hh> | ||
21 | #include <dumux/porousmediumflow/problem.hh> | ||
22 | |||
23 | namespace Dumux { | ||
24 | |||
25 | /*! | ||
26 | * \ingroup OnePTests | ||
27 | * \brief The problem setup for the convergence test with analytic solution | ||
28 | */ | ||
29 | template <class TypeTag> | ||
30 | 20 | class ConvergenceProblem : public PorousMediumFlowProblem<TypeTag> | |
31 | { | ||
32 | using ParentType = PorousMediumFlowProblem<TypeTag>; | ||
33 | using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; | ||
34 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
35 | using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; | ||
36 | using NumEqVector = Dumux::NumEqVector<PrimaryVariables>; | ||
37 | using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; | ||
38 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
39 | using Element = typename GridView::template Codim<0>::Entity; | ||
40 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
41 | |||
42 | static constexpr auto velocityXIdx = 0; | ||
43 | static constexpr auto velocityYIdx = 1; | ||
44 | static constexpr auto pressureIdx = 2; | ||
45 | |||
46 | public: | ||
47 | /*! | ||
48 | * \brief The constructor. | ||
49 | * \param gridGeometry The finite-volume grid geometry | ||
50 | */ | ||
51 | 20 | ConvergenceProblem(std::shared_ptr<const GridGeometry> gridGeometry) | |
52 | : ParentType(gridGeometry) | ||
53 |
7/18✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 20 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 20 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 20 times.
✓ Branch 15 taken 20 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 20 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
|
60 | , c_(getParam<Scalar>("Problem.C")) |
54 | 20 | {} | |
55 | |||
56 | /*! | ||
57 | * \name Problem parameters | ||
58 | */ | ||
59 | // \{ | ||
60 | |||
61 | /*! | ||
62 | * \name Boundary conditions | ||
63 | */ | ||
64 | // \{ | ||
65 | |||
66 | /*! | ||
67 | * \brief Specifies which kind of boundary condition should be | ||
68 | * used for which equation on a given boundary control volume. | ||
69 | * | ||
70 | * \param globalPos The position of the center of the finite volume | ||
71 | */ | ||
72 | ✗ | BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const | |
73 | { | ||
74 |
6/12✓ Branch 0 taken 3144 times.
✓ Branch 1 taken 3240 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 29136 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 33544 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 15120 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 12768 times.
✗ Branch 11 not taken.
|
96952 | BoundaryTypes values; |
75 |
6/12✓ Branch 0 taken 3144 times.
✓ Branch 1 taken 3240 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 29136 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 33544 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 15120 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 12768 times.
✗ Branch 11 not taken.
|
96952 | values.setAllDirichlet(); |
76 | ✗ | return values; | |
77 | } | ||
78 | |||
79 | /*! | ||
80 | * \brief Evaluates Dirichlet boundary conditions. | ||
81 | * \param globalPos The center of the finite volume which ought to be set. | ||
82 | */ | ||
83 | PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const | ||
84 | { | ||
85 |
1/2✓ Branch 2 taken 29344 times.
✗ Branch 3 not taken.
|
55048 | const auto p = analyticalSolution(globalPos)[pressureIdx]; |
86 |
1/2✓ Branch 1 taken 29344 times.
✗ Branch 2 not taken.
|
55048 | return PrimaryVariables(p); |
87 | } | ||
88 | |||
89 | // \} | ||
90 | |||
91 | /*! | ||
92 | * \name Volume terms | ||
93 | */ | ||
94 | // \{ | ||
95 | |||
96 | /*! | ||
97 | * \brief Evaluates the source term for all phases within a given | ||
98 | * sub-control volume. | ||
99 | * | ||
100 | * For this method, the \a priVars parameter stores the rate mass | ||
101 | * of a component is generated or annihilated per volume | ||
102 | * unit. Positive values mean that mass is created, negative ones | ||
103 | * mean that it vanishes. | ||
104 | * | ||
105 | * The units must be according to either using mole or mass fractions (mole/(m^3*s) or kg/(m^3*s)). | ||
106 | */ | ||
107 | 971040 | NumEqVector sourceAtPos(const GlobalPosition& globalPos) const | |
108 | { | ||
109 | 1942080 | const Scalar x = globalPos[0]; | |
110 | 1942080 | const Scalar y = globalPos[1]; | |
111 | using std::exp; | ||
112 | using std::sin; | ||
113 | using std::cos; | ||
114 | 971040 | const Scalar cosOmegaX = cos(omega_*x); | |
115 | static const Scalar expTwo = exp(2); | ||
116 | 971040 | const Scalar expYPlusOne = exp(y+1); | |
117 | |||
118 | 1942080 | const Scalar result = ( -(c_*cosOmegaX + 1)*exp(y - 1) | |
119 | 971040 | + 1.5*c_*expYPlusOne*cosOmegaX | |
120 | 971040 | + omega_*omega_*(expYPlusOne - expTwo + 2)) | |
121 | 971040 | * sin(omega_*x); | |
122 | |||
123 | 971040 | return NumEqVector(result); | |
124 | } | ||
125 | |||
126 | // \} | ||
127 | |||
128 | /*! | ||
129 | * \brief Evaluates the initial value for a control volume. | ||
130 | * \param globalPos The position for which the initial condition should be evaluated | ||
131 | */ | ||
132 | PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const | ||
133 | { return PrimaryVariables(0.0); } | ||
134 | |||
135 | /*! | ||
136 | * \brief Returns the analytical solution of the problem at a given position. | ||
137 | * \param globalPos The global position | ||
138 | */ | ||
139 | 133248 | auto analyticalSolution(const GlobalPosition& globalPos) const | |
140 | { | ||
141 | 133248 | Dune::FieldVector<Scalar, 3> sol(0.0); | |
142 | 266496 | const Scalar x = globalPos[0]; | |
143 | 266496 | const Scalar y = globalPos[1]; | |
144 | using std::exp; using std::sin; using std::cos; | ||
145 | 133248 | const Scalar sinOmegaX = sin(omega_*x); | |
146 | 133248 | const Scalar cosOmegaX = cos(omega_*x); | |
147 | static const Scalar expTwo = exp(2); | ||
148 | 133248 | const Scalar expYPlusOne = exp(y+1); | |
149 | |||
150 | 266496 | sol[pressureIdx] = (expYPlusOne + 2 - expTwo)*sinOmegaX + 10.0; | |
151 | 399744 | sol[velocityXIdx] = c_/(2*omega_)*expYPlusOne*sinOmegaX*sinOmegaX | |
152 | 133248 | -omega_*(expYPlusOne + 2 - expTwo)*cosOmegaX; | |
153 | 399744 | sol[velocityYIdx] = (0.5*c_*(expYPlusOne + 2 - expTwo)*cosOmegaX | |
154 | 133248 | -(c_*cosOmegaX + 1)*exp(y-1))*sinOmegaX; | |
155 | |||
156 | 133248 | return sol; | |
157 | } | ||
158 | |||
159 | private: | ||
160 | static constexpr Scalar eps_ = 1e-7; | ||
161 | static constexpr Scalar omega_ = M_PI; | ||
162 | Scalar c_; | ||
163 | }; | ||
164 | |||
165 | } // end namespace Dumux | ||
166 | |||
167 | #endif | ||
168 |