GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/2pnc/fuelcell/problem.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 76 81 93.8%
Functions: 15 24 62.5%
Branches: 125 223 56.1%
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 TwoPNCTests
10 * \brief Definition of a problem for water management in PEM fuel cells.
11 */
12
13 #ifndef DUMUX_FUELCELL_PROBLEM_HH
14 #define DUMUX_FUELCELL_PROBLEM_HH
15
16 #include <dumux/common/properties.hh>
17 #include <dumux/common/parameters.hh>
18 #include <dumux/common/boundarytypes.hh>
19 #include <dumux/common/numeqvector.hh>
20
21 #include <dumux/porousmediumflow/problem.hh>
22
23 #ifdef NONISOTHERMAL
24 #include <dumux/material/chemistry/electrochemistry/electrochemistryni.hh>
25 #else
26 #include <dumux/material/chemistry/electrochemistry/electrochemistry.hh>
27 #endif
28
29 namespace Dumux {
30
31 /*!
32 * \ingroup TwoPNCTests
33 * \brief Problem for water management in PEM fuel cells.
34 *
35 */
36 template <class TypeTag>
37 class FuelCellProblem : public PorousMediumFlowProblem<TypeTag>
38 {
39 using ParentType = PorousMediumFlowProblem<TypeTag>;
40
41 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
42 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
43 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
44 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
45 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
46 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
47 using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
48 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
49 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
50 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
51 using Element = typename GridView::template Codim<0>::Entity;
52 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
53 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
54 using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
55 // Select the electrochemistry method
56 #ifdef NONISOTHERMAL
57 using ElectroChemistry = typename Dumux::ElectroChemistryNI<Scalar, Indices, FluidSystem, GridGeometry, ElectroChemistryModel::Ochs>;
58 #else
59 using ElectroChemistry = typename Dumux::ElectroChemistry<Scalar, Indices, FluidSystem, GridGeometry, ElectroChemistryModel::Ochs>;
60 #endif
61 static constexpr int dim = GridView::dimension;
62 static constexpr int dimWorld = GridView::dimensionworld;
63 static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethods::box;
64 using GlobalPosition = typename SubControlVolume::GlobalPosition;
65
66 enum { dofCodim = isBox ? dim : 0 };
67 public:
68 3 FuelCellProblem(std::shared_ptr<const GridGeometry> gridGeometry)
69
12/40
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 3 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 3 times.
✓ Branch 15 taken 3 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 3 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 3 times.
✗ Branch 22 not taken.
✓ Branch 24 taken 3 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 3 times.
✗ Branch 28 not taken.
✓ Branch 30 taken 3 times.
✗ Branch 31 not taken.
✓ Branch 33 taken 3 times.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
 9 : ParentType(gridGeometry)
70 {
71
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 nTemperature_ = getParam<int>("Problem.NTemperature");
72
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 nPressure_ = getParam<int>("Problem.NPressure");
73
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 pressureLow_ = getParam<Scalar>("Problem.PressureLow");
74
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 pressureHigh_ = getParam<Scalar>("Problem.PressureHigh");
75
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 temperatureLow_ = getParam<Scalar>("Problem.TemperatureLow");
76
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 temperatureHigh_ = getParam<Scalar>("Problem.TemperatureHigh");
77
78
2/4
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 3 times.
 3 name_ = getParam<std::string>("Problem.Name");
79
80
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 pO2Inlet_ = getParam<Scalar>("ElectroChemistry.pO2Inlet");
81
82 6 FluidSystem::init(/*Tmin=*/temperatureLow_,
83 /*Tmax=*/temperatureHigh_,
84 3 /*nT=*/nTemperature_,
85 /*pmin=*/pressureLow_,
86 /*pmax=*/pressureHigh_,
87
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 /*np=*/nPressure_);
88 3 }
89
90
91 /*!
92 * \brief The problem name.
93 *
94 * This is used as a prefix for files generated by the simulation.
95 */
96 const std::string& name() const
97
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 { return name_; }
98
99 //! \copydoc Dumux::FVProblem::source()
100 501480 NumEqVector source(const Element &element,
101 const FVElementGeometry& fvGeometry,
102 const ElementVolumeVariables& elemVolVars,
103 const SubControlVolume &scv) const
104 {
105 501480 NumEqVector values(0.0);
106
2/2
✓ Branch 0 taken 43092 times.
✓ Branch 1 taken 458388 times.
 501480 const auto& globalPos = scv.dofPosition();
107
108 //reaction sources from electro chemistry
109 1002960 if(inReactionLayer_(globalPos))
110 {
111 43092 const auto& volVars = elemVolVars[scv];
112 43092 auto currentDensity = ElectroChemistry::calculateCurrentDensity(volVars);
113
4/10
✓ Branch 1 taken 43092 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43092 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 43092 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 43092 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 86184 ElectroChemistry::reactionSource(values, currentDensity);
114 }
115
116 501480 return values;
117 }
118
119
120
121 /*!
122 * \brief Specifies which kind of boundary condition should be
123 * used for which equation on a given boundary segment
124 *
125 * \param globalPos The global position
126 */
127 22882 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
128 {
129 22882 BoundaryTypes bcTypes;
130 22882 bcTypes.setAllNeumann();
131
132 45764 if (onUpperBoundary_(globalPos)){
133 bcTypes.setAllDirichlet();
134 }
135 22882 return bcTypes;
136 }
137
138 /*!
139 * \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
140 *
141 * \param globalPos The global position
142 */
143 ✗ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
144 {
145 9198 auto priVars = initial_(globalPos);
146
147 9198 if(onUpperBoundary_(globalPos))
148 {
149 Scalar pn = 1.0e5;
150 priVars[Indices::pressureIdx] = pn;
151 priVars[Indices::switchIdx] = 0.3;//Sw for bothPhases
152 priVars[Indices::switchIdx+1] = pO2Inlet_/4.315e9; //moleFraction xlO2 for bothPhases
153 #ifdef NONISOTHERMAL
154 priVars[Indices::temperatureIdx] = this->spatialParams().temperatureAtPos(globalPos);
155 #endif
156 }
157
158 ✗ return priVars;
159 }
160
161
162
163 /*!
164 * \brief Evaluates the initial values for a control volume.
165 *
166 * \param globalPos The global position
167 */
168 ✗ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
169 868 { return initial_(globalPos); }
170
171
172 /*!
173 * \brief Adds additional VTK output data to the VTKWriter.
174 *
175 * Function is called by the output module on every write.
176 */
177 template<class VTKWriter>
178 3 void addVtkFields(VTKWriter& vtk)
179 {
180 6 const auto& gridView = this->gridGeometry().gridView();
181 3 currentDensity_.resize(gridView.size(dofCodim));
182 3 reactionSourceH2O_.resize(gridView.size(dofCodim));
183 3 reactionSourceO2_.resize(gridView.size(dofCodim));
184 3 Kxx_.resize(gridView.size(dofCodim));
185 3 Kyy_.resize(gridView.size(dofCodim));
186
187
5/12
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 3 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 9 vtk.addField(currentDensity_, "currentDensity [A/cm^2]");
188
5/12
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 3 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 9 vtk.addField(reactionSourceH2O_, "reactionSourceH2O [mol/(sm^2)]");
189
5/12
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 3 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 9 vtk.addField(reactionSourceO2_, "reactionSourceO2 [mol/(sm^2)]");
190
5/12
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 3 times.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 6 vtk.addField(Kxx_, "Kxx");
191
4/10
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 6 vtk.addField(Kyy_, "Kyy");
192 3 }
193
194 45 void updateVtkFields(const SolutionVector& curSol)
195 {
196
2/4
✓ Branch 1 taken 45 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 45 times.
✗ Branch 5 not taken.
 90 auto fvGeometry = localView(this->gridGeometry());
197
5/10
✓ Branch 1 taken 45 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 45 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 45 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 5715 times.
✗ Branch 10 not taken.
✓ Branch 13 taken 3780 times.
✗ Branch 14 not taken.
 11520 for (const auto& element : elements(this->gridGeometry().gridView()))
198 {
199
2/4
✓ Branch 1 taken 3780 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3780 times.
✗ Branch 5 not taken.
 11340 auto elemSol = elementSolution(element, curSol, this->gridGeometry());
200
1/2
✓ Branch 1 taken 5670 times.
✗ Branch 2 not taken.
 5670 fvGeometry.bindElement(element);
201
202
4/4
✓ Branch 0 taken 17010 times.
✓ Branch 1 taken 5670 times.
✓ Branch 2 taken 15120 times.
✓ Branch 3 taken 3780 times.
 41580 for (auto&& scv : scvs(fvGeometry))
203 {
204 17010 VolumeVariables volVars;
205
1/2
✓ Branch 1 taken 17010 times.
✗ Branch 2 not taken.
 17010 volVars.update(elemSol, *this, element, scv);
206
2/2
✓ Branch 0 taken 1575 times.
✓ Branch 1 taken 15435 times.
 17010 const auto& globalPos = scv.dofPosition();
207
2/2
✓ Branch 0 taken 315 times.
✓ Branch 1 taken 1575 times.
 17010 const auto dofIdxGlobal = scv.dofIndex();
208
209 34020 if(inReactionLayer_(globalPos))
210 {
211 //reactionSource Output
212
1/2
✓ Branch 1 taken 1575 times.
✗ Branch 2 not taken.
 1575 PrimaryVariables source;
213
1/2
✓ Branch 1 taken 1575 times.
✗ Branch 2 not taken.
 1575 auto i = ElectroChemistry::calculateCurrentDensity(volVars);
214
4/10
✓ Branch 1 taken 1575 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1575 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1575 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 1575 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 3150 ElectroChemistry::reactionSource(source, i);
215
216 3150 reactionSourceH2O_[dofIdxGlobal] = source[Indices::conti0EqIdx + FluidSystem::H2OIdx];
217 3150 reactionSourceO2_[dofIdxGlobal] = source[Indices::conti0EqIdx + FluidSystem::O2Idx];
218
219 //Current Output in A/cm^2
220 3150 currentDensity_[dofIdxGlobal] = i/10000;
221 }
222 else
223 {
224 15435 reactionSourceH2O_[dofIdxGlobal] = 0.0;
225 15435 reactionSourceO2_[dofIdxGlobal] = 0.0;
226 30870 currentDensity_[dofIdxGlobal] = 0.0;
227 }
228 51030 Kxx_[dofIdxGlobal] = volVars.permeability()[0][0];
229 68040 Kyy_[dofIdxGlobal] = volVars.permeability()[1][1];
230 }
231 }
232 45 }
233
234 private:
235
236 ✗ PrimaryVariables initial_(const GlobalPosition &globalPos) const
237 {
238 5033 PrimaryVariables priVars(0.0);
239 4599 priVars.setState(Indices::bothPhases);
240
241 5033 Scalar pn = 1.0e5;
242 4599 priVars[Indices::pressureIdx] = pn;
243 4599 priVars[Indices::switchIdx] = 0.3;//Sw for bothPhases
244 4599 priVars[Indices::switchIdx+1] = pO2Inlet_/4.315e9; //moleFraction xlO2 for bothPhases
245 #ifdef NONISOTHERMAL
246 4032 priVars[Indices::temperatureIdx] = this->spatialParams().temperatureAtPos(globalPos);
247 #endif
248
249 ✗ return priVars;
250 }
251
252 bool onUpperBoundary_(const GlobalPosition &globalPos) const
253
15/15
✓ Branch 0 taken 3255 times.
✓ Branch 1 taken 5115 times.
✓ Branch 2 taken 3255 times.
✓ Branch 3 taken 5115 times.
✓ Branch 4 taken 3255 times.
✓ Branch 5 taken 11751 times.
✓ Branch 6 taken 11131 times.
✓ Branch 7 taken 11751 times.
✓ Branch 8 taken 11131 times.
✓ Branch 9 taken 11751 times.
✓ Branch 10 taken 7876 times.
✓ Branch 11 taken 6636 times.
✓ Branch 12 taken 7876 times.
✓ Branch 13 taken 6636 times.
✓ Branch 14 taken 7876 times.
 137405 { return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; }
254
255 bool inReactionLayer_(const GlobalPosition& globalPos) const
256
32/32
✓ Branch 0 taken 43092 times.
✓ Branch 1 taken 458388 times.
✓ Branch 2 taken 43092 times.
✓ Branch 3 taken 458388 times.
✓ Branch 4 taken 43092 times.
✓ Branch 5 taken 458388 times.
✓ Branch 6 taken 43092 times.
✓ Branch 7 taken 458388 times.
✓ Branch 8 taken 43092 times.
✓ Branch 9 taken 458388 times.
✓ Branch 10 taken 43092 times.
✓ Branch 11 taken 458388 times.
✓ Branch 12 taken 43092 times.
✓ Branch 13 taken 458388 times.
✓ Branch 14 taken 43092 times.
✓ Branch 15 taken 458388 times.
✓ Branch 16 taken 1575 times.
✓ Branch 17 taken 15435 times.
✓ Branch 18 taken 1575 times.
✓ Branch 19 taken 15435 times.
✓ Branch 20 taken 1575 times.
✓ Branch 21 taken 15435 times.
✓ Branch 22 taken 1575 times.
✓ Branch 23 taken 15435 times.
✓ Branch 24 taken 1575 times.
✓ Branch 25 taken 15435 times.
✓ Branch 26 taken 1575 times.
✓ Branch 27 taken 15435 times.
✓ Branch 28 taken 1575 times.
✓ Branch 29 taken 15435 times.
✓ Branch 30 taken 1575 times.
✓ Branch 31 taken 15435 times.
 4147920 { return globalPos[1] < 0.1*(this->gridGeometry().bBoxMax()[1] - this->gridGeometry().bBoxMin()[1]) + eps_; }
257
258 static constexpr Scalar eps_ = 1e-6;
259 int nTemperature_;
260 int nPressure_;
261 std::string name_ ;
262 Scalar pressureLow_, pressureHigh_;
263 Scalar temperatureLow_, temperatureHigh_;
264 Scalar pO2Inlet_;
265 std::vector<double> currentDensity_;
266 std::vector<double> reactionSourceH2O_;
267 std::vector<double> reactionSourceO2_;
268 std::vector<double> Kxx_;
269 std::vector<double> Kyy_;
270 };
271
272 } // end namespace Dumux
273
274 #endif
275