GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/2p2c/chemicalnonequilibrium/problem.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 63 70 90.0%
Functions: 10 16 62.5%
Branches: 68 124 54.8%
Line Branch Exec Source
1 //
2 // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
3 // SPDX-License-Identifier: GPL-3.0-or-later
4 //
5
6 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
7 // vi: set et ts=4 sw=4 sts=4:
8 /*!
9 * \file
10 * \ingroup TwoPTwoCTests
11 * \brief Problem for the 2p2c chemical nonequilibrium problem.
12 */
13
14 #ifndef DUMUX_TWOPTWOC_NONEQUILIBRIUM_PROBLEM_HH
15 #define DUMUX_TWOPTWOC_NONEQUILIBRIUM_PROBLEM_HH
16
17 #include <dumux/common/properties.hh>
18 #include <dumux/common/parameters.hh>
19 #include <dumux/common/boundarytypes.hh>
20 #include <dumux/common/numeqvector.hh>
21
22 #include <dumux/porousmediumflow/problem.hh>
23
24 namespace Dumux {
25
26 /*!
27 * \ingroup TwoPTwoCTests
28 * \brief Problem for the 2p2c chemical nonequilibrium problem.
29 */
30 template <class TypeTag>
31 class TwoPTwoCChemicalNonequilibriumProblem : public PorousMediumFlowProblem<TypeTag>
32 {
33 using ParentType = PorousMediumFlowProblem<TypeTag>;
34 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
35 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
36 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
37 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
38 using NeumannFluxes = Dumux::NumEqVector<PrimaryVariables>;
39 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
40 using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
41 using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
42 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
43 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
44 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
45 using Element = typename GridView::template Codim<0>::Entity;
46 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
47 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
48 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
49 using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
50
51 using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
52 using Indices = typename ModelTraits::Indices;
53
54 static constexpr int dim = GridView::dimension;
55 static constexpr int dimWorld = GridView::dimensionworld;
56 static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethods::box;
57 enum { dofCodim = isBox ? dim : 0 };
58 public:
59 2 TwoPTwoCChemicalNonequilibriumProblem(std::shared_ptr<const GridGeometry> gridGeometry)
60
10/32
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 2 times.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✓ Branch 24 taken 2 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 2 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
 6 : ParentType(gridGeometry)
61 {
62 // initialize the tables of the fluid system
63
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 Scalar Tmin = this->spatialParams().temperatureAtPos(GlobalPosition()) - 1.0;
64
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 Scalar Tmax = this->spatialParams().temperatureAtPos(GlobalPosition()) + 1.0;
65 2 int nT = 3;
66
67 2 Scalar pmin = 1.0e5 * 0.75;
68 2 Scalar pmax = 2.0e5 * 1.25;
69 2 int np = 1000;
70
71
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np);
72
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 2 name_ = getParam<std::string>("Problem.Name");
73 2 }
74
75 void setGridVariables(std::shared_ptr<GridVariables> gridVariables)
76
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 { gridVariables_ = gridVariables; }
77
78 const GridVariables& gridVariables() const
79 1591072 { return *gridVariables_; }
80
81
82 /*!
83 * \brief Returns the problem name
84 *
85 * This is used as a prefix for files generated by the simulation.
86 */
87 const std::string name() const
88
2/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2 { return name_; }
89
90 /*!
91 * \name Boundary conditions
92 * \brief Specifies which kind of boundary condition should be
93 * used for which equation on a given boundary segment
94 *
95 * \param globalPos The global position
96 */
97 32487 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
98 {
99 32487 BoundaryTypes bcTypes;
100
10/10
✓ Branch 0 taken 26747 times.
✓ Branch 1 taken 5740 times.
✓ Branch 2 taken 26747 times.
✓ Branch 3 taken 5740 times.
✓ Branch 4 taken 26747 times.
✓ Branch 5 taken 5740 times.
✓ Branch 6 taken 26747 times.
✓ Branch 7 taken 5740 times.
✓ Branch 8 taken 26747 times.
✓ Branch 9 taken 5740 times.
 162435 if (globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_)
101 bcTypes.setAllDirichlet();
102 else
103 bcTypes.setAllNeumann();
104 32487 return bcTypes;
105 }
106
107 /*!
108 * \brief Evaluates the boundary conditions for a Dirichlet oundary segment.
109 *
110 * \param globalPos The global position
111 */
112 ✗ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
113 {
114 ✗ PrimaryVariables values;
115 ✗ values = initial_(globalPos);
116 ✗ return values;
117 }
118
119 /*!
120 * \brief Evaluates the boundary conditions for a Neumann boundary segment.
121 *
122 * \param element The finite element
123 * \param fvGeometry The finite volume geometry of the element
124 * \param elemVolVars The element volume variables
125 * \param elemFluxVarsCache Flux variables caches for all faces in stencil
126 * \param scvf The sub-control volume face
127 *
128 * Negative values mean influx.
129 */
130 86030 NeumannFluxes neumann(const Element& element,
131 const FVElementGeometry& fvGeometry,
132 const ElementVolumeVariables& elemVolVars,
133 const ElementFluxVariablesCache& elemFluxVarsCache,
134 const SubControlVolumeFace& scvf) const
135 {
136 86030 NeumannFluxes values(0.0);
137
2/2
✓ Branch 0 taken 14700 times.
✓ Branch 1 taken 57330 times.
 86030 const auto& globalPos = scvf.ipGlobal();
138
4/4
✓ Branch 0 taken 14700 times.
✓ Branch 1 taken 57330 times.
✓ Branch 2 taken 14700 times.
✓ Branch 3 taken 57330 times.
 172060 const auto& volVars = elemVolVars[scvf.insideScvIdx()];
139 86030 Scalar boundaryLayerThickness = 0.0016;
140 //right side
141
10/10
✓ Branch 0 taken 17500 times.
✓ Branch 1 taken 68530 times.
✓ Branch 2 taken 17500 times.
✓ Branch 3 taken 68530 times.
✓ Branch 4 taken 17500 times.
✓ Branch 5 taken 68530 times.
✓ Branch 6 taken 17500 times.
✓ Branch 7 taken 68530 times.
✓ Branch 8 taken 17500 times.
✓ Branch 9 taken 68530 times.
 430150 if (globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_)
142 {
143 17500 Scalar moleFracH2OInside = volVars.moleFraction(FluidSystem::gasPhaseIdx, FluidSystem::H2OIdx);
144 17500 Scalar moleFracRefH2O = 0.0;
145 17500 Scalar evaporationRate = volVars.diffusionCoefficient(FluidSystem::gasPhaseIdx, FluidSystem::AirIdx, FluidSystem::H2OIdx)
146 17500 * (moleFracH2OInside - moleFracRefH2O)
147 17500 / boundaryLayerThickness
148 17500 * volVars.molarDensity(FluidSystem::gasPhaseIdx);
149 17500 values[Indices::conti0EqIdx+2] = evaporationRate;
150
151 17500 values[Indices::energyEqIdx] = FluidSystem::enthalpy(volVars.fluidState(), FluidSystem::gasPhaseIdx) * evaporationRate;
152 52500 values[Indices::energyEqIdx] += FluidSystem::thermalConductivity(volVars.fluidState(), FluidSystem::gasPhaseIdx)
153 35000 * (volVars.temperature() - this->spatialParams().temperatureAtPos(globalPos))
154 17500 /boundaryLayerThickness;
155 }
156 86030 return values;
157 }
158
159
160
161 /*!
162 * \brief Evaluates the initial value for a control volume.
163 *
164 * \param globalPos The center of the finite volume which ought to be set.
165 *
166 * For this method, the \a values parameter stores primary
167 * variables.
168 */
169 ✗ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
170 {
171 862 return initial_(globalPos);
172 }
173
174 /*!
175 * \brief Adds additional VTK output data to the VTKWriter.
176 *
177 * Function is called by the output module on every write.
178 */
179 template<class VTKWriter>
180 2 void addVtkFields(VTKWriter& vtk)
181 {
182
5/12
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 2 times.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 4 vtk.addField(xEquilxwn_, "xEquil^Air_liq");
183
4/10
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 4 vtk.addField(xEquilxnw_, "xEquil^H2O_gas");
184 2 }
185
186 16 void updateVtkFields(const SolutionVector& curSol)
187 {
188 32 const auto& gridView = this->gridGeometry().gridView();
189 16 xEquilxwn_.resize(gridView.size(dofCodim));
190 16 xEquilxnw_.resize(gridView.size(dofCodim));
191
2/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
 32 auto fvGeometry = localView(this->gridGeometry());
192
5/10
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 3216 times.
✗ Branch 10 not taken.
✓ Branch 13 taken 1600 times.
✗ Branch 14 not taken.
 6464 for (const auto& element : elements(this->gridGeometry().gridView()))
193 {
194
2/4
✓ Branch 1 taken 1600 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1600 times.
✗ Branch 5 not taken.
 6400 auto elemSol = elementSolution(element, curSol, this->gridGeometry());
195
1/2
✓ Branch 1 taken 3200 times.
✗ Branch 2 not taken.
 3200 fvGeometry.bindElement(element);
196
197
4/4
✓ Branch 0 taken 8000 times.
✓ Branch 1 taken 3200 times.
✓ Branch 2 taken 6400 times.
✓ Branch 3 taken 1600 times.
 19200 for (auto&& scv : scvs(fvGeometry))
198 {
199 8000 VolumeVariables volVars;
200
1/2
✓ Branch 1 taken 8000 times.
✗ Branch 2 not taken.
 8000 volVars.update(elemSol, *this, element, scv);
201 8000 const auto dofIdxGlobal = scv.dofIndex();
202 16000 xEquilxwn_[dofIdxGlobal] = volVars.xEquil(0,1);
203 24000 xEquilxnw_[dofIdxGlobal] = volVars.xEquil(1,0);
204 }
205 }
206 16 }
207
208 private:
209 // the internal method for the initial condition
210 ✗ PrimaryVariables initial_(const GlobalPosition &globalPos) const
211 {
212 431 PrimaryVariables values(0.0);
213 862 values[Indices::pressureIdx] = 1e5; // water pressure
214 862 values[Indices::switchIdx] = 0.8; // gas saturation
215 862 values[2] = 5e-4; // xwn higher than equilibrium, equilibrium is 3.4e-5
216 862 values[3] = 1e-2; // xnw lower than 1.3e-2
217 1293 values[Indices::temperatureIdx] = this->spatialParams().temperatureAtPos(globalPos);
218 862 values.setState(Indices::bothPhases);
219
220 ✗ return values;
221 }
222
223 static constexpr Scalar eps_ = 1e-6;
224 std::string name_;
225
226 std::shared_ptr<GridVariables> gridVariables_;
227 std::vector<double> xEquilxnw_;
228 std::vector<double> xEquilxwn_;
229 };
230 } // end namespace Dumux
231
232 #endif
233