Reaktoro  v2.11.0
A unified framework for modeling chemically reactive systems
TransportSolver.hpp
1 // // Reaktoro is a unified framework for modeling chemically reactive systems.
2 // //
3 // // Copyright © 2014-2024 Allan Leal
4 // //
5 // // This library is free software; you can redistribute it and/or
6 // // modify it under the terms of the GNU Lesser General Public
7 // // License as published by the Free Software Foundation; either
8 // // version 2.1 of the License, or (at your option) any later version.
9 // //
10 // // This library is distributed in the hope that it will be useful,
11 // // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 // // Lesser General Public License for more details.
14 // //
15 // // You should have received a copy of the GNU Lesser General Public License
16 // // along with this library. If not, see <http://www.gnu.org/licenses/>.
17 
18 // #pragma once
19 
20 // // C++ includes
21 // #include <memory>
22 // #include <vector>
23 
24 // // Reaktoro includes
25 // #include <Reaktoro/Common/Index.hpp>
26 // #include <Reaktoro/Common/Matrix.hpp>
27 // #include <Reaktoro/Common/StringList.hpp>
28 // #include <Reaktoro/Core/ChemicalOutput.hpp>
29 // #include <Reaktoro/Core/ChemicalProps.hpp>
30 // #include <Reaktoro/Core/ChemicalState.hpp>
31 // #include <Reaktoro/Core/ChemicalSystem.hpp>
32 // #include <Reaktoro/Equilibrium/EquilibriumSolver.hpp>
33 
34 // namespace Reaktoro {
35 
36 // // Forward declarations
37 // //class ChemicalProps;
38 // //class ChemicalState;
39 // //class ChemicalSystem;
40 
41 // //class BoundaryState
42 // //{
43 // //public:
44 // // BoundaryState(const ChemicalState& state);
45 // //
46 // //private:
47 // //};
48 // //
49 // //class BoundaryFlux
50 // //{
51 // //public:
52 // // BoundaryFlux(const ChemicalState& state);
53 // //
54 // //private:
55 // //};
56 
57 
58 // class ChemicalField
59 // {
60 // public:
61 // using Iterator = std::vector<ChemicalState>::iterator;
62 
63 // using ConstIterator = std::vector<ChemicalState>::const_iterator;
64 
65 // ChemicalField(Index size, const ChemicalSystem& system);
66 
67 // ChemicalField(Index size, const ChemicalState& state);
68 
69 // auto size() const -> Index { return m_size; }
70 
71 // auto begin() const -> ConstIterator { return m_states.cbegin(); }
72 
73 // auto begin() -> Iterator { return m_states.begin(); }
74 
75 // auto end() const -> ConstIterator { return m_states.cend(); }
76 
77 // auto end() -> Iterator { return m_states.end(); }
78 
79 // auto operator[](Index index) const -> const ChemicalState& { return m_states[index]; }
80 
81 // auto operator[](Index index) -> ChemicalState& { return m_states[index]; }
82 
83 // auto set(const ChemicalState& state) -> void;
84 
85 // auto temperature(VectorXdRef values) -> void;
86 
87 // auto pressure(VectorXdRef values) -> void;
88 
89 // auto elementAmounts(VectorXdRef values) -> void;
90 
91 // auto output(std::string filename, StringList quantities) -> void;
92 
93 // private:
94 // /// The number of degrees of freedom in the chemical field.
95 // Index m_size;
96 
97 // // VectorXd temperatures;
98 // //
99 // // VectorXd pressures;
100 // //
101 // // /// The matrix of amounts for every element (
102 // // Matrix element_amounts;
103 
104 // /// The chemical system common to all degrees of freedom in the chemical field.
105 // ChemicalSystem m_system;
106 
107 // /// The chemical states in the chemical field
108 // std::vector<ChemicalState> m_states;
109 
110 // /// The chemical states in the chemical field
111 // std::vector<ChemicalProps> m_props;
112 // };
113 
114 // /// A class that defines a Tridiagonal Matrix used on TransportSolver.
115 // /// it stores data in a Eigen::VectorXd like, M = {a[0][0], a[0][1], a[0][2],
116 // /// a[1][0], a[1][1], a[1][2],
117 // /// a[2][0], a[2][1], a[2][2]}
118 // class TridiagonalMatrix
119 // {
120 // public:
121 // TridiagonalMatrix() : TridiagonalMatrix(0) {}
122 
123 // TridiagonalMatrix(Index size) : m_size(size), m_data(size * 3) {}
124 
125 // auto size() const -> Index { return m_size; }
126 
127 // auto data() -> VectorXdRef { return m_data; }
128 
129 // auto data() const -> VectorXdConstRef { return m_data; }
130 
131 // auto row(Index index) -> VectorXdRef { return m_data.segment(3 * index, 3); }
132 
133 // auto row(Index index) const -> VectorXdConstRef { return m_data.segment(3 * index, 3); }
134 
135 // auto a() -> VectorXdStridedRef { return VectorXd::Map(m_data.data() + 3, size() - 1, Eigen::InnerStride<3>()); }
136 
137 // auto a() const -> VectorXdConstRef { return VectorXd::Map(m_data.data() + 3, size() - 1, Eigen::InnerStride<3>()); }
138 
139 // auto b() -> VectorXdStridedRef { return VectorXd::Map(m_data.data() + 1, size(), Eigen::InnerStride<3>()); }
140 
141 // auto b() const -> VectorXdConstRef { return VectorXd::Map(m_data.data() + 1, size(), Eigen::InnerStride<3>()); }
142 
143 // auto c() -> VectorXdStridedRef { return VectorXd::Map(m_data.data() + 2, size() - 1, Eigen::InnerStride<3>()); }
144 
145 // auto c() const -> VectorXdConstRef { return VectorXd::Map(m_data.data() + 2, size() - 1, Eigen::InnerStride<3>()); }
146 
147 // auto resize(Index size) -> void;
148 
149 // /// Factorize the tridiagonal matrix to help with linear system A x = b.
150 // auto factorize() -> void;
151 
152 // /// Solve a linear system Ax = b with LU decomposition.
153 // auto solve(VectorXdRef x, VectorXdConstRef b) const -> void;
154 
155 // /// Solve a linear system Ax = b with LU decomposition, using x as the unknown and it's.
156 // /// old values as the vector b.
157 // auto solve(VectorXdRef x) const -> void;
158 
159 // operator MatrixXd() const;
160 
161 // private:
162 // /// The size of the tridiagonal matrix
163 // Index m_size;
164 
165 // /// The coefficients
166 // VectorXd m_data;
167 // };
168 
169 // /// A class that defines the mesh for TransportSolver.
170 // class Mesh
171 // {
172 // public:
173 // Mesh();
174 
175 // Mesh(Index num_cells, double xl = 0.0, double xr = 1.0);
176 
177 // auto setDiscretization(Index num_cells, double xl = 0.0, double xr = 1.0) -> void;
178 
179 // auto numCells() const -> Index { return m_num_cells; }
180 
181 // auto xl() const -> double { return m_xl; }
182 
183 // auto xr() const -> double { return m_xr; }
184 
185 // auto dx() const -> double { return m_dx; }
186 
187 // auto xcells() const -> VectorXdConstRef { return m_xcells; }
188 
189 // private:
190 // /// The number of cells in the discretization.
191 // Index m_num_cells = 10;
192 
193 // /// The x-coordinate of the left boundary (in m).
194 // double m_xl = 0.0;
195 
196 // /// The x-coordinate of the right boundary (in m).
197 // double m_xr = 1.0;
198 
199 // /// The length of the cells (in m).
200 // double m_dx = 0.1;
201 
202 // /// The x-coordinate of the center of the cells.
203 // VectorXd m_xcells;
204 // };
205 
206 // /// A class for solving advection-diffusion problem.
207 // /// Eq: du/dt + v*du/dx = D*d�u/dx�
208 // /// u - amount
209 // /// v - velocity
210 // /// D - diffusion coefficient
211 // class TransportSolver
212 // {
213 // public:
214 // /// Construct a default TransportSolver instance.
215 // TransportSolver();
216 
217 // /// Set the mesh for the numerical solution of the transport problem.
218 // auto setMesh(const Mesh& mesh) -> void { mesh_ = mesh; }
219 
220 // /// Set the velocity for the transport problem.
221 // /// @param val The velocity (in m/s)
222 // auto setVelocity(double val) -> void { velocity = val; }
223 
224 // /// Set the diffusion coefficient for the transport problem.
225 // /// @param val The diffusion coefficient (in m^2/s)
226 // auto setDiffusionCoeff(double val) -> void { diffusion = val; }
227 
228 // /// Set the value of the variable on the boundary.
229 // /// @param val The boundary value for the variable (same unit considered for u).
230 // auto setBoundaryValue(double val) -> void { ul = val; };
231 
232 // /// Set the time step for the numerical solution of the transport problem.
233 // auto setTimeStep(double val) -> void { dt = val; }
234 
235 // /// Return the mesh.
236 // auto mesh() const -> const Mesh& { return mesh_; }
237 
238 // /// Initialize the transport solver before method @ref step is executed.
239 // /// Setup coefficient matrix of the diffusion problem and factorize.
240 // auto initialize() -> void;
241 
242 // /// Step the transport solver.
243 // /// This method solve one step of the transport solver equation, using an explicit approach for
244 // /// advection and total implicit for diffusion. The amount resulted from the advection it is
245 // /// passed to diffusion problem as a "source".
246 // /// @param[in,out] u The solution vector
247 // /// @param q The source rates vector ([same unit considered for u]/m)
248 // auto step(VectorXdRef u, VectorXdConstRef q) -> void;
249 
250 // /// Step the transport solver.
251 // /// @param[in,out] u The solution vector
252 // auto step(VectorXdRef u) -> void;
253 
254 // private:
255 // /// The mesh describing the discretization of the domain.
256 // Mesh mesh_;
257 
258 // /// The time step used to solve the transport problem (in s).
259 // double dt = 0.0;
260 
261 // /// The velocity in the transport problem (in m/s).
262 // double velocity = 0.0;
263 
264 // /// The diffusion coefficient in the transport problem (in m^2/s).
265 // double diffusion = 0.0;
266 
267 // /// The value of the variable on the left boundary.
268 // double ul;
269 
270 // /// The coefficient matrix from the discretized transport equation
271 // TridiagonalMatrix A;
272 
273 // /// The flux limiters at each cell.
274 // VectorXd phi;
275 
276 // /// The previous state of the variables.
277 // VectorXd u0;
278 // };
279 
280 // /// Use this class for solving reactive transport problems.
281 // class ReactiveTransportSolver
282 // {
283 // public:
284 // /// Construct a default ReactiveTransportSolver instance.
285 // ReactiveTransportSolver(const ChemicalSystem& system);
286 
287 // auto setMesh(const Mesh& mesh) -> void;
288 
289 // auto setVelocity(double val) -> void;
290 
291 // auto setDiffusionCoeff(double val) -> void;
292 
293 // auto setBoundaryState(const ChemicalState& state) -> void;
294 
295 // auto setTimeStep(double val) -> void;
296 
297 // auto system() const -> const ChemicalSystem& { return system_; }
298 
299 // auto output() -> ChemicalOutput;
300 
301 // auto initialize(const ChemicalField& field) -> void;
302 
303 // auto step(ChemicalField& field) -> void;
304 
305 // private:
306 // /// The chemical system common to all degrees of freedom in the chemical field.
307 // ChemicalSystem system_;
308 
309 // /// The solver for solving the transport equations
310 // TransportSolver transportsolver;
311 
312 // /// The solver for solving the equilibrium equations
313 // EquilibriumSolver equilibriumsolver;
314 
315 // /// The list of chemical output objects
316 // std::vector<ChemicalOutput> outputs;
317 
318 // /// The amounts of fluid elements on the boundary.
319 // VectorXd bbc;
320 
321 // /// The amounts of a fluid element on each cell of the mesh.
322 // MatrixXd bf;
323 
324 // /// The amounts of a solid element on each cell of the mesh.
325 // MatrixXd bs;
326 
327 // /// The amounts of an element on each cell of the mesh.
328 // MatrixXd b;
329 
330 // /// The current number of steps in the solution of the reactive transport equations.
331 // Index steps = 0;
332 // };
333 
334 // } // namespace Reaktoro