GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: dumux/test/porousmediumflow/2pnc/fuelcell/problem.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 84 84 100.0%
Functions: 15 15 100.0%
Branches: 66 104 63.5%

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-FileCopyrightText: 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
3/6
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 3 times.
✗ Branch 9 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
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
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 3 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 3 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
2/2
✓ Branch 0 taken 43092 times.
✓ Branch 1 taken 458388 times.
501480 NumEqVector values(0.0);
106 501480 const auto& globalPos = scv.dofPosition();
107
108 //reaction sources from electro chemistry
109
2/2
✓ Branch 0 taken 43092 times.
✓ Branch 1 taken 458388 times.
501480 if(inReactionLayer_(globalPos))
110 {
111 43092 const auto& volVars = elemVolVars[scv];
112 43092 auto currentDensity = ElectroChemistry::calculateCurrentDensity(volVars);
113
1/2
✓ Branch 2 taken 43092 times.
✗ Branch 3 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
2/2
✓ Branch 0 taken 76080 times.
✓ Branch 1 taken 22882 times.
98962 BoundaryTypes bcTypes;
130 22882 bcTypes.setAllNeumann();
131
132
2/2
✓ Branch 0 taken 9891 times.
✓ Branch 1 taken 12991 times.
22882 if (onUpperBoundary_(globalPos)){
133 22882 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 4599 PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
144 {
145
1/2
✓ Branch 1 taken 630 times.
✗ Branch 2 not taken.
4599 auto priVars = initial_(globalPos);
146
147 if(onUpperBoundary_(globalPos))
148 {
149
1/2
✓ Branch 1 taken 1281 times.
✗ Branch 2 not taken.
4599 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
1/2
✓ Branch 1 taken 630 times.
✗ Branch 2 not taken.
2016 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 434 PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
169 434 { 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 3 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
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
3 vtk.addField(currentDensity_, "currentDensity [A/cm^2]");
188
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
3 vtk.addField(reactionSourceH2O_, "reactionSourceH2O [mol/(sm^2)]");
189
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
3 vtk.addField(reactionSourceO2_, "reactionSourceO2 [mol/(sm^2)]");
190
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
3 vtk.addField(Kxx_, "Kxx");
191
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
3 vtk.addField(Kyy_, "Kyy");
192 3 }
193
194 45 void updateVtkFields(const SolutionVector& curSol)
195 {
196
1/2
✓ Branch 1 taken 45 times.
✗ Branch 2 not taken.
45 auto fvGeometry = localView(this->gridGeometry());
197
3/6
✓ Branch 1 taken 45 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5715 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3780 times.
✗ Branch 8 not taken.
17055 for (const auto& element : elements(this->gridGeometry().gridView()))
198 {
199
2/3
✓ Branch 1 taken 3780 times.
✓ Branch 2 taken 1890 times.
✗ Branch 3 not taken.
5670 auto elemSol = elementSolution(element, curSol, this->gridGeometry());
200
1/2
✓ Branch 1 taken 3780 times.
✗ Branch 2 not taken.
7560 fvGeometry.bindElement(element);
201
202
2/2
✓ Branch 0 taken 17010 times.
✓ Branch 1 taken 5670 times.
22680 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 17010 const auto dofIdxGlobal = scv.dofIndex();
208
209
2/2
✓ Branch 0 taken 1575 times.
✓ Branch 1 taken 15435 times.
17010 if(inReactionLayer_(globalPos))
210 {
211 //reactionSource Output
212 1575 PrimaryVariables source;
213
2/4
✓ Branch 1 taken 1575 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1575 times.
✗ Branch 5 not taken.
1575 auto i = ElectroChemistry::calculateCurrentDensity(volVars);
214
2/4
✓ Branch 1 taken 1575 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1575 times.
✗ Branch 5 not taken.
3150 ElectroChemistry::reactionSource(source, i);
215
216 1575 reactionSourceH2O_[dofIdxGlobal] = source[Indices::conti0EqIdx + FluidSystem::H2OIdx];
217 1575 reactionSourceO2_[dofIdxGlobal] = source[Indices::conti0EqIdx + FluidSystem::O2Idx];
218
219 //Current Output in A/cm^2
220 1575 currentDensity_[dofIdxGlobal] = i/10000;
221 }
222 else
223 {
224 15435 reactionSourceH2O_[dofIdxGlobal] = 0.0;
225 15435 reactionSourceO2_[dofIdxGlobal] = 0.0;
226 15435 currentDensity_[dofIdxGlobal] = 0.0;
227 }
228 17010 Kxx_[dofIdxGlobal] = volVars.permeability()[0][0];
229 17010 Kyy_[dofIdxGlobal] = volVars.permeability()[1][1];
230 }
231 }
232 45 }
233
234 private:
235
236 5033 PrimaryVariables initial_(const GlobalPosition &globalPos) const
237 {
238
1/2
✓ Branch 1 taken 1911 times.
✗ Branch 2 not taken.
9940 PrimaryVariables priVars(0.0);
239 5033 priVars.setState(Indices::bothPhases);
240
241
1/2
✓ Branch 1 taken 1911 times.
✗ Branch 2 not taken.
4907 Scalar pn = 1.0e5;
242 5033 priVars[Indices::pressureIdx] = pn;
243 5033 priVars[Indices::switchIdx] = 0.3;//Sw for bothPhases
244
1/2
✓ Branch 1 taken 1281 times.
✗ Branch 2 not taken.
4907 priVars[Indices::switchIdx+1] = pO2Inlet_/4.315e9; //moleFraction xlO2 for bothPhases
245 #ifdef NONISOTHERMAL
246
1/2
✓ Branch 1 taken 630 times.
✗ Branch 2 not taken.
2170 priVars[Indices::temperatureIdx] = this->spatialParams().temperatureAtPos(globalPos);
247 #endif
248
249 return priVars;
250 }
251
252
2/2
✓ Branch 0 taken 9891 times.
✓ Branch 1 taken 12991 times.
22882 bool onUpperBoundary_(const GlobalPosition &globalPos) const
253
4/5
✓ Branch 1 taken 6375 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 7287 times.
✓ Branch 4 taken 7876 times.
✓ Branch 0 taken 3255 times.
27481 { return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; }
254
255
4/4
✓ Branch 0 taken 43092 times.
✓ Branch 1 taken 458388 times.
✓ Branch 2 taken 1575 times.
✓ Branch 3 taken 15435 times.
518490 bool inReactionLayer_(const GlobalPosition& globalPos) const
256
4/4
✓ Branch 0 taken 43092 times.
✓ Branch 1 taken 458388 times.
✓ Branch 2 taken 1575 times.
✓ Branch 3 taken 15435 times.
518490 { 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