Reaktoro 
A unified framework for modeling chemically reactive systems
ThermoVector.hpp
1 // Reaktoro is a unified framework for modeling chemically reactive systems.
2 //
3 // Copyright (C) 2014-2015 Allan Leal
4 //
5 // This program is free software: you can redistribute it and/or modify
6 // it under the terms of the GNU General Public License as published by
7 // the Free Software Foundation, either version 3 of the License, or
8 // (at your option) any later version.
9 //
10 // This program 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
13 // GNU General Public License for more details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with this program. If not, see <http://www.gnu.org/licenses/>.
17 
18 #pragma once
19 
20 // Reaktoro includes
21 #include <Reaktoro/Math/Matrix.hpp>
22 #include <Reaktoro/Common/ThermoScalar.hpp>
23 
24 namespace Reaktoro {
25 
26 // Forward declaration
27 template<typename V, typename T, typename P>
28 class ThermoVectorBase;
29 
36 using ThermoVector = ThermoVectorBase<Vector,Vector,Vector>;
37 
40 template<typename V, typename T, typename P>
41 class ThermoVectorBase
42 {
43 public:
45  V val;
46 
48  T ddT;
49 
51  P ddP;
52 
55  {}
56 
59  explicit ThermoVectorBase(Index nrows)
60  : ThermoVectorBase(zeros(nrows), zeros(nrows), zeros(nrows)) {}
61 
65  ThermoVectorBase(Index nrows, double val)
66  : ThermoVectorBase(constants(nrows, val), zeros(nrows), zeros(nrows)) {}
67 
72  ThermoVectorBase(const V& val, const T& ddT, const P& ddP)
73  : val(val), ddT(ddT), ddP(ddP) {}
74 
76  template<typename VR, typename TR, typename PR>
78  : val(other.val), ddT(other.ddT), ddP(other.ddP)
79  {}
80 
82  auto size() const -> Index
83  {
84  return val.size();
85  }
86 
89  auto resize(Index nrows) -> void
90  {
91  val.resize(nrows);
92  ddT.resize(nrows);
93  ddP.resize(nrows);
94  }
95 
97  template<typename VR>
98  auto fill(const ThermoScalarBase<VR>& other) -> void
99  {
100  val.fill(other.val);
101  ddT.fill(other.ddT);
102  ddP.fill(other.ddP);
103  }
104 
106  auto fill(double value) -> void
107  {
108  val.fill(value);
109  ddT.fill(0.0);
110  ddP.fill(0.0);
111  }
112 
114  template<typename VR, typename TR, typename PR>
116  {
117  val = other.val;
118  ddT = other.ddT;
119  ddP = other.ddP;
120  return *this;
121  }
122 
124  template<typename VR>
126  {
127  val.fill(other.val);
128  ddT.fill(other.ddT);
129  ddP.fill(other.ddP);
130  return *this;
131  }
132 
134  auto operator=(double other) -> ThermoVectorBase&
135  {
136  val.fill(other);
137  ddT.fill(0.0);
138  ddP.fill(0.0);
139  return *this;
140  }
141 
143  template<typename VR, typename TR, typename PR>
145  {
146  val += other.val;
147  ddT += other.ddT;
148  ddP += other.ddP;
149  return *this;
150  }
151 
153  template<typename VR>
155  {
156  val.array() += other.val;
157  ddT.array() += other.ddT;
158  ddP.array() += other.ddP;
159  return *this;
160  }
161 
163  auto operator+=(double scalar) -> ThermoVectorBase&
164  {
165  val.array() += scalar;
166  return *this;
167  }
168 
170  template<typename VR, typename TR, typename PR>
172  {
173  val -= other.val;
174  ddT -= other.ddT;
175  ddP -= other.ddP;
176  return *this;
177  }
178 
180  template<typename VR>
182  {
183  val.array() -= other.val;
184  ddT.array() -= other.ddT;
185  ddP.array() -= other.ddP;
186  return *this;
187  }
188 
190  auto operator-=(double other) -> ThermoVectorBase&
191  {
192  val.array() -= other;
193  return *this;
194  }
195 
197  template<typename VR, typename TR, typename PR>
199  {
200  ddT = diag(ddT) * other.val + diag(val) * other.ddT;
201  ddP = diag(ddP) * other.val + diag(val) * other.ddP;
202  val = diag(val) * other.val;
203  return *this;
204  }
205 
207  auto operator*=(double scalar) -> ThermoVectorBase&
208  {
209  val *= scalar;
210  ddT *= scalar;
211  ddP *= scalar;
212  return *this;
213  }
214 
216  auto operator/=(double scalar) -> ThermoVectorBase&
217  {
218  *this *= 1.0/scalar;
219  return *this;
220  }
221 
224  {
225  return {val[irow], ddT[irow], ddP[irow]};
226  }
227 
230  {
231  return {val[irow], ddT[irow], ddP[irow]};
232  }
233 
235  explicit operator Vector() const
236  {
237  return val;
238  }
239 };
240 
241 template<typename V, typename T, typename P>
242 auto operator<<(std::ostream& out, const ThermoVectorBase<V,T,P>& a) -> std::ostream&
243 {
244  out << a.val;
245  return out;
246 }
247 
248 template<typename V, typename T, typename P>
249 auto operator+(const ThermoVectorBase<V,T,P>& l) -> ThermoVectorBase<V,T,P>
250 {
251  return l;
252 }
253 
254 template<typename V, typename T, typename P>
255 auto operator-(const ThermoVectorBase<V,T,P>& l) -> ThermoVectorBase<decltype(-l.val), decltype(-l.ddT), decltype(-l.ddP)>
256 {
257  return {-l.val, -l.ddT, -l.ddP};
258 }
259 
260 template<typename VL, typename TL, typename PL, typename VR, typename TR, typename PR>
261 auto operator+(const ThermoVectorBase<VL,TL,PL>& l, const ThermoVectorBase<VR,TR,PR>& r) -> ThermoVectorBase<decltype(l.val + r.val), decltype(l.ddT + r.ddT), decltype(l.ddP + r.ddP)>
262 {
263  return {l.val + r.val, l.ddT + r.ddT, l.ddP + r.ddP};
264 }
265 
266 template<typename V, typename T, typename P>
267 auto operator+(const ThermoVectorBase<V,T,P>& l, const Vector& r) -> ThermoVectorBase<decltype(l.val + r),T,P>
268 {
269  return {l.val + r, l.ddT, l.ddP};
270 }
271 
272 template<typename V, typename T, typename P>
273 auto operator+(const Vector& l, const ThermoVectorBase<V,T,P>& r) -> decltype(r + l)
274 {
275  return r + l;
276 }
277 
278 template<typename VL, typename TL, typename PL, typename VR, typename TR, typename PR>
279 auto operator-(const ThermoVectorBase<VL,TL,PL>& l, const ThermoVectorBase<VR,TR,PR>& r) -> ThermoVectorBase<decltype(l.val - r.val), decltype(l.ddT - r.ddT), decltype(l.ddP - r.ddP)>
280 {
281  return {l.val - r.val, l.ddT - r.ddT, l.ddP - r.ddP};
282 }
283 
284 template<typename V, typename T, typename P>
285 auto operator-(const ThermoVectorBase<V,T,P>& l, const Vector& r) -> ThermoVectorBase<decltype(l.val - r),T,P>
286 {
287  return {l.val - r, l.ddT, l.ddP};
288 }
289 
290 template<typename V, typename T, typename P>
291 auto operator-(const Vector& l, const ThermoVectorBase<V,T,P>& r) -> decltype(-(r - l))
292 {
293  return -(r - l);
294 }
295 
296 template<typename V, typename T, typename P>
297 auto operator*(double l, const ThermoVectorBase<V,T,P>& r) -> ThermoVectorBase<decltype(l * r.val), decltype(l * r.ddT), decltype(l * r.ddP)>
298 {
299  return {l * r.val, l * r.ddT, l * r.ddP};
300 }
301 
302 template<typename V, typename T, typename P>
303 auto operator*(const ThermoVectorBase<V,T,P>& l, double r) -> decltype(r * l)
304 {
305  return r * l;
306 }
307 
308 template<typename VL, typename V, typename T, typename P>
309 auto operator*(const ThermoScalarBase<VL>& l, const ThermoVectorBase<V,T,P>& r) -> ThermoVectorBase<decltype(l.val * r.val), decltype(l.val * r.ddT + l.ddT * r.val), decltype(l.val * r.ddP + l.ddP * r.val)>
310 {
311  return {l.val * r.val, l.val * r.ddT + l.ddT * r.val, l.val * r.ddP + l.ddP * r.val};
312 }
313 
314 template<typename V, typename T, typename P, typename VR>
315 auto operator*(const ThermoVectorBase<V,T,P>& l, const ThermoScalarBase<VR>& r) -> decltype(r * l)
316 {
317  return r * l;
318 }
319 
320 template<typename VL, typename TL, typename PL, typename VR, typename TR, typename PR>
321 auto operator%(const ThermoVectorBase<VL,TL,PL>& l, const ThermoVectorBase<VR,TR,PR>& r) -> ThermoVectorBase<decltype(diag(l.val) * r.val), decltype(diag(l.val) * r.ddT + diag(r.val) * l.ddT), decltype(diag(l.val) * r.ddP + diag(r.val) * l.ddP)>
322 {
323  return {diag(l.val) * r.val,
324  diag(l.val) * r.ddT + diag(r.val) * l.ddT,
325  diag(l.val) * r.ddP + diag(r.val) * l.ddP};
326 }
327 
328 template<typename V, typename T, typename P>
329 auto operator%(const Vector& l, const ThermoVectorBase<V,T,P>& r) -> ThermoVectorBase<decltype(diag(l) * r.val), decltype(diag(l) * r.ddT), decltype(diag(l) * r.ddP)>
330 {
331  return {diag(l) * r.val,
332  diag(l) * r.ddT,
333  diag(l) * r.ddP};
334 }
335 
336 template<typename VL, typename TL, typename PL, typename VR, typename TR, typename PR>
337 auto operator%(const ThermoVectorBase<VL,TL,PL>& l, const Vector& r) -> decltype(r % l)
338 {
339  return r % l;
340 }
341 
342 template<typename VL, typename TL, typename PL, typename VR, typename TR, typename PR>
343 auto operator/(const ThermoVectorBase<VL,TL,PL>& l, const ThermoVectorBase<VR,TR,PR>& r) -> ThermoVector
344 {
345  const Vector tmp = 1.0/(r.val % r.val);
346  return {l.val/r.val,
347  diag(tmp) * (diag(r.val) * l.ddT - diag(l.val) * r.ddT),
348  diag(tmp) * (diag(r.val) * l.ddP - diag(l.val) * r.ddP)};
349 }
350 
351 template<typename V, typename T, typename P, typename VR>
352 auto operator/(const ThermoVectorBase<V,T,P>& l, const ThermoScalarBase<VR>& r) -> ThermoVectorBase<decltype(l.val/r.val), decltype(double() * (l.ddT * r.val - l.val * r.ddT)), decltype(double() * (l.ddP * r.val - l.val * r.ddP))>
353 {
354  const double tmp = 1.0/(r.val * r.val);
355  return {l.val/r.val, tmp * (l.ddT * r.val - l.val * r.ddT), tmp * (l.ddP * r.val - l.val * r.ddP)};
356 }
357 
358 template<typename V, typename T, typename P>
359 auto operator/(const ThermoVectorBase<V,T,P>& l, double r) -> decltype((1.0/r) * l)
360 {
361  return (1.0/r) * l;
362 }
363 
364 template<typename VL, typename TL, typename PL, typename VR, typename TR, typename PR>
365 auto operator<(const ThermoVectorBase<VL,TL,PL>& l, const ThermoVectorBase<VR,TR,PR>& r) -> bool
366 {
367  return l.val < r.val;
368 }
369 
370 template<typename VL, typename TL, typename PL, typename VR, typename TR, typename PR>
371 auto operator<=(const ThermoVectorBase<VL,TL,PL>& l, const ThermoVectorBase<VR,TR,PR>& r) -> bool
372 {
373  return l.val <= r.val;
374 }
375 
376 template<typename VL, typename TL, typename PL, typename VR, typename TR, typename PR>
377 auto operator>(const ThermoVectorBase<VL,TL,PL>& l, const ThermoVectorBase<VR,TR,PR>& r) -> bool
378 {
379  return l.val > r.val;
380 }
381 
382 template<typename VL, typename TL, typename PL, typename VR, typename TR, typename PR>
383 auto operator>=(const ThermoVectorBase<VL,TL,PL>& l, const ThermoVectorBase<VR,TR,PR>& r) -> bool
384 {
385  return l.val >= r.val;
386 }
387 
388 template<typename VL, typename TL, typename PL, typename VR, typename TR, typename PR>
389 auto operator==(const ThermoVectorBase<VL,TL,PL>& l, const ThermoVectorBase<VR,TR,PR>& r) -> bool
390 {
391  return l.val == r.val;
392 }
393 
394 template<typename VL, typename TL, typename PL, typename VR, typename TR, typename PR>
395 auto operator!=(const ThermoVectorBase<VL,TL,PL>& l, const ThermoVectorBase<VR,TR,PR>& r) -> bool
396 {
397  return l.val == r.val;
398 }
399 
401 template<typename V, typename T, typename P>
403 {
404  return {vec.val[irow], vec.ddT[irow], vec.ddP[irow]};
405 }
406 
408 template<typename V, typename T, typename P>
410 {
411  return {vec.val[irow], vec.ddT[irow], vec.ddP[irow]};
412 }
413 
415 template<typename V, typename T, typename P>
416 auto rows(ThermoVectorBase<V,T,P>& vec, Index irow, Index nrows) -> ThermoVectorBase<decltype(vec.val.segment(irow, nrows)), decltype(vec.ddT.segment(irow, nrows)), decltype(vec.ddP.segment(irow, nrows))>
417 {
418  return {vec.val.segment(irow, nrows), vec.ddT.segment(irow, nrows), vec.ddP.segment(irow, nrows)};
419 }
420 
422 template<typename V, typename T, typename P>
423 auto rows(const ThermoVectorBase<V,T,P>& vec, Index irow, Index nrows) -> ThermoVectorBase<decltype(vec.val.segment(irow, nrows)), decltype(vec.ddT.segment(irow, nrows)), decltype(vec.ddP.segment(irow, nrows))>
424 {
425  return {vec.val.segment(irow, nrows), vec.ddT.segment(irow, nrows), vec.ddP.segment(irow, nrows)};
426 }
427 
429 template<typename V, typename T, typename P>
430 auto rows(ThermoVectorBase<V,T,P>& vec, const Indices& irows) -> ThermoVectorBase<decltype(Reaktoro::rows(vec.val, irows)), decltype(Reaktoro::rows(vec.ddT, irows)), decltype(Reaktoro::rows(vec.ddP, irows))>
431 {
432  return {Reaktoro::rows(vec.val, irows), Reaktoro::rows(vec.ddT, irows), Reaktoro::rows(vec.ddP, irows)};
433 }
434 
436 template<typename V, typename T, typename P>
437 auto rows(const ThermoVectorBase<V,T,P>& vec, const Indices& irows) -> ThermoVectorBase<decltype(Reaktoro::rows(vec.val, irows)), decltype(Reaktoro::rows(vec.ddT, irows)), decltype(Reaktoro::rows(vec.ddP, irows))>
438 {
439  return {Reaktoro::rows(vec.val, irows), Reaktoro::rows(vec.ddT, irows), Reaktoro::rows(vec.ddP, irows)};
440 }
441 
442 template<typename V, typename T, typename P>
443 auto abs(const ThermoVectorBase<V,T,P>& l) -> ThermoVector
444 {
445  const Vector tmp1 = abs(l.val);
446  const Vector tmp2 = l.val/tmp1;
447  return {tmp1, diag(tmp2) * l.ddT, diag(tmp2) * l.ddP};
448 }
449 
450 template<typename V, typename T, typename P>
451 auto sqrt(const ThermoVectorBase<V,T,P>& l) -> ThermoVector
452 {
453  const Vector tmp1 = sqrt(l.val);
454  const Vector tmp2 = 0.5 * tmp1/l.val;
455  return {tmp1, diag(tmp2) * l.ddT, diag(tmp2) * l.ddP};
456 }
457 
458 template<typename V, typename T, typename P>
459 auto pow(const ThermoVectorBase<V,T,P>& l, double power) -> ThermoVector
460 {
461  const Vector tmp1 = pow(l.val, power);
462  const Vector tmp2 = power * tmp1/l.val;
463  return {tmp1, diag(tmp2) * l.ddT, diag(tmp2) * l.ddP};
464 }
465 
466 template<typename V, typename T, typename P>
467 auto exp(const ThermoVectorBase<V,T,P>& l) -> ThermoVector
468 {
469  const Vector tmp = exp(l.val);
470  return {tmp, diag(tmp) * l.ddT, diag(tmp) * l.ddP};
471 }
472 
473 template<typename V, typename T, typename P>
474 auto log(const ThermoVectorBase<V,T,P>& l) -> ThermoVector
475 {
476  const Vector tmp1 = log(l.val);
477  const Vector tmp2 = 1.0/l.val;
478  return {tmp1, diag(tmp2) * l.ddT, diag(tmp2) * l.ddP};
479 }
480 
481 template<typename V, typename T, typename P>
482 auto log10(const ThermoVectorBase<V,T,P>& l) -> ThermoVector
483 {
484  const double log10e = 0.4342944819032518166679324;
485  const Vector tmp1 = log10e*log(l.val);
486  const Vector tmp2 = log10e/l.val;
487  return {tmp1, diag(tmp2) * l.ddT, diag(tmp2) * l.ddP, diag(tmp2) * l.ddn};
488 }
489 
490 } // namespace Reaktoro
auto row(ChemicalVectorBase< V, T, P, N > &vec, Index irow) -> ChemicalScalarBase< double &, decltype(tr(vec.ddn.row(irow)))>
Return a ChemicalScalarBase with reference to the chemical scalar in a given row. ...
Definition: ChemicalVector.hpp:590
ThermoVectorBase(Index nrows, double val)
Construct a ThermoVectorBase instance with given number of rows and value.
Definition: ThermoVector.hpp:65
V val
The vector of values of the thermodynamic properties.
Definition: ThermoVector.hpp:45
V val
The value of the thermodynamic property.
Definition: ThermoScalar.hpp:49
auto operator*=(const ThermoVectorBase< VR, TR, PR > &other) -> ThermoVectorBase &
Assign-multiplication of a ThermoVectorBase instance.
Definition: ThermoVector.hpp:198
ThermoVectorBase< Vector, Vector, Vector > ThermoVector
A type that defines a vector thermo property.
Definition: ScalarTypes.hpp:41
auto operator[](Index irow) -> ThermoScalarBase< double & >
Return a ChemicalScalarBase with reference to the thermo scalar in a given row.
Definition: ThermoVector.hpp:223
auto fill(double value) -> void
Assign a scalarsto this.
Definition: ThermoVector.hpp:106
std::vector< Index > Indices
Define a type that represents a collection of indices.
Definition: Index.hpp:29
auto operator*=(double scalar) -> ThermoVectorBase &
Assign-multiplication of a ThermoVectorBase instance.
Definition: ThermoVector.hpp:207
V ddP
The partial pressure derivative of the thermodynamic property.
Definition: ThermoScalar.hpp:55
auto operator+=(double scalar) -> ThermoVectorBase &
Assign-addition of a scalar.
Definition: ThermoVector.hpp:163
ThermoVectorBase()
Construct a default ThermoVectorBase instance.
Definition: ThermoVector.hpp:54
ThermoVectorBase(const ThermoVectorBase< VR, TR, PR > &other)
Construct a ChemicalVectorBase instance from another.
Definition: ThermoVector.hpp:77
auto diag(const Eigen::MatrixBase< Derived > &vec) -> decltype(vec.asDiagonal())
Return a diagonal matrix representation of a vector.
Definition: Matrix.hxx:175
T ddT
The vector of partial temperature derivatives of the thermodynamic properties.
Definition: ThermoVector.hpp:48
P ddP
The vector of partial pressure derivatives of the thermodynamic properties.
Definition: ThermoVector.hpp:51
V ddT
The partial temperature derivative of the thermodynamic property.
Definition: ThermoScalar.hpp:52
auto rows(ChemicalVectorBase< V, T, P, N > &vec, Index irow, Index nrows) -> ChemicalVectorBase< decltype(rows(vec.val, irow, nrows)), decltype(rows(vec.ddT, irow, nrows)), decltype(rows(vec.ddP, irow, nrows)), decltype(rows(vec.ddn, irow, nrows))>
Return a reference of a sequence of rows of this ChemicalVectorBase instance.
Definition: ChemicalVector.hpp:618
auto operator-=(const ThermoVectorBase< VR, TR, PR > &other) -> ThermoVectorBase &
Assign-subtraction of a ThermoVectorBase instance.
Definition: ThermoVector.hpp:171
auto size() const -> Index
Return the number of rows in this ChemicalVectorBase instance.
Definition: ThermoVector.hpp:82
auto operator[](Index irow) const -> ThermoScalarBase< const double & >
Return a ChemicalScalarBase with const reference to the thermo scalar in a given row.
Definition: ThermoVector.hpp:229
auto operator+=(const ThermoVectorBase< VR, TR, PR > &other) -> ThermoVectorBase &
Assign-addition of a ThermoVectorBase instance.
Definition: ThermoVector.hpp:144
auto operator+=(const ThermoScalarBase< VR > &other) -> ThermoVectorBase &
Assign-addition of a ThermoScalarBase instance.
Definition: ThermoVector.hpp:154
Eigen::VectorXd Vector
Define an alias to the vector type of the Eigen library.
Definition: Matrix.hpp:384
A template base class to represent a thermodynamic scalar and its partial derivatives.
Definition: ScalarTypes.hpp:29
auto operator=(const ThermoVectorBase< VR, TR, PR > &other) -> ThermoVectorBase &
Assign a ThermoVectorBase instance to this ThermoVectorBase instance.
Definition: ThermoVector.hpp:115
A template base class to represent a vector of thermodynamic scalars and their partial derivatives...
Definition: ScalarTypes.hpp:32
auto fill(const ThermoScalarBase< VR > &other) -> void
Assign a ThermoScalarBase instance to this.
Definition: ThermoVector.hpp:98
auto operator/=(double scalar) -> ThermoVectorBase &
Assign-division of a ThermoVectorBase instance.
Definition: ThermoVector.hpp:216
ThermoVectorBase(const V &val, const T &ddT, const P &ddP)
Construct a ChemicalVectorBase instance with given values and derivatives.
Definition: ThermoVector.hpp:72
std::size_t Index
Define a type that represents an index.
Definition: Index.hpp:26
auto operator=(const ThermoScalarBase< VR > &other) -> ThermoVectorBase &
Assign a ThermoScalarBase instance to this ThermoVectorBase instance.
Definition: ThermoVector.hpp:125
auto resize(Index nrows) -> void
Resize this ChemicalVectorBase instance with new number of rows.
Definition: ThermoVector.hpp:89
auto operator=(double other) -> ThermoVectorBase &
Assign a scalar to this ThermoVectorBase instance.
Definition: ThermoVector.hpp:134
The namespace containing all components of the Reaktoro library.
Definition: ChemicalScalar.hpp:24
ThermoVectorBase(Index nrows)
Construct a ThermoVectorBase instance with given number of rows.
Definition: ThermoVector.hpp:59
auto operator-=(double other) -> ThermoVectorBase &
Assign-subtraction of a scalar.
Definition: ThermoVector.hpp:190
auto operator-=(const ThermoScalarBase< VR > &other) -> ThermoVectorBase &
Assign-subtraction of a ThermoVectorBase instance.
Definition: ThermoVector.hpp:181
auto zeros(Index rows) -> decltype(Vector::Zero(rows))
Return an expression of a zero vector.
Definition: Matrix.hxx:22