mirror of
https://github.com/Cantera/cantera.git
synced 2025-02-25 18:55:29 -06:00
Implement Steady-State, Adiabatic, Cylindrical (by default) Plug-Flow Reactor & Add python example
This commit is contained in:
Nick Curtis
Ingmar Schoegl
parent
5973dbd88b
commit
568cf48233
@ -2241,7 +2241,7 @@ if addInstallActions:
|
||||
|
||||
|
||||
if env['system_sundials'] == 'y':
|
||||
env['sundials_libs'] = ['sundials_cvodes', 'sundials_ida', 'sundials_nvecserial']
|
||||
env['sundials_libs'] = ['sundials_cvodes', 'sundials_idas', 'sundials_nvecserial']
|
||||
if env['use_lapack']:
|
||||
if env.get('has_sundials_lapack'):
|
||||
env['sundials_libs'].extend(('sundials_sunlinsollapackdense',
|
||||
|
||||
523
data/SiF4_NH3_mec.yaml
Normal file
523
data/SiF4_NH3_mec.yaml
Normal file
@ -0,0 +1,523 @@
|
||||
description: |-
|
||||
nits(length='cm', time='s', quantity='mol', act_energy='cal/mol')
|
||||
|
||||
generator: cti2yaml
|
||||
cantera-version: 3.0.0a5
|
||||
date: Sun, 23 Apr 2023 21:03:15 -0400
|
||||
input-files: [SiF4_NH3_mec.cti]
|
||||
|
||||
units: {length: cm, quantity: mol, activation-energy: cal/mol}
|
||||
|
||||
phases:
|
||||
- name: SI3N4
|
||||
thermo: ideal-surface
|
||||
elements: [H, N, Si, F]
|
||||
species: [HN_SIF(S), HN_NH2(S), F3SI_NH2(S), F2SINH(S), H2NFSINH(S), HN(FSINH)2(S)]
|
||||
kinetics: surface
|
||||
reactions: [SI3N4-reactions]
|
||||
state:
|
||||
T: 300.0
|
||||
P: 1.01325e+05
|
||||
adjacent-phases: [gas, SiBulk, NBulk]
|
||||
site-density: 4.1683e-09
|
||||
- name: gas
|
||||
thermo: ideal-gas
|
||||
elements: [H, N, Si, F]
|
||||
species: [H2, H, N2, N, NH, NH2, NNH, N2H2, N2H3, N2H4, HF, F, SIF4, SIF3,
|
||||
SIHF3, SIF3NH2, NH3]
|
||||
kinetics: gas
|
||||
reactions: [gas-reactions]
|
||||
transport: mixture-averaged
|
||||
state:
|
||||
T: 300.0
|
||||
P: 1.01325e+05
|
||||
- name: SiBulk
|
||||
thermo: fixed-stoichiometry
|
||||
elements: [Si]
|
||||
species: [SI(D)]
|
||||
- name: NBulk
|
||||
thermo: fixed-stoichiometry
|
||||
elements: [N]
|
||||
species: [N(D)]
|
||||
|
||||
species:
|
||||
- name: H2
|
||||
composition: {H: 2}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [200.0, 1000.0, 6000.0]
|
||||
data:
|
||||
- [2.34433112, 7.98052075e-03, -1.9478151e-05, 2.01572094e-08, -7.37611761e-12,
|
||||
-917.935173, 0.683010238]
|
||||
- [2.93286575, 8.26608026e-04, -1.46402364e-07, 1.54100414e-11, -6.888048e-16,
|
||||
-813.065581, -1.02432865]
|
||||
transport:
|
||||
model: gas
|
||||
geometry: linear
|
||||
diameter: 2.92
|
||||
well-depth: 38.0
|
||||
polarizability: 0.79
|
||||
rotational-relaxation: 280.0
|
||||
note: TPIS78
|
||||
- name: H
|
||||
composition: {H: 1}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [200.0, 1000.0, 6000.0]
|
||||
data:
|
||||
- [2.5, 0.0, 0.0, 0.0, 0.0, 2.547366e+04, -0.44668285]
|
||||
- [2.5, 0.0, 0.0, 0.0, 0.0, 2.547366e+04, -0.44668285]
|
||||
transport:
|
||||
model: gas
|
||||
geometry: atom
|
||||
diameter: 2.05
|
||||
well-depth: 145.0
|
||||
note: L6/94
|
||||
- name: N2
|
||||
composition: {N: 2}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [200.0, 1000.0, 6000.0]
|
||||
data:
|
||||
- [3.53100528, -1.23660988e-04, -5.02999433e-07, 2.43530612e-09, -1.40881235e-12,
|
||||
-1046.97628, 2.96747038]
|
||||
- [2.95257637, 1.3969004e-03, -4.92631603e-07, 7.86010195e-11, -4.60755204e-15,
|
||||
-923.948688, 5.87188762]
|
||||
transport:
|
||||
model: gas
|
||||
geometry: linear
|
||||
diameter: 3.621
|
||||
well-depth: 97.53
|
||||
polarizability: 1.76
|
||||
rotational-relaxation: 4.0
|
||||
note: G8/02
|
||||
- name: N
|
||||
composition: {N: 1}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 5000.0]
|
||||
data:
|
||||
- [2.503071, -2.180018e-05, 5.420529e-08, -5.64756e-11, 2.099904e-14,
|
||||
5.60989e+04, 4.167566]
|
||||
- [2.450268, 1.066146e-04, -7.465337e-08, 1.879652e-11, -1.025984e-15,
|
||||
5.611604e+04, 4.448758]
|
||||
transport:
|
||||
model: gas
|
||||
geometry: atom
|
||||
diameter: 3.298
|
||||
well-depth: 71.4
|
||||
note: '120186'
|
||||
- name: NH
|
||||
composition: {H: 1, N: 1}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [200.0, 1000.0, 6000.0]
|
||||
data:
|
||||
- [3.4929084, 3.1179197e-04, -1.4890484e-06, 2.4816442e-09, -1.0356967e-12,
|
||||
4.1894294e+04, 1.8483277]
|
||||
- [2.7836929, 1.3298429e-03, -4.2478047e-07, 7.8348504e-11, -5.504447e-15,
|
||||
4.2134514e+04, 5.7407798]
|
||||
transport:
|
||||
model: gas
|
||||
geometry: linear
|
||||
diameter: 2.65
|
||||
well-depth: 80.0
|
||||
rotational-relaxation: 4.0
|
||||
- name: NH2
|
||||
composition: {H: 2, N: 1}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 5000.0]
|
||||
data:
|
||||
- [3.432493, 3.29954e-03, -6.6136e-06, 8.590947e-09, -3.572047e-12,
|
||||
2.177228e+04, 3.090111]
|
||||
- [2.961311, 2.932699e-03, -9.0636e-07, 1.617257e-10, -1.2042e-14, 2.191977e+04,
|
||||
5.777878]
|
||||
transport:
|
||||
model: gas
|
||||
geometry: nonlinear
|
||||
diameter: 2.65
|
||||
well-depth: 80.0
|
||||
polarizability: 2.26
|
||||
rotational-relaxation: 4.0
|
||||
note: '121686'
|
||||
- name: NNH
|
||||
composition: {H: 1, N: 2}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [250.0, 1000.0, 4000.0]
|
||||
data:
|
||||
- [3.501344, 2.053587e-03, 7.17041e-07, 4.921348e-10, -9.67117e-13,
|
||||
2.833347e+04, 6.391837]
|
||||
- [4.415342, 1.614388e-03, -1.632894e-07, -8.559846e-11, 1.614791e-14,
|
||||
2.788029e+04, 0.9042888]
|
||||
transport:
|
||||
model: gas
|
||||
geometry: nonlinear
|
||||
diameter: 3.798
|
||||
well-depth: 71.4
|
||||
rotational-relaxation: 1.0
|
||||
note: '120186'
|
||||
- name: N2H2
|
||||
composition: {H: 2, N: 2}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 5000.0]
|
||||
data:
|
||||
- [1.617999, 0.01306312, -1.715712e-05, 1.605608e-08, -6.093639e-12,
|
||||
2.467526e+04, 13.79467]
|
||||
- [3.371185, 6.039968e-03, -2.303854e-06, 4.062789e-10, -2.713144e-14,
|
||||
2.418172e+04, 4.980585]
|
||||
transport:
|
||||
model: gas
|
||||
geometry: nonlinear
|
||||
diameter: 3.798
|
||||
well-depth: 71.4
|
||||
rotational-relaxation: 1.0
|
||||
note: '121286'
|
||||
- name: N2H3
|
||||
composition: {H: 3, N: 2}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 5000.0]
|
||||
data:
|
||||
- [3.174204, 4.715907e-03, 1.334867e-05, -1.919685e-08, 7.487564e-12,
|
||||
1.72727e+04, 7.557224]
|
||||
- [4.441846, 7.214271e-03, -2.495684e-06, 3.920565e-10, -2.29895e-14,
|
||||
1.664221e+04, -0.4275205]
|
||||
transport:
|
||||
model: gas
|
||||
geometry: nonlinear
|
||||
diameter: 3.9
|
||||
well-depth: 200.0
|
||||
rotational-relaxation: 1.0
|
||||
note: '120186'
|
||||
- name: N2H4
|
||||
composition: {H: 4, N: 2}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 5000.0]
|
||||
data:
|
||||
- [0.06442606, 0.0274973, -2.899451e-05, 1.74524e-08, -4.422282e-12,
|
||||
1.045192e+04, 21.27789]
|
||||
- [4.977317, 9.595519e-03, -3.547639e-06, 6.124299e-10, -4.029795e-14,
|
||||
9341.219, -2.96299]
|
||||
transport:
|
||||
model: gas
|
||||
geometry: nonlinear
|
||||
diameter: 4.23
|
||||
well-depth: 205.0
|
||||
polarizability: 4.26
|
||||
rotational-relaxation: 1.5
|
||||
note: '121286'
|
||||
- name: HF
|
||||
composition: {F: 1, H: 1}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 5000.0]
|
||||
data:
|
||||
- [3.4379986, 5.3571598e-04, -1.5229655e-06, 1.7564491e-09, -5.786994e-13,
|
||||
-3.3818972e+04, 1.1930153]
|
||||
- [2.991911, 7.1489475e-04, -6.8630973e-08, -1.161713e-11, 1.9412375e-15,
|
||||
-3.3621364e+04, 3.8123288]
|
||||
transport:
|
||||
model: gas
|
||||
geometry: linear
|
||||
diameter: 3.148
|
||||
well-depth: 330.0
|
||||
dipole: 1.92
|
||||
polarizability: 2.46
|
||||
rotational-relaxation: 1.0
|
||||
note: J6/77
|
||||
- name: F
|
||||
composition: {F: 1}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 5000.0]
|
||||
data:
|
||||
- [2.812874, -3.3023098e-06, -1.289731e-06, 1.6837365e-09, -6.4587833e-13,
|
||||
8660.4019, 3.0984198]
|
||||
- [2.7004353, -2.2293182e-04, 9.7941385e-08, -1.9123038e-11, 1.3768154e-15,
|
||||
8716.3617, 3.8067182]
|
||||
transport:
|
||||
model: gas
|
||||
geometry: atom
|
||||
diameter: 2.75
|
||||
well-depth: 80.0
|
||||
note: J9/65
|
||||
- name: SIF4
|
||||
composition: {F: 4, Si: 1}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 5000.0]
|
||||
data:
|
||||
- [2.1893068, 0.033702007, -4.6723179e-05, 3.1584638e-08, -8.4506114e-12,
|
||||
-1.9603289e+05, 13.287308]
|
||||
- [10.478473, 2.8586756e-03, -1.2646314e-06, 2.4746863e-10, -1.7824296e-14,
|
||||
-1.979055e+05, -27.520641]
|
||||
transport:
|
||||
model: gas
|
||||
geometry: nonlinear
|
||||
diameter: 4.88
|
||||
well-depth: 171.9
|
||||
rotational-relaxation: 1.0
|
||||
note: J6/76
|
||||
- name: SIF3
|
||||
composition: {F: 3, Si: 1}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 3000.0]
|
||||
data:
|
||||
- [4.6628685, 0.010087878, -1.8055442e-06, -7.769299e-09, 4.3778518e-12,
|
||||
-1.2129652e+05, 4.672966]
|
||||
- [8.5247898, 1.3237924e-03, -2.1042787e-07, -1.149504e-10, 3.0553014e-14,
|
||||
-1.2235223e+05, -15.502343]
|
||||
transport:
|
||||
model: gas
|
||||
geometry: nonlinear
|
||||
diameter: 4.359
|
||||
well-depth: 309.6
|
||||
rotational-relaxation: 1.0
|
||||
note: '41889'
|
||||
- name: SIHF3
|
||||
composition: {F: 3, H: 1, Si: 1}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 3000.0]
|
||||
data:
|
||||
- [3.9180529, 0.014639172, -1.8560698e-06, -1.0582003e-08, 5.6175433e-12,
|
||||
-1.4704386e+05, 7.0242615]
|
||||
- [9.3635674, 2.9475559e-03, -3.577633e-07, -2.8582245e-10, 6.9157286e-14,
|
||||
-1.4860736e+05, -21.694529]
|
||||
transport:
|
||||
model: gas
|
||||
geometry: nonlinear
|
||||
diameter: 4.681
|
||||
well-depth: 180.8
|
||||
rotational-relaxation: 1.0
|
||||
note: '41889'
|
||||
- name: SIF3NH2
|
||||
composition: {F: 3, H: 2, N: 1, Si: 1}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 3000.0]
|
||||
data:
|
||||
- [6.229403, 0.017780151, -2.6123043e-06, -1.2672435e-08, 7.0445559e-12,
|
||||
-1.6258489e+05, 0.20454407]
|
||||
- [12.109636, 4.3832823e-03, -4.1422453e-07, -3.9890902e-10, 8.9589543e-14,
|
||||
-1.6417678e+05, -30.469284]
|
||||
transport:
|
||||
model: gas
|
||||
geometry: nonlinear
|
||||
diameter: 4.975
|
||||
well-depth: 231.0
|
||||
rotational-relaxation: 1.0
|
||||
note: '41889'
|
||||
- name: NH3
|
||||
composition: {H: 3, N: 1}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 5000.0]
|
||||
data:
|
||||
- [2.204352, 0.01011476, -1.465265e-05, 1.447235e-08, -5.328509e-12,
|
||||
-6525.488, 8.127138]
|
||||
- [2.461904, 6.059166e-03, -2.004977e-06, 3.136003e-10, -1.938317e-14,
|
||||
-6493.27, 7.472097]
|
||||
transport:
|
||||
model: gas
|
||||
geometry: nonlinear
|
||||
diameter: 2.92
|
||||
well-depth: 481.0
|
||||
dipole: 1.47
|
||||
rotational-relaxation: 10.0
|
||||
note: '121386'
|
||||
- name: F3SI_NH2(S)
|
||||
composition: {F: 3, H: 2, N: 1, Si: 1}
|
||||
sites: 2.0
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 1685.0]
|
||||
data:
|
||||
- [0.84197538, 8.3710416e-03, -1.307703e-05, 9.7593603e-09, -2.727938e-12,
|
||||
-524.86288, -4.5272678]
|
||||
- [2.4753989, 8.8112187e-04, -2.0939481e-07, 4.2757187e-12, 1.6006564e-14,
|
||||
-812.5562, -12.188747]
|
||||
note: J3/67
|
||||
- name: F2SINH(S)
|
||||
composition: {F: 2, H: 1, N: 1, Si: 1}
|
||||
sites: 2.0
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 1685.0]
|
||||
data:
|
||||
- [0.84197538, 8.3710416e-03, -1.307703e-05, 9.7593603e-09, -2.727938e-12,
|
||||
-524.86288, -4.5272678]
|
||||
- [2.4753989, 8.8112187e-04, -2.0939481e-07, 4.2757187e-12, 1.6006564e-14,
|
||||
-812.5562, -12.188747]
|
||||
note: J3/67
|
||||
- name: H2NFSINH(S)
|
||||
composition: {F: 1, H: 3, N: 2, Si: 1}
|
||||
sites: 2.0
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 1685.0]
|
||||
data:
|
||||
- [0.84197538, 8.3710416e-03, -1.307703e-05, 9.7593603e-09, -2.727938e-12,
|
||||
-524.86288, -4.5272678]
|
||||
- [2.4753989, 8.8112187e-04, -2.0939481e-07, 4.2757187e-12, 1.6006564e-14,
|
||||
-812.5562, -12.188747]
|
||||
note: J3/67
|
||||
- name: HN(FSINH)2(S)
|
||||
composition: {F: 2, H: 3, N: 3, Si: 2}
|
||||
sites: 4.0
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 1685.0]
|
||||
data:
|
||||
- [0.84197538, 8.3710416e-03, -1.307703e-05, 9.7593603e-09, -2.727938e-12,
|
||||
-524.86288, -4.5272678]
|
||||
- [2.4753989, 8.8112187e-04, -2.0939481e-07, 4.2757187e-12, 1.6006564e-14,
|
||||
-812.5562, -12.188747]
|
||||
note: J3/67
|
||||
- name: HN_SIF(S)
|
||||
composition: {F: 1, H: 1, N: 1, Si: 1}
|
||||
sites: 2.0
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 1685.0]
|
||||
data:
|
||||
- [0.84197538, 8.3710416e-03, -1.307703e-05, 9.7593603e-09, -2.727938e-12,
|
||||
-524.86288, -4.5272678]
|
||||
- [2.4753989, 8.8112187e-04, -2.0939481e-07, 4.2757187e-12, 1.6006564e-14,
|
||||
-812.5562, -12.188747]
|
||||
note: J3/67
|
||||
- name: HN_NH2(S)
|
||||
composition: {H: 3, N: 2}
|
||||
sites: 2.0
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 1685.0]
|
||||
data:
|
||||
- [0.84197538, 8.3710416e-03, -1.307703e-05, 9.7593603e-09, -2.727938e-12,
|
||||
-524.86288, -4.5272678]
|
||||
- [2.4753989, 8.8112187e-04, -2.0939481e-07, 4.2757187e-12, 1.6006564e-14,
|
||||
-812.5562, -12.188747]
|
||||
note: J3/67
|
||||
- name: SI(D)
|
||||
composition: {Si: 1}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 1685.0]
|
||||
data:
|
||||
- [0.84197538, 8.3710416e-03, -1.307703e-05, 9.7593603e-09, -2.727938e-12,
|
||||
-524.86288, -4.5272678]
|
||||
- [2.4753989, 8.8112187e-04, -2.0939481e-07, 4.2757187e-12, 1.6006564e-14,
|
||||
-812.5562, -12.188747]
|
||||
equation-of-state:
|
||||
model: constant-volume
|
||||
density: 2.066 g/cm^3
|
||||
note: J3/67
|
||||
- name: N(D)
|
||||
composition: {N: 1}
|
||||
thermo:
|
||||
model: NASA7
|
||||
temperature-ranges: [300.0, 1000.0, 1685.0]
|
||||
data:
|
||||
- [0.84197538, 8.3710416e-03, -1.307703e-05, 9.7593603e-09, -2.727938e-12,
|
||||
-524.86288, -4.5272678]
|
||||
- [2.4753989, 8.8112187e-04, -2.0939481e-07, 4.2757187e-12, 1.6006564e-14,
|
||||
-812.5562, -12.188747]
|
||||
equation-of-state:
|
||||
model: constant-volume
|
||||
density: 1.374 g/cm^3
|
||||
note: J3/67
|
||||
|
||||
gas-reactions:
|
||||
- equation: H + H + M <=> H2 + M # Reaction 1
|
||||
type: three-body
|
||||
rate-constant: {A: 1.0e+18, b: -1.0, Ea: 0.0}
|
||||
efficiencies: {H2: 0.0}
|
||||
- equation: H + H + H2 <=> H2 + H2 # Reaction 2
|
||||
rate-constant: {A: 9.2e+16, b: -0.6, Ea: 0.0}
|
||||
- equation: NH + N <=> N2 + H # Reaction 3
|
||||
rate-constant: {A: 3.0e+13, b: 0.0, Ea: 0.0}
|
||||
- equation: NH + H <=> N + H2 # Reaction 4
|
||||
rate-constant: {A: 1.0e+14, b: 0.0, Ea: 0.0}
|
||||
- equation: NH2 + H <=> NH + H2 # Reaction 5
|
||||
rate-constant: {A: 6.92e+13, b: 0.0, Ea: 3650.0}
|
||||
- equation: NH3 + H <=> NH2 + H2 # Reaction 6
|
||||
rate-constant: {A: 6.36e+05, b: 2.39, Ea: 1.0171e+04}
|
||||
- equation: NNH <=> N2 + H # Reaction 7
|
||||
rate-constant: {A: 1.0e+04, b: 0.0, Ea: 0.0}
|
||||
- equation: NNH + H <=> N2 + H2 # Reaction 8
|
||||
rate-constant: {A: 1.0e+14, b: 0.0, Ea: 0.0}
|
||||
- equation: NNH + NH2 <=> N2 + NH3 # Reaction 9
|
||||
rate-constant: {A: 5.0e+13, b: 0.0, Ea: 0.0}
|
||||
- equation: NNH + NH <=> N2 + NH2 # Reaction 10
|
||||
rate-constant: {A: 5.0e+13, b: 0.0, Ea: 0.0}
|
||||
- equation: NH2 + NH <=> N2H2 + H # Reaction 11
|
||||
rate-constant: {A: 5.0e+13, b: 0.0, Ea: 0.0}
|
||||
- equation: NH + NH <=> N2 + H + H # Reaction 12
|
||||
rate-constant: {A: 2.54e+13, b: 0.0, Ea: 0.0}
|
||||
- equation: NH2 + N <=> N2 + H + H # Reaction 13
|
||||
rate-constant: {A: 7.2e+13, b: 0.0, Ea: 0.0}
|
||||
- equation: N2H2 + M <=> NNH + H + M # Reaction 14
|
||||
type: three-body
|
||||
rate-constant: {A: 5.0e+16, b: 0.0, Ea: 5.0e+04}
|
||||
efficiencies: {H2: 2.0, N2: 2.0}
|
||||
- equation: N2H2 + H <=> NNH + H2 # Reaction 15
|
||||
rate-constant: {A: 5.0e+13, b: 0.0, Ea: 1000.0}
|
||||
- equation: N2H2 + NH <=> NNH + NH2 # Reaction 16
|
||||
rate-constant: {A: 1.0e+13, b: 0.0, Ea: 1000.0}
|
||||
- equation: N2H2 + NH2 <=> NH3 + NNH # Reaction 17
|
||||
rate-constant: {A: 1.0e+13, b: 0.0, Ea: 1000.0}
|
||||
- equation: NH2 + NH2 <=> N2H2 + H2 # Reaction 18
|
||||
rate-constant: {A: 5.0e+11, b: 0.0, Ea: 0.0}
|
||||
- equation: NH3 + M <=> NH2 + H + M # Reaction 19
|
||||
type: three-body
|
||||
rate-constant: {A: 1.4e+16, b: 0.0, Ea: 9.06e+04}
|
||||
- equation: N2H3 + H <=> NH2 + NH2 # Reaction 20
|
||||
rate-constant: {A: 1.6e+12, b: 0.0, Ea: 0.0}
|
||||
- equation: N2H3 + M <=> N2H2 + H + M # Reaction 21
|
||||
type: three-body
|
||||
rate-constant: {A: 3.5e+16, b: 0.0, Ea: 4.6e+04}
|
||||
- equation: N2H3 + NH <=> NH2 + N2H2 # Reaction 22
|
||||
rate-constant: {A: 2.0e+13, b: 0.0, Ea: 0.0}
|
||||
- equation: NH2 + NH2 + M <=> N2H4 + M # Reaction 23
|
||||
type: three-body
|
||||
rate-constant: {A: 3.0e+20, b: -1.0, Ea: 0.0}
|
||||
- equation: H + N2H4 <=> H2 + N2H3 # Reaction 24
|
||||
rate-constant: {A: 1.3e+13, b: 0.0, Ea: 2500.0}
|
||||
- equation: NH2 + N2H4 <=> NH3 + N2H3 # Reaction 25
|
||||
rate-constant: {A: 3.9e+12, b: 0.0, Ea: 1500.0}
|
||||
- equation: NH + H + M <=> NH2 + M # Reaction 26
|
||||
type: three-body
|
||||
rate-constant: {A: 2.0e+16, b: -0.5, Ea: 0.0}
|
||||
- equation: NH2 + NH2 <=> NH3 + NH # Reaction 27
|
||||
rate-constant: {A: 5.0e+12, b: 0.0, Ea: 1.0e+04}
|
||||
- equation: F + NH3 <=> NH2 + HF # Reaction 28
|
||||
rate-constant: {A: 4.27e+11, b: 0.5, Ea: 800.0}
|
||||
- equation: SIF4 <=> SIF3 + F # Reaction 29
|
||||
rate-constant: {A: 3.0e+12, b: 0.0, Ea: 1.4717e+05}
|
||||
- equation: H + SIF4 <=> HF + SIF3 # Reaction 30
|
||||
rate-constant: {A: 1.0e+13, b: 0.0, Ea: 5.0e+04}
|
||||
- equation: NH2 + SIF4 <=> SIF3NH2 + F # Reaction 31
|
||||
rate-constant: {A: 1.0e+11, b: 0.0, Ea: 4.095e+04}
|
||||
- equation: NH3 + SIF3 <=> SIF3NH2 + H # Reaction 32
|
||||
rate-constant: {A: 1.0e+11, b: 0.0, Ea: 5000.0}
|
||||
- equation: NH3 + SIF3 <=> SIHF3 + NH2 # Reaction 33
|
||||
rate-constant: {A: 1.0e+11, b: 0.0, Ea: 1.0e+04}
|
||||
|
||||
SI3N4-reactions:
|
||||
- equation: F3SI_NH2(S) => F2SINH(S) + HF # Reaction 34
|
||||
rate-constant: {A: 1.0e+05, b: 0.0, Ea: 0.0}
|
||||
- equation: NH3 + F2SINH(S) => H2NFSINH(S) + HF # Reaction 35
|
||||
rate-constant: {A: 7.562e+08, b: 0.5, Ea: 0.0}
|
||||
- equation: H2NFSINH(S) + F2SINH(S) => HN(FSINH)2(S) + HF # Reaction 36
|
||||
rate-constant: {A: 1.0e+15, b: 0.0, Ea: 0.0}
|
||||
- equation: NH3 + HN_SIF(S) => HN_NH2(S) + SI(D) + HF # Reaction 37
|
||||
rate-constant: {A: 7.56e+08, b: 0.5, Ea: 0.0}
|
||||
- equation: SIF4 + HN_NH2(S) => F3SI_NH2(S) + N(D) + HF # Reaction 38
|
||||
rate-constant: {A: 3.1e+08, b: 0.5, Ea: 0.0}
|
||||
- equation: HN(FSINH)2(S) + F2SINH(S) => 3 HN_SIF(S) + N(D) + HF # Reaction 39
|
||||
rate-constant: {A: 1.0e+15, b: 0.0, Ea: 0.0}
|
||||
@ -73,7 +73,7 @@ if env['system_sundials'] == 'n':
|
||||
)
|
||||
|
||||
# Copy sundials header files into common include directory
|
||||
for subdir in ('sundials', 'nvector', 'cvodes', 'ida', 'sunmatrix',
|
||||
for subdir in ('sundials', 'nvector', 'cvodes', 'idas', 'sunmatrix',
|
||||
'sunlinsol', 'sunnonlinsol'):
|
||||
for header in multi_glob(env, 'sundials/include/'+subdir, 'h'):
|
||||
ext_copies.extend(
|
||||
@ -84,7 +84,7 @@ if env['system_sundials'] == 'n':
|
||||
|
||||
# Compile Sundials source files. Skip files related to the Sundials Fortran
|
||||
# interface, which start with 'fsun'.
|
||||
subdirs = ['sundials', 'nvector/serial', 'cvodes', 'ida', 'sunmatrix/band',
|
||||
subdirs = ['sundials', 'nvector/serial', 'cvodes', 'idas', 'sunmatrix/band',
|
||||
'sunmatrix/dense', 'sunmatrix/sparse', 'sunlinsol/dense',
|
||||
'sunlinsol/band','sunlinsol/spgmr', 'sunnonlinsol/newton']
|
||||
if env['use_lapack']:
|
||||
|
||||
@ -1,270 +0,0 @@
|
||||
/**
|
||||
* @file DAE_Solver.h
|
||||
* Header file for class DAE_Solver
|
||||
*/
|
||||
|
||||
// This file is part of Cantera. See License.txt in the top-level directory or
|
||||
// at https://cantera.org/license.txt for license and copyright information.
|
||||
|
||||
#ifndef CT_DAE_Solver_H
|
||||
#define CT_DAE_Solver_H
|
||||
|
||||
#include "ResidJacEval.h"
|
||||
#include "cantera/base/global.h"
|
||||
|
||||
namespace Cantera
|
||||
{
|
||||
|
||||
class Jacobian
|
||||
{
|
||||
public:
|
||||
Jacobian() = default;
|
||||
virtual ~Jacobian() = default;
|
||||
virtual bool supplied() {
|
||||
return false;
|
||||
}
|
||||
virtual bool isBanded() {
|
||||
return false;
|
||||
}
|
||||
virtual int lowerBandWidth() {
|
||||
return 0;
|
||||
}
|
||||
virtual int upperBandWidth() {
|
||||
return 0;
|
||||
}
|
||||
};
|
||||
|
||||
class BandedJacobian : public Jacobian
|
||||
{
|
||||
public:
|
||||
BandedJacobian(int ml, int mu) {
|
||||
m_ml = ml;
|
||||
m_mu = mu;
|
||||
}
|
||||
virtual bool supplied() {
|
||||
return false;
|
||||
}
|
||||
virtual bool isBanded() {
|
||||
return true;
|
||||
}
|
||||
virtual int lowerBandWidth() {
|
||||
return m_ml;
|
||||
}
|
||||
virtual int upperBandWidth() {
|
||||
return m_mu;
|
||||
}
|
||||
protected:
|
||||
int m_ml, m_mu;
|
||||
};
|
||||
|
||||
const int cDirect = 0;
|
||||
const int cKrylov = 1;
|
||||
|
||||
|
||||
/**
|
||||
* Wrapper for DAE solvers
|
||||
*
|
||||
* @attention This class currently does not have any test cases or examples. Its
|
||||
* implementation may be incomplete, and future changes to Cantera may
|
||||
* unexpectedly cause this class to stop working. If you use this class,
|
||||
* please consider contributing examples or test cases. In the absence of
|
||||
* new tests or examples, this class may be deprecated and removed in a
|
||||
* future version of Cantera. See
|
||||
* https://github.com/Cantera/cantera/issues/267 for additional information.
|
||||
*/
|
||||
class DAE_Solver
|
||||
{
|
||||
public:
|
||||
DAE_Solver(ResidJacEval& f) :
|
||||
m_resid(f),
|
||||
m_neq(f.nEquations())
|
||||
{
|
||||
}
|
||||
|
||||
virtual ~DAE_Solver() = default;
|
||||
|
||||
/**
|
||||
* Set error tolerances. This version specifies a scalar
|
||||
* relative tolerance, and a vector absolute tolerance.
|
||||
*/
|
||||
virtual void setTolerances(doublereal reltol,
|
||||
doublereal* abstol) {
|
||||
warn("setTolerances");
|
||||
}
|
||||
|
||||
/**
|
||||
* Set error tolerances. This version specifies a scalar
|
||||
* relative tolerance, and a scalar absolute tolerance.
|
||||
*/
|
||||
virtual void setTolerances(doublereal reltol, doublereal abstol) {
|
||||
warn("setTolerances");
|
||||
}
|
||||
|
||||
/**
|
||||
* Specify a Jacobian evaluator. If this method is not called,
|
||||
* the Jacobian will be computed by finite difference.
|
||||
*/
|
||||
void setJacobian(Jacobian& jac) {
|
||||
warn("setJacobian");
|
||||
}
|
||||
|
||||
virtual void setLinearSolverType(int solverType) {
|
||||
warn("setLinearSolverType");
|
||||
}
|
||||
|
||||
virtual void setDenseLinearSolver() {
|
||||
warn("setDenseLinearSolver");
|
||||
}
|
||||
|
||||
virtual void setBandedLinearSolver(int m_upper, int m_lower) {
|
||||
warn("setBandedLinearSolver");
|
||||
}
|
||||
virtual void setMaxStepSize(doublereal dtmax) {
|
||||
warn("setMaxStepSize");
|
||||
}
|
||||
virtual void setMaxOrder(int n) {
|
||||
warn("setMaxOrder");
|
||||
}
|
||||
virtual void setMaxNumSteps(int n) {
|
||||
warn("setMaxNumSteps");
|
||||
}
|
||||
virtual void setInitialStepSize(doublereal h0) {
|
||||
warn("setInitialStepSize");
|
||||
}
|
||||
virtual void setStopTime(doublereal tstop) {
|
||||
warn("setStopTime");
|
||||
}
|
||||
virtual void setMaxErrTestFailures(int n) {
|
||||
warn("setMaxErrTestFailures");
|
||||
}
|
||||
virtual void setMaxNonlinIterations(int n) {
|
||||
warn("setMaxNonlinIterations");
|
||||
}
|
||||
virtual void setMaxNonlinConvFailures(int n) {
|
||||
warn("setMaxNonlinConvFailures");
|
||||
}
|
||||
virtual void inclAlgebraicInErrorTest(bool yesno) {
|
||||
warn("inclAlgebraicInErrorTest");
|
||||
}
|
||||
|
||||
//! Calculate consistent value of the starting solution given the starting
|
||||
//! solution derivatives
|
||||
/**
|
||||
* This method may be called if the initial conditions do not satisfy the
|
||||
* residual equation F = 0. Given the derivatives of all variables, this
|
||||
* method computes the initial y values.
|
||||
*/
|
||||
virtual void correctInitial_Y_given_Yp(doublereal* y, doublereal* yp,
|
||||
doublereal tout) {
|
||||
warn("correctInitial_Y_given_Yp");
|
||||
}
|
||||
|
||||
//! Calculate consistent value of the algebraic constraints and derivatives
|
||||
//! at the start of the problem
|
||||
/**
|
||||
* This method may be called if the initial conditions do not satisfy the
|
||||
* residual equation F = 0. Given the initial values of all differential
|
||||
* variables, it computes the initial values of all algebraic variables and
|
||||
* the initial derivatives of all differential variables.
|
||||
* @param y Calculated value of the solution vector after the procedure ends
|
||||
* @param yp Calculated value of the solution derivative after the procedure
|
||||
* @param tout The first value of t at which a soluton will be
|
||||
* requested (from IDASolve). (This is needed here to
|
||||
* determine the direction of integration and rough scale
|
||||
* in the independent variable t.
|
||||
*/
|
||||
virtual void correctInitial_YaYp_given_Yd(doublereal* y, doublereal* yp,
|
||||
doublereal tout) {
|
||||
warn("correctInitial_YaYp_given_Yd");
|
||||
}
|
||||
|
||||
/**
|
||||
* Solve the system of equations up to time tout.
|
||||
*/
|
||||
virtual int solve(doublereal tout) {
|
||||
warn("solve");
|
||||
return 0;
|
||||
}
|
||||
|
||||
/**
|
||||
* Take one internal step.
|
||||
*/
|
||||
virtual doublereal step(doublereal tout) {
|
||||
warn("step");
|
||||
return 0;
|
||||
}
|
||||
|
||||
//! Number of equations.
|
||||
int nEquations() const {
|
||||
return m_resid.nEquations();
|
||||
}
|
||||
|
||||
/**
|
||||
* initialize. Base class method does nothing.
|
||||
*/
|
||||
virtual void init(doublereal t0) {}
|
||||
|
||||
/**
|
||||
* Set a solver-specific input parameter.
|
||||
*/
|
||||
virtual void setInputParameter(int flag, doublereal value) {
|
||||
warn("setInputParameter");
|
||||
}
|
||||
|
||||
/**
|
||||
* Get the value of a solver-specific output parameter.
|
||||
*/
|
||||
virtual doublereal getOutputParameter(int flag) const {
|
||||
warn("getOutputParameter");
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
//! the current value of solution component k.
|
||||
virtual doublereal solution(int k) const {
|
||||
warn("solution");
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
virtual const doublereal* solutionVector() const {
|
||||
warn("solutionVector");
|
||||
return &m_dummy;
|
||||
}
|
||||
|
||||
//! the current value of the derivative of solution component k.
|
||||
virtual doublereal derivative(int k) const {
|
||||
warn("derivative");
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
virtual const doublereal* derivativeVector() const {
|
||||
warn("derivativeVector");
|
||||
return &m_dummy;
|
||||
}
|
||||
|
||||
protected:
|
||||
doublereal m_dummy;
|
||||
ResidJacEval& m_resid;
|
||||
|
||||
//! Number of total equations in the system
|
||||
integer m_neq;
|
||||
double m_time = 0.0;
|
||||
|
||||
private:
|
||||
void warn(const std::string& msg) const {
|
||||
writelog(">>>> Warning: method "+msg+" of base class "
|
||||
+"DAE_Solver called. Nothing done.\n");
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
//! Factor method for choosing a DAE solver
|
||||
/*!
|
||||
* @param itype String identifying the type (IDA is the only option)
|
||||
* @param f Residual function to be solved by the DAE algorithm
|
||||
* @returns a point to the instantiated DAE_Solver object
|
||||
*/
|
||||
DAE_Solver* newDAE_Solver(const std::string& itype, ResidJacEval& f);
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
@ -16,7 +16,7 @@ namespace Cantera
|
||||
{
|
||||
|
||||
/**
|
||||
* Virtual base class for ODE right-hand-side function evaluators.
|
||||
* Virtual base class for ODE/DAE right-hand-side function evaluators.
|
||||
* Classes derived from FuncEval evaluate the right-hand-side function
|
||||
* \f$ \vec{F}(t,\vec{y})\f$ in
|
||||
* \f[
|
||||
@ -31,13 +31,30 @@ public:
|
||||
virtual ~FuncEval() = default;
|
||||
|
||||
/**
|
||||
* Evaluate the right-hand-side function. Called by the integrator.
|
||||
* Evaluate the right-hand-side ODE function. Called by the integrator.
|
||||
* @param[in] t time.
|
||||
* @param[in] y solution vector, length neq()
|
||||
* @param[out] ydot rate of change of solution vector, length neq()
|
||||
* @param[in] p sensitivity parameter vector, length nparams()
|
||||
*/
|
||||
virtual void eval(double t, double* y, double* ydot, double* p)=0;
|
||||
virtual void eval(double t, double* y, double* ydot, double* p);
|
||||
|
||||
/**
|
||||
* Evaluate the right-hand-side DAE function. Called by the integrator.
|
||||
* @param[in] t time.
|
||||
* @param[in] y solution vector, length neq()
|
||||
* @param[out] ydot rate of change of solution vector, length neq()
|
||||
* @param[in] p sensitivity parameter vector, length nparams()
|
||||
* @param[out] residual the DAE residuals, length nparams()
|
||||
*/
|
||||
virtual void eval(double t, double* y, double* ydot, double* p,
|
||||
double* residual);
|
||||
|
||||
//! Given a vector of length neq(), mark which variables should be
|
||||
//! considered algebraic constraints
|
||||
virtual void getConstraints(double* constraints) {
|
||||
throw NotImplementedError("FuncEval::getConstraints");
|
||||
}
|
||||
|
||||
//! Evaluate the right-hand side using return code to indicate status.
|
||||
/*!
|
||||
@ -50,6 +67,17 @@ public:
|
||||
*/
|
||||
int eval_nothrow(double t, double* y, double* ydot);
|
||||
|
||||
//! Evaluate the right-hand side using return code to indicate status.
|
||||
/*!
|
||||
* Errors are indicated using the return value, rather than by throwing
|
||||
* exceptions. This method is used when calling from a C-based integrator
|
||||
* such as CVODES. Exceptions may either be stored or printed, based on the
|
||||
* setting of suppressErrors().
|
||||
* @returns 0 for a successful evaluation; 1 after a potentially-
|
||||
* recoverable error; -1 after an unrecoverable error.
|
||||
*/
|
||||
int eval_nothrow(double t, double* y, double* ydot, double* residaul);
|
||||
|
||||
/*! Evaluate the setup processes for the Jacobian preconditioner.
|
||||
* @param[in] t time.
|
||||
* @param[in] y solution vector, length neq()
|
||||
@ -102,6 +130,11 @@ public:
|
||||
throw NotImplementedError("FuncEval::getState");
|
||||
}
|
||||
|
||||
//! Fill in the vectors *y* and *ydot* with the current state of the system
|
||||
virtual void getState(double* y, double* ydot) {
|
||||
throw NotImplementedError("FuncEval::getState");
|
||||
}
|
||||
|
||||
//! Number of equations.
|
||||
virtual size_t neq()=0;
|
||||
|
||||
|
||||
133
include/cantera/numerics/IDAIntegrator.h
Normal file
133
include/cantera/numerics/IDAIntegrator.h
Normal file
@ -0,0 +1,133 @@
|
||||
/**
|
||||
* @file IDAIntegrator.h
|
||||
* Header file for class IDAIntegrator
|
||||
*/
|
||||
|
||||
// This file is part of Cantera. See License.txt in the top-level directory or
|
||||
// at http://www.cantera.org/license.txt for license and copyright information.
|
||||
|
||||
#ifndef CT_IDAIntegrator_H
|
||||
#define CT_IDAIntegrator_H
|
||||
|
||||
#include "cantera/numerics/Integrator.h"
|
||||
#include "cantera/base/ctexceptions.h"
|
||||
#include "sundials/sundials_nvector.h"
|
||||
|
||||
namespace Cantera
|
||||
{
|
||||
|
||||
|
||||
/**
|
||||
* Wrapper for Sundials IDA solver
|
||||
* @see FuncEval.h. Classes that use IDAIntegrator:
|
||||
* FlowReactor
|
||||
*/
|
||||
class IDAIntegrator : public Integrator
|
||||
{
|
||||
public:
|
||||
/**
|
||||
* Constructor. Default settings: dense Jacobian, no user-supplied
|
||||
* Jacobian function, Newton iteration.
|
||||
*/
|
||||
IDAIntegrator();
|
||||
virtual ~IDAIntegrator();
|
||||
virtual void setTolerances(double reltol, size_t n, double* abstol);
|
||||
virtual void setTolerances(double reltol, double abstol);
|
||||
virtual void setSensitivityTolerances(double reltol, double abstol);
|
||||
virtual void setProblemType(int probtype);
|
||||
virtual void initialize(double t0, FuncEval& func);
|
||||
virtual void reinitialize(double t0, FuncEval& func);
|
||||
virtual void integrate(double tout);
|
||||
virtual double step(double tout);
|
||||
virtual double& solution(size_t k);
|
||||
virtual double* solution();
|
||||
virtual int nEquations() const {
|
||||
return static_cast<int>(m_neq);
|
||||
}
|
||||
virtual int nEvals() const;
|
||||
virtual void setMaxOrder(int n) {
|
||||
m_maxord = n;
|
||||
}
|
||||
virtual void setMaxStepSize(double hmax);
|
||||
virtual void setMaxSteps(int nmax);
|
||||
virtual int maxSteps();
|
||||
virtual void setMaxErrTestFails(int n);
|
||||
virtual void setBandwidth(int N_Upper, int N_Lower) {
|
||||
m_mupper = N_Upper;
|
||||
m_mlower = N_Lower;
|
||||
}
|
||||
virtual int nSensParams() {
|
||||
return static_cast<int>(m_np);
|
||||
}
|
||||
virtual double sensitivity(size_t k, size_t p);
|
||||
|
||||
//! Returns a string listing the weighted error estimates associated
|
||||
//! with each solution component.
|
||||
//! This information can be used to identify which variables are
|
||||
//! responsible for integrator failures or unexpected small timesteps.
|
||||
virtual std::string getErrorInfo(int N);
|
||||
|
||||
//! Error message information provide by CVodes
|
||||
std::string m_error_message;
|
||||
|
||||
virtual void setMaxNonlinIterations(int n);
|
||||
virtual void setMaxNonlinConvFailures(int n);
|
||||
virtual void inclAlgebraicInErrorTest(bool yesno);
|
||||
virtual void setMethod(MethodType t);
|
||||
virtual void setIterator(IterType t);
|
||||
|
||||
protected:
|
||||
protected:
|
||||
//! Applies user-specified options to the underlying CVODES solver. Called
|
||||
//! during integrator initialization or reinitialization.
|
||||
void applyOptions();
|
||||
|
||||
private:
|
||||
void sensInit(double t0, FuncEval& func);
|
||||
size_t m_neq;
|
||||
void* m_ida_mem; //!< Pointer to the IDA memory for the problem
|
||||
void* m_linsol; //!< Sundials linear solver object
|
||||
void* m_linsol_matrix; //!< matrix used by Sundials
|
||||
FuncEval* m_func;
|
||||
double m_t0;
|
||||
double m_time; //!< The current integrator time
|
||||
N_Vector m_y, m_ydot, m_abstol;
|
||||
int m_type;
|
||||
int m_itol;
|
||||
int m_maxord;
|
||||
double m_reltol;
|
||||
double m_abstols;
|
||||
double m_reltolsens, m_abstolsens;
|
||||
size_t m_nabs;
|
||||
double m_hmax, m_hmin;
|
||||
int m_maxsteps;
|
||||
int m_maxErrTestFails;
|
||||
N_Vector* m_yS;
|
||||
N_Vector* m_ySdot;
|
||||
size_t m_np;
|
||||
int m_mupper, m_mlower;
|
||||
N_Vector m_constraints;
|
||||
|
||||
//! Indicates whether the sensitivities stored in m_yS have been updated
|
||||
//! for at the current integrator time.
|
||||
bool m_sens_ok;
|
||||
|
||||
//! Maximum number of nonlinear solver iterations at one solution
|
||||
/*!
|
||||
* If zero, this is the default of 4.
|
||||
*/
|
||||
int m_maxNonlinIters;
|
||||
|
||||
//! Maximum number of nonlinear convergence failures
|
||||
int m_maxNonlinConvFails;
|
||||
|
||||
//! If true, the algebraic variables don't contribute to error tolerances
|
||||
bool m_setSuppressAlg;
|
||||
|
||||
//! Initial IDA stepsize
|
||||
double m_init_step;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
@ -1,311 +0,0 @@
|
||||
/**
|
||||
* @file IDA_Solver.h
|
||||
* Header file for class IDA_Solver
|
||||
*/
|
||||
|
||||
// This file is part of Cantera. See License.txt in the top-level directory or
|
||||
// at https://cantera.org/license.txt for license and copyright information.
|
||||
|
||||
#ifndef CT_IDA_SOLVER_H
|
||||
#define CT_IDA_SOLVER_H
|
||||
|
||||
#include "DAE_Solver.h"
|
||||
#include "SundialsContext.h"
|
||||
|
||||
#include "sundials/sundials_nvector.h"
|
||||
|
||||
// These constants are defined internally in the IDA package, ida.c
|
||||
#define IDA_NN 0
|
||||
#define IDA_SS 1
|
||||
#define IDA_SV 2
|
||||
#define IDA_WF 3
|
||||
|
||||
#define REAL_WORKSPACE_SIZE 0
|
||||
|
||||
namespace Cantera
|
||||
{
|
||||
|
||||
class ResidData;
|
||||
|
||||
/**
|
||||
* Wrapper for Sundials IDA solver
|
||||
*
|
||||
* @attention This class currently does not have any test cases or examples. Its
|
||||
* implementation may be incomplete, and future changes to Cantera may
|
||||
* unexpectedly cause this class to stop working. If you use this class,
|
||||
* please consider contributing examples or test cases. In the absence of
|
||||
* new tests or examples, this class may be deprecated and removed in a
|
||||
* future version of Cantera. See
|
||||
* https://github.com/Cantera/cantera/issues/267 for additional information.
|
||||
*/
|
||||
class IDA_Solver : public DAE_Solver
|
||||
{
|
||||
public:
|
||||
//! Constructor.
|
||||
/*!
|
||||
* Default settings: dense Jacobian, no user-supplied Jacobian function,
|
||||
* Newton iteration.
|
||||
*
|
||||
* @param f Function that will supply the time dependent residual to be solved
|
||||
*/
|
||||
IDA_Solver(ResidJacEval& f);
|
||||
|
||||
virtual ~IDA_Solver();
|
||||
|
||||
virtual void setTolerances(doublereal reltol,
|
||||
doublereal* abstol);
|
||||
|
||||
virtual void setTolerances(doublereal reltol, doublereal abstol);
|
||||
|
||||
virtual void setLinearSolverType(int solverType);
|
||||
|
||||
//! Set up the problem to use a dense linear direct solver
|
||||
virtual void setDenseLinearSolver();
|
||||
|
||||
//! Set up the problem to use a band solver
|
||||
/*!
|
||||
* @param m_upper upper band width of the matrix
|
||||
* @param m_lower lower band width of the matrix
|
||||
*/
|
||||
virtual void setBandedLinearSolver(int m_upper, int m_lower);
|
||||
|
||||
virtual void setMaxOrder(int n);
|
||||
|
||||
//! Set the maximum number of time steps
|
||||
/*!
|
||||
* @param n input of maximum number of time steps
|
||||
*/
|
||||
virtual void setMaxNumSteps(int n);
|
||||
|
||||
//! Set the initial step size
|
||||
/*!
|
||||
* @param h0 initial step size value
|
||||
*/
|
||||
virtual void setInitialStepSize(doublereal h0);
|
||||
|
||||
//! Set the stop time
|
||||
/*!
|
||||
* @param tstop the independent variable value past which the solution is
|
||||
* not to proceed.
|
||||
*/
|
||||
virtual void setStopTime(doublereal tstop);
|
||||
|
||||
//! Get the current step size from IDA via a call
|
||||
/*!
|
||||
* @returns the current step size.
|
||||
*/
|
||||
virtual double getCurrentStepFromIDA();
|
||||
|
||||
//! Set the form of the Jacobian
|
||||
/*!
|
||||
* @param formJac Form of the Jacobian
|
||||
* 0 numerical Jacobian
|
||||
* 1 analytical Jacobian given by the evalJacobianDP() function
|
||||
*/
|
||||
virtual void setJacobianType(int formJac);
|
||||
|
||||
|
||||
virtual void setMaxErrTestFailures(int n);
|
||||
|
||||
//! Set the maximum number of nonlinear iterations on a timestep
|
||||
/*!
|
||||
* @param n Set the max iterations. The default is 4, which seems awfully low to me.
|
||||
*/
|
||||
virtual void setMaxNonlinIterations(int n);
|
||||
|
||||
//! Set the maximum number of nonlinear solver convergence failures
|
||||
/*!
|
||||
* @param n Value of nonlin failures. If value is exceeded, the calculation terminates.
|
||||
*/
|
||||
virtual void setMaxNonlinConvFailures(int n);
|
||||
|
||||
|
||||
virtual void inclAlgebraicInErrorTest(bool yesno);
|
||||
|
||||
/**
|
||||
* Get the value of a solver-specific output parameter.
|
||||
*/
|
||||
virtual doublereal getOutputParameter(int flag) const;
|
||||
|
||||
virtual void correctInitial_Y_given_Yp(doublereal* y, doublereal* yp,
|
||||
doublereal tout);
|
||||
|
||||
virtual void correctInitial_YaYp_given_Yd(doublereal* y, doublereal* yp, doublereal tout);
|
||||
|
||||
//! Step the system to a final value of the time
|
||||
/*!
|
||||
* @param tout Final value of the time
|
||||
* @returns the IDASolve() return flag
|
||||
*
|
||||
* The return values for IDASolve are described below. (The numerical return
|
||||
* values are defined above in this file.) All unsuccessful returns give a
|
||||
* negative return value.
|
||||
*
|
||||
* IDA_SUCCESS
|
||||
* IDASolve succeeded and no roots were found.
|
||||
*
|
||||
* IDA_ROOT_RETURN: IDASolve succeeded, and found one or more roots.
|
||||
* If nrtfn > 1, call IDAGetRootInfo to see which g_i were found
|
||||
* to have a root at (*tret).
|
||||
*
|
||||
* IDA_TSTOP_RETURN:
|
||||
* IDASolve returns computed results for the independent variable
|
||||
* value tstop. That is, tstop was reached.
|
||||
*
|
||||
* IDA_MEM_NULL:
|
||||
* The IDA_mem argument was NULL.
|
||||
*
|
||||
* IDA_ILL_INPUT:
|
||||
* One of the inputs to IDASolve is illegal. This includes the
|
||||
* situation when a component of the error weight vectors
|
||||
* becomes < 0 during internal stepping. It also includes the
|
||||
* situation where a root of one of the root functions was found
|
||||
* both at t0 and very near t0. The ILL_INPUT flag
|
||||
* will also be returned if the linear solver function IDA---
|
||||
* (called by the user after calling IDACreate) failed to set one
|
||||
* of the linear solver-related fields in ida_mem or if the linear
|
||||
* solver's init routine failed. In any case, the user should see
|
||||
* the printed error message for more details.
|
||||
*
|
||||
* IDA_TOO_MUCH_WORK:
|
||||
* The solver took mxstep internal steps but could not reach tout.
|
||||
* The default value for mxstep is MXSTEP_DEFAULT = 500.
|
||||
*
|
||||
* IDA_TOO_MUCH_ACC:
|
||||
* The solver could not satisfy the accuracy demanded by the user
|
||||
* for some internal step.
|
||||
*
|
||||
* IDA_ERR_FAIL:
|
||||
* Error test failures occurred too many times (=MXETF = 10) during
|
||||
* one internal step.
|
||||
*
|
||||
* IDA_CONV_FAIL:
|
||||
* Convergence test failures occurred too many times (= MXNCF = 10)
|
||||
* during one internal step.
|
||||
*
|
||||
* IDA_LSETUP_FAIL:
|
||||
* The linear solver's setup routine failed
|
||||
* in an unrecoverable manner.
|
||||
*
|
||||
* IDA_LSOLVE_FAIL:
|
||||
* The linear solver's solve routine failed
|
||||
* in an unrecoverable manner.
|
||||
*
|
||||
* IDA_CONSTR_FAIL:
|
||||
* The inequality constraints were violated,
|
||||
* and the solver was unable to recover.
|
||||
*
|
||||
* IDA_REP_RES_ERR:
|
||||
* The user's residual function repeatedly returned a recoverable
|
||||
* error flag, but the solver was unable to recover.
|
||||
*
|
||||
* IDA_RES_FAIL:
|
||||
* The user's residual function returned a nonrecoverable error
|
||||
* flag.
|
||||
*/
|
||||
virtual int solve(doublereal tout);
|
||||
|
||||
virtual doublereal step(doublereal tout);
|
||||
|
||||
virtual void init(doublereal t0);
|
||||
|
||||
virtual doublereal solution(int k) const;
|
||||
|
||||
virtual const doublereal* solutionVector() const;
|
||||
|
||||
virtual doublereal derivative(int k) const;
|
||||
|
||||
virtual const doublereal* derivativeVector() const;
|
||||
|
||||
void* IDAMemory() {
|
||||
return m_ida_mem;
|
||||
}
|
||||
|
||||
protected:
|
||||
//! Pointer to the IDA memory for the problem
|
||||
void* m_ida_mem = nullptr;
|
||||
void* m_linsol = nullptr; //!< Sundials linear solver object
|
||||
void* m_linsol_matrix = nullptr; //!< matrix used by Sundials
|
||||
SundialsContext m_sundials_ctx; //!< SUNContext object for Sundials>=6.0
|
||||
|
||||
//! Initial value of the time
|
||||
double m_t0 = 0.0;
|
||||
|
||||
//! Current value of the solution vector
|
||||
N_Vector m_y = nullptr;
|
||||
|
||||
//! Current value of the derivative of the solution vector
|
||||
N_Vector m_ydot = nullptr;
|
||||
N_Vector m_id = nullptr;
|
||||
N_Vector m_constraints = nullptr;
|
||||
N_Vector m_abstol = nullptr;
|
||||
int m_type = 0;
|
||||
|
||||
int m_itol = IDA_SS;
|
||||
int m_iter = 0;
|
||||
double m_reltol = 1e-9;
|
||||
double m_abstols = 1e-15;
|
||||
int m_nabs = 0;
|
||||
|
||||
//! Maximum value of the timestep allowed
|
||||
double m_hmax = 0.0;
|
||||
|
||||
//! Minimum value of the timestep allowed
|
||||
double m_hmin = 0.0;
|
||||
|
||||
//! Value of the initial time step
|
||||
double m_h0 = 0.0;
|
||||
|
||||
//! Maximum number of time steps allowed
|
||||
int m_maxsteps = 20000;
|
||||
|
||||
//! maximum time step order of the method
|
||||
int m_maxord = 0;
|
||||
|
||||
//! Form of the Jacobian
|
||||
/*!
|
||||
* 0 numerical Jacobian created by IDA
|
||||
* 1 analytical Jacobian. Must have populated the evalJacobianDP()
|
||||
* function in the ResidJacEval class.
|
||||
* 2 numerical Jacobian formed by the ResidJacEval class (unimplemented)
|
||||
*/
|
||||
int m_formJac = 0;
|
||||
|
||||
//! maximum time
|
||||
double m_tstop = 0.0;
|
||||
|
||||
//! Value of the previous, previous time
|
||||
double m_told_old = 0.0;
|
||||
|
||||
//! Value of the previous time
|
||||
double m_told = 0;
|
||||
|
||||
//! Value of the current time
|
||||
double m_tcurrent = 0;
|
||||
|
||||
//! Value of deltaT for the current step
|
||||
double m_deltat = 0.0;
|
||||
|
||||
//! maximum number of error test failures
|
||||
int m_maxErrTestFails = -1;
|
||||
|
||||
//! Maximum number of nonlinear solver iterations at one solution
|
||||
/*!
|
||||
* If zero, this is the default of 4.
|
||||
*/
|
||||
int m_maxNonlinIters = 0;
|
||||
|
||||
//! Maximum number of nonlinear convergence failures
|
||||
int m_maxNonlinConvFails = -1;
|
||||
|
||||
//! If true, the algebraic variables don't contribute to error tolerances
|
||||
int m_setSuppressAlg = 0;
|
||||
|
||||
std::unique_ptr<ResidData> m_fdata;
|
||||
int m_mupper = 0;
|
||||
int m_mlower = 0;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
@ -3,7 +3,7 @@
|
||||
*/
|
||||
|
||||
/**
|
||||
* @defgroup odeGroup ODE Integrators
|
||||
* @defgroup intGroup Zero-dimensional Integrators
|
||||
*/
|
||||
|
||||
// This file is part of Cantera. See License.txt in the top-level directory or
|
||||
@ -288,6 +288,16 @@ public:
|
||||
return stats;
|
||||
}
|
||||
|
||||
virtual void setMaxNonlinIterations(int n) {
|
||||
warn("setMaxNonlinIterations");
|
||||
}
|
||||
virtual void setMaxNonlinConvFailures(int n) {
|
||||
warn("setMaxNonlinConvFailures");
|
||||
}
|
||||
virtual void inclAlgebraicInErrorTest(bool yesno) {
|
||||
warn("inclAlgebraicInErrorTest");
|
||||
}
|
||||
|
||||
protected:
|
||||
//! Pointer to preconditioner object used in integration which is
|
||||
//! set by setPreconditioner and initialized inside of
|
||||
@ -295,6 +305,7 @@ protected:
|
||||
shared_ptr<PreconditionerBase> m_preconditioner;
|
||||
//! Type of preconditioning used in applyOptions
|
||||
PreconditionerSide m_prec_side = PreconditionerSide::NO_PRECONDITION;
|
||||
// methods for DAE solvers
|
||||
|
||||
private:
|
||||
doublereal m_dummy;
|
||||
@ -304,7 +315,7 @@ private:
|
||||
}
|
||||
};
|
||||
|
||||
// defined in ODE_integrators.cpp
|
||||
// defined in Integrators.cpp
|
||||
Integrator* newIntegrator(const std::string& itype);
|
||||
|
||||
} // namespace
|
||||
|
||||
50
include/cantera/numerics/sundials_headers.h
Normal file
50
include/cantera/numerics/sundials_headers.h
Normal file
@ -0,0 +1,50 @@
|
||||
#ifndef CT_SUNDIALS_HEADERS
|
||||
#define CT_SUNDIALS_HEADERS
|
||||
|
||||
#include "sundials/sundials_types.h"
|
||||
#include "sundials/sundials_math.h"
|
||||
#include "sundials/sundials_nvector.h"
|
||||
#include "nvector/nvector_serial.h"
|
||||
#include "cvodes/cvodes.h"
|
||||
#include "idas/idas.h"
|
||||
#if CT_SUNDIALS_VERSION >= 30
|
||||
#if CT_SUNDIALS_USE_LAPACK
|
||||
#include "sunlinsol/sunlinsol_lapackdense.h"
|
||||
#include "sunlinsol/sunlinsol_lapackband.h"
|
||||
#else
|
||||
#include "sunlinsol/sunlinsol_dense.h"
|
||||
#include "sunlinsol/sunlinsol_band.h"
|
||||
#endif
|
||||
#include "sunlinsol/sunlinsol_spgmr.h"
|
||||
#include "cvodes/cvodes_direct.h"
|
||||
#include "cvodes/cvodes_diag.h"
|
||||
#include "cvodes/cvodes_spils.h"
|
||||
#include "idas/idas_direct.h"
|
||||
#include "idas/idas_spils.h"
|
||||
#else
|
||||
#if CT_SUNDIALS_USE_LAPACK
|
||||
#include "cvodes/cvodes_lapack.h"
|
||||
#include "idas/idas_lapack.h"
|
||||
#else
|
||||
#include "cvodes/cvodes_dense.h"
|
||||
#include "cvodes/cvodes_band.h"
|
||||
#include "idas/idas_dense.h"
|
||||
#include "idas/idas_band.h"
|
||||
#endif
|
||||
#include "cvodes/cvodes_diag.h"
|
||||
#include "cvodes/cvodes_spgmr.h"
|
||||
#include "idas/idas_spgmr.h"
|
||||
#endif
|
||||
|
||||
#define CV_SS 1
|
||||
#define IDA_SS 1
|
||||
#define CV_SV 2
|
||||
#define IDA_SV 2
|
||||
|
||||
#if CT_SUNDIALS_VERSION < 25
|
||||
typedef int sd_size_t;
|
||||
#else
|
||||
typedef long int sd_size_t;
|
||||
#endif
|
||||
|
||||
#endif
|
||||
@ -6,55 +6,167 @@
|
||||
#ifndef CT_FLOWREACTOR_H
|
||||
#define CT_FLOWREACTOR_H
|
||||
|
||||
#include "Reactor.h"
|
||||
#include "IdealGasReactor.h"
|
||||
|
||||
namespace Cantera
|
||||
{
|
||||
|
||||
//! Adiabatic flow in a constant-area duct.
|
||||
class FlowReactor : public Reactor
|
||||
class FlowReactor : public IdealGasReactor
|
||||
{
|
||||
public:
|
||||
FlowReactor() = default;
|
||||
//! Note: currently assumes a cylinder
|
||||
FlowReactor(double area=1.0,
|
||||
double sa_to_vol=-1,
|
||||
double ss_atol=1e-14,
|
||||
double ss_rtol=1e-7,
|
||||
int max_ss_steps = 20000,
|
||||
int max_ss_error_fails=10);
|
||||
|
||||
virtual std::string type() const {
|
||||
return "FlowReactor";
|
||||
}
|
||||
|
||||
virtual void getState(doublereal* y);
|
||||
virtual void getState(double* y)
|
||||
{
|
||||
throw CanteraError("FlowReactor::getState", "Not Implemented!");
|
||||
}
|
||||
virtual void getState(double* y, double* ydot);
|
||||
|
||||
virtual void initialize(doublereal t0 = 0.0);
|
||||
virtual void eval(double t, double* LHS, double* RHS);
|
||||
virtual void updateState(doublereal* y);
|
||||
void setMassFlowRate(double mdot);
|
||||
/*!
|
||||
* Evaluate the reactor governing equations. Called by ReactorNet::eval.
|
||||
* @param[in] t time.
|
||||
* @param[in] y solution vector, length neq()
|
||||
* @param[out] ydot rate of change of solution vector, length neq()
|
||||
* @param[in] params sensitivity parameter vector, length ReactorNet::nparams()
|
||||
*/
|
||||
virtual void eval(double t, double* y,
|
||||
double* ydot, double* params)
|
||||
{
|
||||
throw CanteraError("FlowReactor::eval", "Not Implemented!");
|
||||
}
|
||||
|
||||
/*!
|
||||
* Evaluate the reactor governing equations. Called by ReactorNet::eval.
|
||||
* @param[in] t time.
|
||||
* @param[in] y solution vector, length neq()
|
||||
* @param[in] ydot rate of change of solution vector, length neq()
|
||||
* @param[in] params sensitivity parameter vector, length ReactorNet::nparams()
|
||||
* @param[out] residual resisduals vector, length neq()
|
||||
*/
|
||||
virtual void eval(double t, double* y,
|
||||
double* ydot, double* params,
|
||||
double* residual);
|
||||
|
||||
//! Given a vector of length neq(), mark which variables should be
|
||||
//! considered algebraic constraints
|
||||
virtual void getConstraints(double* constraints) {
|
||||
// mark all variables differential equations unless otherwise specified
|
||||
std::fill(constraints, constraints + m_nv, 1.0);
|
||||
// the species coverages are algebraic constraints
|
||||
std::fill(constraints + m_non_spec_eq + m_nsp, constraints + m_nv, 0.0);
|
||||
}
|
||||
|
||||
|
||||
virtual void syncState();
|
||||
virtual void updateState(doublereal* y);
|
||||
|
||||
void setMassFlowRate(double mdot);
|
||||
|
||||
void setTimeConstant(doublereal tau) {
|
||||
m_fctr = 1.0/tau;
|
||||
void setMassFlowRate(doublereal mdot) {
|
||||
m_rho = m_thermo->density();
|
||||
m_u = mdot/(m_rho * m_area);
|
||||
}
|
||||
|
||||
//! The current gas speed in the reactor [m/s]
|
||||
double speed() const {
|
||||
return m_speed;
|
||||
return m_u;
|
||||
}
|
||||
double distance() const {
|
||||
return m_dist;
|
||||
|
||||
//! The area of the reactor [m^2]
|
||||
double area() const {
|
||||
return m_area;
|
||||
}
|
||||
|
||||
//! Sets the area of the reactor [m^2]
|
||||
void setArea(double area) {
|
||||
m_area = area;
|
||||
setMassFlowRate(m_rho * m_u * area);
|
||||
}
|
||||
|
||||
//! The ratio of the reactor's surface area to volume ratio [m^-1]
|
||||
//! @note If the surface area to volume ratio is unspecified by the user,
|
||||
//! this will be calculated assuming the reactor is a cylinder.
|
||||
double surfaceAreaToVolumeRatio() const {
|
||||
if (m_sa_to_vol > 0)
|
||||
return m_sa_to_vol;
|
||||
|
||||
// assuming a cylinder, volume = Pi * r^2 * L, and perimeter = 2 * Pi * r * L
|
||||
// where L is the length, and r is the radius of the reactor
|
||||
// hence, perimeter / area = 2 * Pi * r * L / (Pi * L * r^2) = 2 / r
|
||||
return 2.0 / sqrt(m_area / Pi);
|
||||
}
|
||||
|
||||
//! Set the reactor's surface area to volume ratio [m^-1]
|
||||
void setSurfaceAreaToVolumeRatio(double sa_to_vol) {
|
||||
m_sa_to_vol = sa_to_vol;
|
||||
}
|
||||
|
||||
//! Set the steady state tolerances used to determine the initial state for
|
||||
//! surface coverages
|
||||
void setSteadyStateAtol(double atol) {
|
||||
m_ss_atol = atol;
|
||||
}
|
||||
|
||||
//! Set the steady state tolerances used to determine the initial state for
|
||||
//! surface coverages
|
||||
void setSteadyStateRtol(double rtol) {
|
||||
m_ss_rtol = rtol;
|
||||
}
|
||||
|
||||
//! Set the steady state tolerances used to determine the initial state for
|
||||
//! surface coverages
|
||||
void setSteadyStateMaxSteps(int max_steps) {
|
||||
m_max_ss_steps = max_steps;
|
||||
}
|
||||
|
||||
//! Set the steady state tolerances used to determine the initial state for
|
||||
//! surface coverages
|
||||
void setSteadyStateMaxErrorFailures(int max_fails) {
|
||||
m_max_ss_error_fails = max_fails;
|
||||
}
|
||||
|
||||
|
||||
//! Return the index in the solution vector for this reactor of the
|
||||
//! component named *nm*. Possible values for *nm* are "X" (position),
|
||||
//! "U", the name of a homogeneous phase species, or the name of a surface
|
||||
//! species.
|
||||
virtual size_t componentIndex(const std::string& nm) const;
|
||||
|
||||
virtual void updateSurfaceState(double* y);
|
||||
|
||||
protected:
|
||||
double m_speed = 0.0;
|
||||
double m_dist = 0.0;
|
||||
double m_T = 0.0;
|
||||
double m_fctr = 1.0e10;
|
||||
double m_rho0 = 0.0;
|
||||
double m_speed0 = 0.0;
|
||||
double m_P0 = 0.0;
|
||||
double m_h0 = 0.0;
|
||||
double m_u, m_T, m_P, m_rho;
|
||||
//! offset to the species equations
|
||||
const size_t m_non_spec_eq = 4;
|
||||
//! reactor area [m^2]
|
||||
double m_area;
|
||||
//! reactor surface area to volume ratio [m^-1]
|
||||
double m_sa_to_vol;
|
||||
//! temporary storage for surface species production rates
|
||||
vector_fp m_sdot_temp;
|
||||
//! temporary storage for species partial molar enthalipes
|
||||
vector_fp m_hk;
|
||||
|
||||
//! steady-state relative tolerance, used to determine initial surface coverages
|
||||
double m_ss_rtol;
|
||||
//! steady-state absolute tolerance, used to determine initial surface coverages
|
||||
double m_ss_atol;
|
||||
//! maximum number of steady-state coverage integrator-steps
|
||||
int m_max_ss_steps;
|
||||
//! maximum number of steady-state integrator error test failures
|
||||
int m_max_ss_error_fails;
|
||||
};
|
||||
}
|
||||
|
||||
|
||||
@ -99,6 +99,17 @@ public:
|
||||
*/
|
||||
virtual void getState(doublereal* y);
|
||||
|
||||
//! Get the current state and derivative vector of the reactor for a DAE solver
|
||||
/*!
|
||||
* @param[out] y state vector representing the initial state of the reactor
|
||||
* @param[out] ydot state vector representing the initial derivatives of the
|
||||
* reactor
|
||||
*/
|
||||
virtual void getState(doublereal* y, doublereal* ydot)
|
||||
{
|
||||
throw CanteraError("Reactor::getState", "Not Implemented");
|
||||
}
|
||||
|
||||
virtual void initialize(doublereal t0 = 0.0);
|
||||
|
||||
//! Evaluate the reactor governing equations. Called by ReactorNet::eval.
|
||||
@ -109,6 +120,28 @@ public:
|
||||
//! coefficients for governing equations, length m_nv, default values 0
|
||||
virtual void eval(double t, double* LHS, double* RHS);
|
||||
|
||||
/*!
|
||||
* Evaluate the reactor governing equations. Called by ReactorNet::eval.
|
||||
* @param[in] t time.
|
||||
* @param[in] y solution vector, length neq()
|
||||
* @param[in] ydot rate of change of solution vector, length neq()
|
||||
* @param[in] params sensitivity parameter vector, length ReactorNet::nparams()
|
||||
* @param[out] residual resisduals vector, length neq()
|
||||
*/
|
||||
virtual void evalEqs(doublereal t, doublereal* y,
|
||||
doublereal* ydot, doublereal* params,
|
||||
doublereal* residual)
|
||||
{
|
||||
throw CanteraError("Reactor::evalEqs", "Not Implemented");
|
||||
}
|
||||
|
||||
//! Given a vector of length neq(), mark which variables should be
|
||||
//! considered algebraic constraints
|
||||
virtual void getConstraints(double* constraints) {
|
||||
throw CanteraError("Reactor::getConstraints", "Not Implemented");
|
||||
}
|
||||
|
||||
|
||||
virtual void syncState();
|
||||
|
||||
//! Set the state of the reactor to correspond to the state vector *y*.
|
||||
|
||||
@ -138,7 +138,7 @@ public:
|
||||
//! Return a reference to the *n*-th Wall connected to this reactor.
|
||||
WallBase& wall(size_t n);
|
||||
|
||||
void addSurface(ReactorSurface* surf);
|
||||
virtual void addSurface(ReactorSurface* surf);
|
||||
|
||||
//! Return a reference to the *n*-th ReactorSurface connected to this
|
||||
//! reactor
|
||||
|
||||
@ -68,7 +68,11 @@ public:
|
||||
//! sensitivity equations.
|
||||
void setSensitivityTolerances(double rtol, double atol);
|
||||
|
||||
//! @}
|
||||
//! Set the maximum number of internal integration time-steps the
|
||||
//! integrator will take before reaching the next output time
|
||||
//! @param nmax The maximum number of steps, setting this value
|
||||
//! to zero disables this option.
|
||||
virtual void setMaxSteps(int nmax);
|
||||
|
||||
//! Current value of the simulation time.
|
||||
doublereal time() {
|
||||
@ -98,6 +102,13 @@ public:
|
||||
//! Problem type of integrator
|
||||
std::string linearSolverType() const;
|
||||
|
||||
//! Returns the maximum number of internal integration time-steps the
|
||||
//! integrator will take before reaching the next output time
|
||||
//!
|
||||
virtual int maxSteps() {
|
||||
return m_nmax;
|
||||
}
|
||||
|
||||
/**
|
||||
* Advance the state of all reactors in time. Take as many internal
|
||||
* timesteps as necessary to reach *time*.
|
||||
@ -200,10 +211,18 @@ public:
|
||||
virtual void eval(doublereal t, doublereal* y,
|
||||
doublereal* ydot, doublereal* p);
|
||||
|
||||
//! eval coupling for IDA / DAEs
|
||||
virtual void eval(doublereal t, doublereal* y,
|
||||
doublereal* ydot, doublereal* p,
|
||||
doublereal* residaul);
|
||||
|
||||
virtual void getState(doublereal* y);
|
||||
|
||||
//! Return k-th derivative at the current time
|
||||
virtual void getDerivative(int k, double* dky);
|
||||
virtual void getState(doublereal* y, doublereal* ydot);
|
||||
|
||||
virtual void getConstraints(doublereal* constraints);
|
||||
|
||||
virtual size_t nparams() {
|
||||
return m_sens_params.size();
|
||||
@ -322,6 +341,12 @@ protected:
|
||||
int m_maxErrTestFails = 0;
|
||||
bool m_verbose = false;
|
||||
|
||||
// flag indicating whether this is a flow reactor net or a regular reactor net
|
||||
bool m_is_dae;
|
||||
|
||||
//! The maximum number of steps allowed in the integrator
|
||||
size_t m_nmax;
|
||||
|
||||
//! Names corresponding to each sensitivity parameter
|
||||
std::vector<std::string> m_paramNames;
|
||||
|
||||
|
||||
@ -55,6 +55,14 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
|
||||
void setMassFlowRate(double) except +translate_exception
|
||||
double speed()
|
||||
double distance()
|
||||
void setArea(double) except +translate_exception
|
||||
double area() except +translate_exception
|
||||
void setSurfaceAreaToVolumeRatio(double) except +translate_exception
|
||||
double surfaceAreaToVolumeRatio() except +translate_exception
|
||||
void setSteadyStateAtol(double) except +translate_exception
|
||||
void setSteadyStateRtol(double) except +translate_exception
|
||||
void setSteadyStateMaxSteps(int) except +translate_exception
|
||||
void setSteadyStateMaxErrorFailures(int) except +translate_exception
|
||||
|
||||
# walls
|
||||
cdef cppclass CxxWallBase "Cantera::WallBase":
|
||||
|
||||
@ -445,21 +445,75 @@ cdef class FlowReactor(Reactor):
|
||||
"""
|
||||
reactor_type = "FlowReactor"
|
||||
|
||||
def __init__(self, *args, **kwargs):
|
||||
area = kwargs.pop('area', None)
|
||||
sa_to_vol = kwargs.pop('surface_area_to_volume_ratio', None)
|
||||
ss_atol = kwargs.pop('steady_state_atol', None)
|
||||
ss_rtol = kwargs.pop('steady_state_rtol', None)
|
||||
ss_max_steps = kwargs.pop('steady_state_max_steps', None)
|
||||
ss_max_fail = kwargs.pop('steady_state_max_error_failures', None)
|
||||
super(FlowReactor, self).__init__(*args, **kwargs)
|
||||
if area is not None:
|
||||
self.area = area
|
||||
if sa_to_vol is not None:
|
||||
self.surface_area_to_volume_ratio = sa_to_vol
|
||||
if ss_atol is not None:
|
||||
self.steady_state_atol = ss_atol
|
||||
if ss_rtol is not None:
|
||||
self.steady_state_rtol = ss_rtol
|
||||
if ss_max_steps is not None:
|
||||
self.steady_state_max_steps = ss_max_steps
|
||||
if ss_max_fail is not None:
|
||||
self.steady_state_max_error_failures = ss_max_fail
|
||||
|
||||
property mass_flow_rate:
|
||||
""" Mass flow rate per unit area [kg/m^2*s] """
|
||||
def __set__(self, double value):
|
||||
(<CxxFlowReactor*>self.reactor).setMassFlowRate(value)
|
||||
|
||||
property area:
|
||||
""" Get/set the area of the reactor [m^2] """
|
||||
def __get__(self):
|
||||
return (<CxxFlowReactor*>self.reactor).area()
|
||||
def __set__(self, area):
|
||||
(<CxxFlowReactor*>self.reactor).setArea(area)
|
||||
|
||||
property steady_state_atol:
|
||||
""" Set the steady-state tolerances used to determine the initial surface
|
||||
species coverages"""
|
||||
def __set__(self, atol):
|
||||
(<CxxFlowReactor*>self.reactor).setSteadyStateAtol(atol)
|
||||
|
||||
property steady_state_rtol:
|
||||
""" Set the steady-state tolerances used to determine the initial surface
|
||||
species coverages"""
|
||||
def __set__(self, rtol):
|
||||
(<CxxFlowReactor*>self.reactor).setSteadyStateRtol(rtol)
|
||||
|
||||
property steady_state_max_steps:
|
||||
""" Set the maximum number of integrator steps used to determine the initial surface
|
||||
species coverages"""
|
||||
def __set__(self, nsteps):
|
||||
(<CxxFlowReactor*>self.reactor).setSteadyStateMaxSteps(nsteps)
|
||||
|
||||
property steady_state_max_error_failures:
|
||||
""" Set the maximum number of integrator error failures allowed when determining
|
||||
the initial surface species coverages"""
|
||||
def __set__(self, nsteps):
|
||||
(<CxxFlowReactor*>self.reactor).setSteadyStateMaxErrorFailures(nsteps)
|
||||
|
||||
property surface_area_to_volume_ratio:
|
||||
""" Get/set the surface area to volume ratio of the reactor [m^-1] """
|
||||
def __get__(self):
|
||||
return (<CxxFlowReactor*>self.reactor).surfaceAreaToVolumeRatio()
|
||||
def __set__(self, sa_to_vol):
|
||||
(<CxxFlowReactor*>self.reactor).setSurfaceAreaToVolumeRatio(sa_to_vol)
|
||||
|
||||
property speed:
|
||||
""" Speed [m/s] of the flow in the reactor at the current position """
|
||||
def __get__(self):
|
||||
return (<CxxFlowReactor*>self.reactor).speed()
|
||||
|
||||
property distance:
|
||||
""" The distance of the fluid element from the inlet of the reactor."""
|
||||
def __get__(self):
|
||||
return (<CxxFlowReactor*>self.reactor).distance()
|
||||
|
||||
|
||||
cdef class ExtensibleReactor(Reactor):
|
||||
"""
|
||||
|
||||
115
samples/python/surface_chemistry/1D_pfr_surfchem.py
Normal file
115
samples/python/surface_chemistry/1D_pfr_surfchem.py
Normal file
@ -0,0 +1,115 @@
|
||||
"""
|
||||
Solution of a 1-D steady state plug-flow reactor with surface chemistry.
|
||||
Assumes:
|
||||
- adiabatic
|
||||
- frictionless
|
||||
- constant area, cylindrical reactor
|
||||
|
||||
Based off the jupyter notebook created by Yuanjie Jiang
|
||||
"""
|
||||
|
||||
|
||||
import numpy as np
|
||||
import cantera as ct
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
|
||||
mech = 'SiF4_NH3_mec.yaml'
|
||||
# import the models for gas and bulk
|
||||
gas, bulk_si, bulk_n = ct.import_phases(mech, ['gas', 'SiBulk', 'NBulk'])
|
||||
# import the model for gas-Si-N interface
|
||||
gas_si_n_interface = ct.Interface(mech, 'SI3N4', [gas, bulk_si, bulk_n])
|
||||
|
||||
# Set the initial conditions
|
||||
T0 = 1713 # K
|
||||
p0 = 2 * ct.one_atm / 760.0 # Pa ~2Torr
|
||||
gas.TPX = T0, p0, "NH3:6, SiF4:1"
|
||||
bulk_si.TP = T0, p0
|
||||
bulk_n.TP = T0, p0
|
||||
gas_si_n_interface.TP = T0, p0
|
||||
D = 5.08e-2 # diameter of the tube [m]
|
||||
Ac = np.pi * D**2 / 4 # cross section of the tube [m]
|
||||
mu = 5.7e-5 # kg/(m-s) dynamic viscosity
|
||||
perim = np.pi * D # perimeter of the tube
|
||||
# calculate the site fractions of surface species at the entrance of the tube
|
||||
# at steady state
|
||||
u0 = 11.53 # m/s initial velocity of the flow
|
||||
|
||||
|
||||
reactor = ct.FlowReactor(gas, area=Ac)
|
||||
# set the mass flow-rate
|
||||
reactor.mass_flow_rate = gas.density * u0 * Ac
|
||||
|
||||
rsurf = ct.ReactorSurface(gas_si_n_interface, reactor)
|
||||
net = ct.ReactorNet([reactor])
|
||||
soln = ct.SolutionArray(gas, extra=['t', 'speed', 'coverages'])
|
||||
|
||||
while net.time < 0.7:
|
||||
net.step()
|
||||
soln.append(T=reactor.thermo.T,
|
||||
P=reactor.thermo.P,
|
||||
X=reactor.thermo.X,
|
||||
t=net.time,
|
||||
speed=reactor.speed,
|
||||
coverages=rsurf.coverages)
|
||||
|
||||
f, ax = plt.subplots(4, 2, figsize=(9, 16), dpi=96)
|
||||
# plot gas velocity along the flow direction
|
||||
ax[0, 0].plot(soln.t, soln.speed, color='C0')
|
||||
ax[0, 0].set_xlabel('Distance (m)')
|
||||
ax[0, 0].set_ylabel('Velocity (m/s)')
|
||||
|
||||
# plot gas density along the flow direction
|
||||
ax[0, 1].plot(soln.t, soln.density, color='C1')
|
||||
ax[0, 1].set_xlabel('Distance (m)')
|
||||
ax[0, 1].set_ylabel(r'Density ($\mathregular{kg/m^3}$)')
|
||||
ax[0, 1].ticklabel_format(axis='y', style='sci',
|
||||
scilimits=(-2, 2)) # scientific notation
|
||||
|
||||
# plot major and minor gas species separately
|
||||
minor_idx = []
|
||||
major_idx = []
|
||||
for i, name in enumerate(gas.species_names):
|
||||
mean = np.mean(soln(name).Y)
|
||||
if mean <= 0.01:
|
||||
minor_idx.append(i)
|
||||
else:
|
||||
major_idx.append(i)
|
||||
|
||||
# plot minor gas species along the flow direction
|
||||
for i in minor_idx:
|
||||
style = '-' if i < 10 else '--'
|
||||
ax[1, 0].plot(soln.t, soln(gas.species_names[i]).Y,
|
||||
label=gas.species_names[i], linestyle=style)
|
||||
ax[1, 0].legend(fontsize=7.5, loc='best')
|
||||
ax[1, 0].set_xlabel('Distance (m)')
|
||||
ax[1, 0].set_ylabel('Mass Fraction')
|
||||
ax[1, 0].ticklabel_format(axis='y', style='sci',
|
||||
scilimits=(-2, 2)) # scientific notation
|
||||
|
||||
# plot major gas species along the flow direction
|
||||
for j in major_idx:
|
||||
ax[1, 1].plot(soln.t, soln(gas.species_names[j]).Y,
|
||||
label=gas.species_names[j])
|
||||
ax[1, 1].legend(fontsize=8, loc='best')
|
||||
ax[1, 1].set_xlabel('Distance (m)')
|
||||
ax[1, 1].set_ylabel('Mass Fraction')
|
||||
|
||||
# plot the pressure of the gas along the flow direction
|
||||
ax[2, 0].plot(soln.t, soln.P, color='C2')
|
||||
ax[2, 0].set_xlabel('Distance (m)')
|
||||
ax[2, 0].set_ylabel('Pressure (Pa)')
|
||||
|
||||
# plot the site fraction of the surface species along the flow direction
|
||||
for i, name in enumerate(gas_si_n_interface.species_names):
|
||||
ax[2, 1].plot(soln.t, soln.coverages[:, i], label=name)
|
||||
ax[2, 1].legend(fontsize=8)
|
||||
ax[2, 1].set_xlabel('Distance (m)')
|
||||
ax[2, 1].set_ylabel('Site Fraction')
|
||||
|
||||
# plot the temperature profile along the flow direction
|
||||
ax[3, 0].plot(soln.t, soln.T[:], color='C3')
|
||||
ax[3, 0].set_xlabel('Distance (m)')
|
||||
ax[3, 0].set_ylabel('Temperature (K)')
|
||||
f.tight_layout(pad=0.5)
|
||||
plt.show()
|
||||
@ -9,32 +9,7 @@
|
||||
#include <iostream>
|
||||
using namespace std;
|
||||
|
||||
#include "sundials/sundials_types.h"
|
||||
#include "sundials/sundials_math.h"
|
||||
#include "sundials/sundials_nvector.h"
|
||||
#include "nvector/nvector_serial.h"
|
||||
#include "cvodes/cvodes.h"
|
||||
#if CT_SUNDIALS_USE_LAPACK
|
||||
#include "sunlinsol/sunlinsol_lapackdense.h"
|
||||
#include "sunlinsol/sunlinsol_lapackband.h"
|
||||
#else
|
||||
#include "sunlinsol/sunlinsol_dense.h"
|
||||
#include "sunlinsol/sunlinsol_band.h"
|
||||
#endif
|
||||
#include "sunlinsol/sunlinsol_spgmr.h"
|
||||
#include "cvodes/cvodes_direct.h"
|
||||
#include "cvodes/cvodes_diag.h"
|
||||
#include "cvodes/cvodes_spils.h"
|
||||
|
||||
#define CV_SS 1
|
||||
#define CV_SV 2
|
||||
|
||||
#ifdef SUNDIALS_INDEX_TYPE
|
||||
// For SUNDIALS >= 3.2.0
|
||||
typedef SUNDIALS_INDEX_TYPE sd_size_t;
|
||||
#else
|
||||
typedef long int sd_size_t;
|
||||
#endif
|
||||
#include "cantera/numerics/sundials_headers.h"
|
||||
|
||||
namespace {
|
||||
|
||||
|
||||
@ -1,27 +0,0 @@
|
||||
//! @file DAE_solvers.cpp Factory routine for picking the DAE solver package
|
||||
|
||||
// This file is part of Cantera. See License.txt in the top-level directory or
|
||||
// at https://cantera.org/license.txt for license and copyright information.
|
||||
|
||||
#include "cantera/base/ct_defs.h"
|
||||
#include "cantera/numerics/DAE_Solver.h"
|
||||
#include "cantera/numerics/IDA_Solver.h"
|
||||
|
||||
// DAE_DEVEL is turned off at the current time
|
||||
#define DAE_DEVEL
|
||||
#ifdef DAE_DEVEL
|
||||
|
||||
namespace Cantera
|
||||
{
|
||||
DAE_Solver* newDAE_Solver(const std::string& itype, ResidJacEval& f)
|
||||
{
|
||||
if (itype == "IDA") {
|
||||
return new IDA_Solver(f);
|
||||
} else {
|
||||
throw CanteraError("newDAE_Solver",
|
||||
"unknown DAE solver: "+itype);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#endif
|
||||
@ -37,6 +37,39 @@ int FuncEval::eval_nothrow(double t, double* y, double* ydot)
|
||||
return 0; // successful evaluation
|
||||
}
|
||||
|
||||
int FuncEval::eval_nothrow(double t, double* y, double* ydot, double* r)
|
||||
{
|
||||
try {
|
||||
eval(t, y, ydot, m_sens_params.data(), r);
|
||||
} catch (CanteraError& err) {
|
||||
if (suppressErrors()) {
|
||||
m_errors.push_back(err.what());
|
||||
} else {
|
||||
writelog(err.what());
|
||||
}
|
||||
return 1; // possibly recoverable error
|
||||
} catch (std::exception& err) {
|
||||
if (suppressErrors()) {
|
||||
m_errors.push_back(err.what());
|
||||
} else {
|
||||
writelog("FuncEval::eval_nothrow: unhandled exception:\n");
|
||||
writelog(err.what());
|
||||
writelogendl();
|
||||
}
|
||||
return -1; // unrecoverable error
|
||||
} catch (...) {
|
||||
std::string msg = "FuncEval::eval_nothrow: unhandled exception"
|
||||
" of unknown type\n";
|
||||
if (suppressErrors()) {
|
||||
m_errors.push_back(msg);
|
||||
} else {
|
||||
writelog(msg);
|
||||
}
|
||||
return -1; // unrecoverable error
|
||||
}
|
||||
return 0; // successful evaluation
|
||||
}
|
||||
|
||||
std::string FuncEval::getErrors() const {
|
||||
std::stringstream errs;
|
||||
for (const auto& err : m_errors) {
|
||||
@ -112,4 +145,14 @@ int FuncEval::preconditioner_solve_nothrow(double* rhs, double* output)
|
||||
return 0; // successful evaluation
|
||||
}
|
||||
|
||||
void FuncEval::eval(double t, double* y, double* ydot, double* p)
|
||||
{
|
||||
throw CanteraError("FuncEval::eval", "Not implemented!");
|
||||
}
|
||||
|
||||
void FuncEval::eval(double t, double* y, double* ydot, double* p, double* residual)
|
||||
{
|
||||
throw CanteraError("FuncEval::eval", "Not implemented!");
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
588
src/numerics/IDAIntegrator.cpp
Normal file
588
src/numerics/IDAIntegrator.cpp
Normal file
@ -0,0 +1,588 @@
|
||||
//! @file IDAIntegrator.cpp
|
||||
|
||||
// This file is part of Cantera. See License.txt in the top-level directory or
|
||||
// at http://www.cantera.org/license.txt for license and copyright information.
|
||||
|
||||
#include "cantera/numerics/IDAIntegrator.h"
|
||||
#include "cantera/base/stringUtils.h"
|
||||
|
||||
#include "cantera/numerics/sundials_headers.h"
|
||||
|
||||
using namespace std;
|
||||
namespace Cantera
|
||||
{
|
||||
|
||||
extern "C" {
|
||||
//! Function called by IDA to evaluate the residual, given y and ydot.
|
||||
/*!
|
||||
* IDA allows passing in a void* pointer to access external data. Instead of
|
||||
* requiring the user to provide a residual function directly to IDA (which
|
||||
* would require using the sundials data types N_Vector, etc.), we define
|
||||
* this function as the single function that IDA always calls. The real
|
||||
* evaluation of the residual is done by an instance of a subclass of
|
||||
* ResidEval, passed in to this function as a pointer in the parameters.
|
||||
*
|
||||
* FROM IDA WRITEUP -> What the IDA solver expects as a return flag from its
|
||||
* residual routines:
|
||||
*
|
||||
* A IDAResFn res should return a value of 0 if successful, a positive value
|
||||
* if a recoverable error occured (e.g. yy has an illegal value), or a
|
||||
* negative value if a nonrecoverable error occured. In the latter case, the
|
||||
* program halts. If a recoverable error occured, the integrator will
|
||||
* attempt to correct and retry.
|
||||
*/
|
||||
static int ida_rhs(realtype t, N_Vector y, N_Vector ydot, N_Vector r, void* f_data)
|
||||
{
|
||||
FuncEval* f = (FuncEval*) f_data;
|
||||
return f->eval_nothrow(t, NV_DATA_S(y), NV_DATA_S(ydot), NV_DATA_S(r));
|
||||
}
|
||||
|
||||
//! Function called by IDA when an error is encountered instead of
|
||||
//! writing to stdout. Here, save the error message provided by IDA so
|
||||
//! that it can be included in the subsequently raised CanteraError.
|
||||
static void ida_err(int error_code, const char* module,
|
||||
const char* function, char* msg, void* eh_data)
|
||||
{
|
||||
IDAIntegrator* integrator = (IDAIntegrator*) eh_data;
|
||||
integrator->m_error_message = msg;
|
||||
integrator->m_error_message += "\n";
|
||||
}
|
||||
}
|
||||
|
||||
IDAIntegrator::IDAIntegrator() :
|
||||
m_ida_mem(0),
|
||||
m_linsol(0),
|
||||
m_linsol_matrix(0),
|
||||
m_func(0),
|
||||
m_t0(0.0),
|
||||
m_y(0),
|
||||
m_ydot(0),
|
||||
m_abstol(0),
|
||||
m_type(DENSE+NOJAC),
|
||||
m_itol(IDA_SS),
|
||||
m_maxord(0),
|
||||
m_reltol(1.e-9),
|
||||
m_abstols(1.e-15),
|
||||
m_nabs(0),
|
||||
m_hmax(0.0),
|
||||
m_hmin(0.0),
|
||||
m_maxsteps(20000),
|
||||
m_maxErrTestFails(-1),
|
||||
m_mupper(0),
|
||||
m_mlower(0),
|
||||
m_constraints(0),
|
||||
m_maxNonlinIters(0),
|
||||
m_maxNonlinConvFails(-1),
|
||||
m_setSuppressAlg(0),
|
||||
m_init_step(1.e-14)
|
||||
{
|
||||
}
|
||||
|
||||
IDAIntegrator::~IDAIntegrator()
|
||||
{
|
||||
if (m_ida_mem) {
|
||||
IDAFree(&m_ida_mem);
|
||||
}
|
||||
if (m_y) {
|
||||
N_VDestroy_Serial(m_y);
|
||||
}
|
||||
if (m_ydot) {
|
||||
N_VDestroy_Serial(m_ydot);
|
||||
}
|
||||
if (m_abstol) {
|
||||
N_VDestroy_Serial(m_abstol);
|
||||
}
|
||||
if (m_constraints) {
|
||||
N_VDestroy_Serial(m_constraints);
|
||||
}
|
||||
}
|
||||
|
||||
double& IDAIntegrator::solution(size_t k)
|
||||
{
|
||||
return NV_Ith_S(m_y,k);
|
||||
}
|
||||
|
||||
double* IDAIntegrator::solution()
|
||||
{
|
||||
return NV_DATA_S(m_y);
|
||||
}
|
||||
|
||||
void IDAIntegrator::setTolerances(double reltol, size_t n, double* abstol)
|
||||
{
|
||||
m_itol = IDA_SV;
|
||||
m_nabs = n;
|
||||
if (n != m_neq) {
|
||||
if (m_abstol) {
|
||||
N_VDestroy_Serial(m_abstol);
|
||||
}
|
||||
m_abstol = N_VNew_Serial(static_cast<sd_size_t>(n));
|
||||
}
|
||||
for (size_t i=0; i<n; i++) {
|
||||
NV_Ith_S(m_abstol, i) = abstol[i];
|
||||
}
|
||||
m_reltol = reltol;
|
||||
}
|
||||
|
||||
void IDAIntegrator::setTolerances(double reltol, double abstol)
|
||||
{
|
||||
m_itol = IDA_SS;
|
||||
m_reltol = reltol;
|
||||
m_abstols = abstol;
|
||||
}
|
||||
|
||||
void IDAIntegrator::setSensitivityTolerances(double reltol, double abstol)
|
||||
{
|
||||
m_reltolsens = reltol;
|
||||
m_abstolsens = abstol;
|
||||
}
|
||||
|
||||
|
||||
void IDAIntegrator::setProblemType(int probtype)
|
||||
{
|
||||
m_type = probtype;
|
||||
}
|
||||
|
||||
void IDAIntegrator::setMaxStepSize(double hmax)
|
||||
{
|
||||
m_hmax = hmax;
|
||||
if (m_ida_mem) {
|
||||
IDASetMaxStep(m_ida_mem, hmax);
|
||||
}
|
||||
}
|
||||
|
||||
void IDAIntegrator::setMaxSteps(int nmax)
|
||||
{
|
||||
m_maxsteps = nmax;
|
||||
if (m_ida_mem) {
|
||||
IDASetMaxNumSteps(m_ida_mem, m_maxsteps);
|
||||
}
|
||||
}
|
||||
|
||||
int IDAIntegrator::maxSteps()
|
||||
{
|
||||
return m_maxsteps;
|
||||
}
|
||||
|
||||
void IDAIntegrator::setMaxErrTestFails(int n)
|
||||
{
|
||||
m_maxErrTestFails = n;
|
||||
if (m_ida_mem) {
|
||||
IDASetMaxErrTestFails(m_ida_mem, n);
|
||||
}
|
||||
}
|
||||
|
||||
void IDAIntegrator::setMaxNonlinIterations(int n)
|
||||
{
|
||||
m_maxNonlinIters = n;
|
||||
if (m_ida_mem)
|
||||
{
|
||||
int flag = IDASetMaxNonlinIters(m_ida_mem, m_maxNonlinIters);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDAIntegrator::init",
|
||||
"IDASetmaxNonlinIters failed.");
|
||||
}
|
||||
}
|
||||
|
||||
if (m_setSuppressAlg != 0) {
|
||||
int flag = IDASetSuppressAlg(m_ida_mem, m_setSuppressAlg);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDAIntegrator::init", "IDASetSuppressAlg failed.");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void IDAIntegrator::setMaxNonlinConvFailures(int n)
|
||||
{
|
||||
m_maxNonlinIters = n;
|
||||
if (m_ida_mem)
|
||||
{
|
||||
int flag = IDASetMaxConvFails(m_ida_mem, m_maxNonlinConvFails);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDAIntegrator::init",
|
||||
"IDASetMaxConvFails failed.");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void IDAIntegrator::inclAlgebraicInErrorTest(bool yesno)
|
||||
{
|
||||
if (yesno) {
|
||||
m_setSuppressAlg = 0;
|
||||
} else {
|
||||
m_setSuppressAlg = 1;
|
||||
}
|
||||
|
||||
if (m_ida_mem)
|
||||
{
|
||||
int flag = IDASetSuppressAlg(m_ida_mem, m_setSuppressAlg);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDAIntegrator::init", "IDASetSuppressAlg failed.");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void IDAIntegrator::initialize(double t0, FuncEval& func)
|
||||
{
|
||||
m_neq = func.neq();
|
||||
m_t0 = t0;
|
||||
m_time = t0;
|
||||
m_func = &func;
|
||||
func.clearErrors();
|
||||
|
||||
if (m_y) {
|
||||
N_VDestroy_Serial(m_y); // free solution vector if already allocated
|
||||
}
|
||||
m_y = N_VNew_Serial(static_cast<sd_size_t>(m_neq)); // allocate solution vector
|
||||
N_VConst(0.0, m_y);
|
||||
|
||||
if (m_ydot)
|
||||
{
|
||||
N_VDestroy_Serial(m_ydot); // free derivative vector if already allocated
|
||||
}
|
||||
m_ydot = N_VNew_Serial(m_neq);
|
||||
N_VConst(0.0, m_ydot);
|
||||
|
||||
// check abs tolerance array size
|
||||
if (m_itol == IDA_SV && m_nabs < m_neq) {
|
||||
throw CanteraError("IDAIntegrator::initialize",
|
||||
"not enough absolute tolerance values specified.");
|
||||
}
|
||||
|
||||
if (m_constraints) {
|
||||
N_VDestroy_Serial(m_constraints);
|
||||
}
|
||||
m_constraints = N_VNew_Serial(static_cast<sd_size_t>(m_neq));
|
||||
// set the constraints
|
||||
func.getConstraints(NV_DATA_S(m_constraints));
|
||||
|
||||
|
||||
// get the initial conditions
|
||||
func.getState(NV_DATA_S(m_y), NV_DATA_S(m_ydot));
|
||||
|
||||
if (m_ida_mem) {
|
||||
IDAFree(&m_ida_mem);
|
||||
}
|
||||
|
||||
//! Create the IDA solver
|
||||
m_ida_mem = IDACreate();
|
||||
if (!m_ida_mem) {
|
||||
throw CanteraError("IDAIntegrator::initialize",
|
||||
"IDACreate failed.");
|
||||
}
|
||||
|
||||
int flag = IDAInit(m_ida_mem, ida_rhs, m_t0, m_y, m_ydot);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
if (flag == IDA_MEM_FAIL) {
|
||||
throw CanteraError("IDAIntegrator::initialize",
|
||||
"Memory allocation failed.");
|
||||
} else if (flag == IDA_ILL_INPUT) {
|
||||
throw CanteraError("IDAIntegrator::initialize",
|
||||
"Illegal value for IDAInit input argument.");
|
||||
} else {
|
||||
throw CanteraError("IDAIntegrator::initialize",
|
||||
"IDAInit failed.");
|
||||
}
|
||||
}
|
||||
|
||||
// set constraints
|
||||
flag = IDASetId(m_ida_mem, m_constraints);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
if (flag == IDA_MEM_NULL) {
|
||||
throw CanteraError("IDAIntegrator::initialize",
|
||||
"Memory allocation failed.");
|
||||
} else {
|
||||
throw CanteraError("IDAIntegrator::initialize",
|
||||
"IDASetId failed.");
|
||||
}
|
||||
}
|
||||
|
||||
flag = IDASetErrHandlerFn(m_ida_mem, &ida_err, this);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDAIntegrator::initialize",
|
||||
"IDASetErrHandlerFn failed.");
|
||||
}
|
||||
|
||||
if (m_itol == IDA_SV) {
|
||||
flag = IDASVtolerances(m_ida_mem, m_reltol, m_abstol);
|
||||
} else {
|
||||
flag = IDASStolerances(m_ida_mem, m_reltol, m_abstols);
|
||||
}
|
||||
if (flag != IDA_SUCCESS) {
|
||||
if (flag == IDA_MEM_FAIL) {
|
||||
throw CanteraError("IDAIntegrator::initialize",
|
||||
"Memory allocation failed.");
|
||||
} else if (flag == IDA_ILL_INPUT) {
|
||||
throw CanteraError("IDAIntegrator::initialize",
|
||||
"Illegal value for IDAInit input argument.");
|
||||
} else {
|
||||
throw CanteraError("IDAIntegrator::initialize",
|
||||
"IDAInit failed.");
|
||||
}
|
||||
}
|
||||
|
||||
flag = IDASetUserData(m_ida_mem, &func);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDAIntegrator::initialize",
|
||||
"IDASetUserData failed.");
|
||||
}
|
||||
applyOptions();
|
||||
}
|
||||
|
||||
void IDAIntegrator::reinitialize(double t0, FuncEval& func)
|
||||
{
|
||||
m_t0 = t0;
|
||||
m_time = t0;
|
||||
func.getState(NV_DATA_S(m_y));
|
||||
m_func = &func;
|
||||
func.clearErrors();
|
||||
|
||||
int result = IDAReInit(m_ida_mem, m_t0, m_y, m_ydot);
|
||||
if (result != IDA_SUCCESS) {
|
||||
throw CanteraError("IDAIntegrator::reinitialize",
|
||||
"IDAReInit failed. result = {}", result);
|
||||
}
|
||||
applyOptions();
|
||||
}
|
||||
|
||||
void IDAIntegrator::applyOptions()
|
||||
{
|
||||
if (m_type == DENSE + NOJAC) {
|
||||
sd_size_t N = static_cast<sd_size_t>(m_neq);
|
||||
#if CT_SUNDIALS_VERSION >= 30
|
||||
SUNLinSolFree((SUNLinearSolver) m_linsol);
|
||||
SUNMatDestroy((SUNMatrix) m_linsol_matrix);
|
||||
m_linsol_matrix = SUNDenseMatrix(N, N);
|
||||
#if CT_SUNDIALS_USE_LAPACK
|
||||
m_linsol = SUNLapackDense(m_y, (SUNMatrix) m_linsol_matrix);
|
||||
#else
|
||||
m_linsol = SUNDenseLinearSolver(m_y, (SUNMatrix) m_linsol_matrix);
|
||||
#endif
|
||||
IDADlsSetLinearSolver(m_ida_mem, (SUNLinearSolver) m_linsol,
|
||||
(SUNMatrix) m_linsol_matrix);
|
||||
#else
|
||||
#if CT_SUNDIALS_USE_LAPACK
|
||||
IDALapackDense(m_ida_mem, N);
|
||||
#else
|
||||
IDADense(m_ida_mem, N);
|
||||
#endif
|
||||
#endif
|
||||
} else if (m_type == DIAG) {
|
||||
throw CanteraError("IDAIntegrator::applyOptions",
|
||||
"Cannot use a diagonal matrix with IDA.");
|
||||
} else if (m_type == GMRES) {
|
||||
#if CT_SUNDIALS_VERSION >= 30
|
||||
m_linsol = SUNSPGMR(m_y, PREC_NONE, 0);
|
||||
IDASpilsSetLinearSolver(m_ida_mem, (SUNLinearSolver) m_linsol);
|
||||
#else
|
||||
IDASpgmr(m_ida_mem, PREC_NONE, 0);
|
||||
#endif
|
||||
} else if (m_type == BAND + NOJAC) {
|
||||
sd_size_t N = static_cast<sd_size_t>(m_neq);
|
||||
long int nu = m_mupper;
|
||||
long int nl = m_mlower;
|
||||
#if CT_SUNDIALS_VERSION >= 30
|
||||
SUNLinSolFree((SUNLinearSolver) m_linsol);
|
||||
SUNMatDestroy((SUNMatrix) m_linsol_matrix);
|
||||
m_linsol_matrix = SUNBandMatrix(N, nu, nl, nu+nl);
|
||||
#if CT_SUNDIALS_USE_LAPACK
|
||||
m_linsol = SUNLapackBand(m_y, (SUNMatrix) m_linsol_matrix);
|
||||
#else
|
||||
m_linsol = SUNBandLinearSolver(m_y, (SUNMatrix) m_linsol_matrix);
|
||||
#endif
|
||||
IDADlsSetLinearSolver(m_ida_mem, (SUNLinearSolver) m_linsol,
|
||||
(SUNMatrix) m_linsol_matrix);
|
||||
#else
|
||||
#if CT_SUNDIALS_USE_LAPACK
|
||||
IDALapackBand(m_ida_mem, N, nu, nl);
|
||||
#else
|
||||
IDABand(m_ida_mem, N, nu, nl);
|
||||
#endif
|
||||
#endif
|
||||
} else {
|
||||
throw CanteraError("IDAIntegrator::applyOptions",
|
||||
"unsupported option");
|
||||
}
|
||||
|
||||
if (m_init_step > 0) {
|
||||
IDASetInitStep(m_ida_mem, m_init_step);
|
||||
}
|
||||
|
||||
if (m_maxord > 0) {
|
||||
IDASetMaxOrd(m_ida_mem, m_maxord);
|
||||
}
|
||||
if (m_maxsteps > 0) {
|
||||
IDASetMaxNumSteps(m_ida_mem, m_maxsteps);
|
||||
}
|
||||
if (m_hmax > 0) {
|
||||
IDASetMaxStep(m_ida_mem, m_hmax);
|
||||
}
|
||||
if (m_maxNonlinIters > 0) {
|
||||
int flag = IDASetMaxNonlinIters(m_ida_mem, m_maxNonlinIters);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDAIntegrator::init",
|
||||
"IDASetmaxNonlinIters failed.");
|
||||
}
|
||||
}
|
||||
if (m_maxNonlinConvFails > 0) {
|
||||
int flag = IDASetMaxConvFails(m_ida_mem, m_maxNonlinConvFails);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDAIntegrator::init",
|
||||
"IDASetMaxConvFails failed.");
|
||||
}
|
||||
}
|
||||
if (m_setSuppressAlg != 0) {
|
||||
int flag = IDASetSuppressAlg(m_ida_mem, m_setSuppressAlg);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDAIntegrator::init", "IDASetSuppressAlg failed.");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void IDAIntegrator::sensInit(double t0, FuncEval& func)
|
||||
{
|
||||
m_np = func.nparams();
|
||||
m_sens_ok = false;
|
||||
|
||||
N_Vector y = N_VNew_Serial(static_cast<sd_size_t>(func.neq()));
|
||||
m_yS = N_VCloneVectorArray_Serial(static_cast<sd_size_t>(m_np), y);
|
||||
for (size_t n = 0; n < m_np; n++) {
|
||||
N_VConst(0.0, m_yS[n]);
|
||||
}
|
||||
N_VDestroy_Serial(y);
|
||||
N_Vector ydot = N_VNew_Serial(static_cast<sd_size_t>(func.neq()));
|
||||
m_ySdot = N_VCloneVectorArray_Serial(static_cast<sd_size_t>(m_np), ydot);
|
||||
for (size_t n = 0; n < m_np; n++) {
|
||||
N_VConst(0.0, m_ySdot[n]);
|
||||
}
|
||||
|
||||
int flag = IDASensInit(m_ida_mem, static_cast<sd_size_t>(m_np),
|
||||
IDA_STAGGERED, IDASensResFn(0), m_yS, m_ySdot);
|
||||
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDAIntegrator::sensInit", "Error in IDASensInit");
|
||||
}
|
||||
vector_fp atol(m_np);
|
||||
for (size_t n = 0; n < m_np; n++) {
|
||||
// This scaling factor is tuned so that reaction and species enthalpy
|
||||
// sensitivities can be computed simultaneously with the same abstol.
|
||||
atol[n] = m_abstolsens / func.m_paramScales[n];
|
||||
}
|
||||
flag = IDASensSStolerances(m_ida_mem, m_reltolsens, atol.data());
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDAIntegrator::sensInit", "Error in IDASensSStolerances");
|
||||
}
|
||||
}
|
||||
|
||||
void IDAIntegrator::integrate(double tout)
|
||||
{
|
||||
if (tout == m_time) {
|
||||
return;
|
||||
}
|
||||
int flag = IDASolve(m_ida_mem, tout, &m_time, m_y, m_ydot, IDA_NORMAL);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
string f_errs = m_func->getErrors();
|
||||
if (!f_errs.empty()) {
|
||||
f_errs = "Exceptions caught during RHS evaluation:\n" + f_errs;
|
||||
}
|
||||
throw CanteraError("IDAIntegrator::integrate",
|
||||
"IDA error encountered. Error code: {}\n{}\n"
|
||||
"{}"
|
||||
"Components with largest weighted error estimates:\n{}",
|
||||
flag, m_error_message, f_errs, getErrorInfo(10));
|
||||
}
|
||||
}
|
||||
|
||||
double IDAIntegrator::step(double tout)
|
||||
{
|
||||
int flag = IDASolve(m_ida_mem, tout, &m_time, m_y, m_ydot, IDA_ONE_STEP);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
string f_errs = m_func->getErrors();
|
||||
if (!f_errs.empty()) {
|
||||
f_errs = "Exceptions caught during RHS evaluation:\n" + f_errs;
|
||||
}
|
||||
throw CanteraError("IDAIntegrator::step",
|
||||
"IDA error encountered. Error code: {}\n{}\n"
|
||||
"{}"
|
||||
"Components with largest weighted error estimates:\n{}",
|
||||
flag, f_errs, m_error_message, getErrorInfo(10));
|
||||
|
||||
}
|
||||
return m_time;
|
||||
}
|
||||
|
||||
int IDAIntegrator::nEvals() const
|
||||
{
|
||||
// neither IDAGetNumRhsEvals, IDADlsGetNumRhsEvals, IDASpilsGetNumRhsEvals,
|
||||
// or IDAGetNumLinResEvals seem to exist
|
||||
throw CanteraError("IDAIntegrator::nEvals",
|
||||
"Not implemented.");
|
||||
}
|
||||
|
||||
double IDAIntegrator::sensitivity(size_t k, size_t p)
|
||||
{
|
||||
if (m_time == m_t0) {
|
||||
// calls to IDAGetSens are only allowed after a successful time step.
|
||||
return 0.0;
|
||||
}
|
||||
if (!m_sens_ok && m_np) {
|
||||
int flag = IDAGetSens(m_ida_mem, &m_time, m_yS);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDAIntegrator::sensitivity",
|
||||
"IDAGetSens failed. Error code: {}", flag);
|
||||
}
|
||||
m_sens_ok = true;
|
||||
}
|
||||
|
||||
if (k >= m_neq) {
|
||||
throw CanteraError("IDAIntegrator::sensitivity",
|
||||
"sensitivity: k out of range ({})", k);
|
||||
}
|
||||
if (p >= m_np) {
|
||||
throw CanteraError("IDAIntegrator::sensitivity",
|
||||
"sensitivity: p out of range ({})", p);
|
||||
}
|
||||
return NV_Ith_S(m_yS[p],k);
|
||||
}
|
||||
|
||||
string IDAIntegrator::getErrorInfo(int N)
|
||||
{
|
||||
N_Vector errs = N_VNew_Serial(static_cast<sd_size_t>(m_neq));
|
||||
N_Vector errw = N_VNew_Serial(static_cast<sd_size_t>(m_neq));
|
||||
IDAGetErrWeights(m_ida_mem, errw);
|
||||
IDAGetEstLocalErrors(m_ida_mem, errs);
|
||||
|
||||
vector<tuple<double, double, size_t> > weightedErrors;
|
||||
for (size_t i=0; i<m_neq; i++) {
|
||||
double err = NV_Ith_S(errs, i) * NV_Ith_S(errw, i);
|
||||
weightedErrors.emplace_back(-abs(err), err, i);
|
||||
}
|
||||
N_VDestroy(errs);
|
||||
N_VDestroy(errw);
|
||||
|
||||
N = std::min(N, static_cast<int>(m_neq));
|
||||
sort(weightedErrors.begin(), weightedErrors.end());
|
||||
fmt::memory_buffer s;
|
||||
for (int i=0; i<N; i++) {
|
||||
fmt_append(s, "{}: {}\n",
|
||||
get<2>(weightedErrors[i]), get<1>(weightedErrors[i]));
|
||||
}
|
||||
return to_string(s);
|
||||
}
|
||||
|
||||
void IDAIntegrator::setMethod(MethodType t)
|
||||
{
|
||||
if (t != BDF_Method) {
|
||||
// IDA only has the BDF method
|
||||
throw CanteraError("IDAIntegrator::setMethod", "unknown method");
|
||||
}
|
||||
}
|
||||
|
||||
void IDAIntegrator::setIterator(IterType t)
|
||||
{
|
||||
if (t != Newton_Iter) {
|
||||
// IDA only has the newton iterator
|
||||
throw CanteraError("IDAIntegrator::setIterator", "unknown iterator");
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
@ -1,656 +0,0 @@
|
||||
//! @file IDA_Solver.cpp
|
||||
|
||||
// This file is part of Cantera. See License.txt in the top-level directory or
|
||||
// at https://cantera.org/license.txt for license and copyright information.
|
||||
|
||||
#include "cantera/numerics/IDA_Solver.h"
|
||||
#include "cantera/base/stringUtils.h"
|
||||
|
||||
#include "sundials/sundials_types.h"
|
||||
#include "sundials/sundials_math.h"
|
||||
#include "ida/ida.h"
|
||||
|
||||
#if CT_SUNDIALS_USE_LAPACK
|
||||
#include "sunlinsol/sunlinsol_lapackdense.h"
|
||||
#include "sunlinsol/sunlinsol_lapackband.h"
|
||||
#else
|
||||
#include "sunlinsol/sunlinsol_dense.h"
|
||||
#include "sunlinsol/sunlinsol_band.h"
|
||||
#endif
|
||||
|
||||
#include "sunlinsol/sunlinsol_spgmr.h"
|
||||
#include "ida/ida_direct.h"
|
||||
#include "ida/ida_spils.h"
|
||||
#include "nvector/nvector_serial.h"
|
||||
|
||||
typedef long int sd_size_t;
|
||||
|
||||
namespace Cantera
|
||||
{
|
||||
|
||||
//! A simple class to hold an array of parameter values and a pointer to an
|
||||
//! instance of a subclass of ResidEval.
|
||||
class ResidData
|
||||
{
|
||||
public:
|
||||
ResidData(ResidJacEval* f, IDA_Solver* s, int npar = 0) {
|
||||
m_func = f;
|
||||
m_solver = s;
|
||||
}
|
||||
|
||||
virtual ~ResidData() {
|
||||
}
|
||||
|
||||
ResidJacEval* m_func;
|
||||
IDA_Solver* m_solver;
|
||||
};
|
||||
}
|
||||
|
||||
extern "C" {
|
||||
//! Function called by IDA to evaluate the residual, given y and ydot.
|
||||
/*!
|
||||
* IDA allows passing in a void* pointer to access external data. Instead of
|
||||
* requiring the user to provide a residual function directly to IDA (which
|
||||
* would require using the sundials data types N_Vector, etc.), we define
|
||||
* this function as the single function that IDA always calls. The real
|
||||
* evaluation of the residual is done by an instance of a subclass of
|
||||
* ResidEval, passed in to this function as a pointer in the parameters.
|
||||
*
|
||||
* FROM IDA WRITEUP -> What the IDA solver expects as a return flag from its
|
||||
* residual routines:
|
||||
*
|
||||
* A IDAResFn res should return a value of 0 if successful, a positive value
|
||||
* if a recoverable error occurred (for example, yy has an illegal value), or a
|
||||
* negative value if a nonrecoverable error occurred. In the latter case, the
|
||||
* program halts. If a recoverable error occurred, the integrator will
|
||||
* attempt to correct and retry.
|
||||
*/
|
||||
static int ida_resid(realtype t, N_Vector y, N_Vector ydot, N_Vector r, void* f_data)
|
||||
{
|
||||
Cantera::ResidData* d = (Cantera::ResidData*) f_data;
|
||||
Cantera::ResidJacEval* f = d->m_func;
|
||||
Cantera::IDA_Solver* s = d->m_solver;
|
||||
double delta_t = s->getCurrentStepFromIDA();
|
||||
// TODO evaluate evalType. Assumed to be Base_ResidEval
|
||||
int flag = f->evalResidNJ(t, delta_t, NV_DATA_S(y), NV_DATA_S(ydot),
|
||||
NV_DATA_S(r));
|
||||
if (flag < 0) {
|
||||
// This signals to IDA that a nonrecoverable error has occurred.
|
||||
return flag;
|
||||
} else {
|
||||
return 0;
|
||||
}
|
||||
}
|
||||
|
||||
//! Function called by by IDA to evaluate the Jacobian, given y and ydot.
|
||||
/*!
|
||||
* typedef int (*IDADlsDenseJacFn)(sd_size_t N, realtype t, realtype c_j,
|
||||
* N_Vector y, N_Vector yp, N_Vector r,
|
||||
* DlsMat Jac, void *user_data,
|
||||
* N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
|
||||
*
|
||||
* A IDADlsDenseJacFn should return
|
||||
* - 0 if successful,
|
||||
* - a positive int if a recoverable error occurred, or
|
||||
* - a negative int if a nonrecoverable error occurred.
|
||||
*
|
||||
* In the case of a recoverable error return, the integrator will attempt to
|
||||
* recover by reducing the stepsize (which changes cj).
|
||||
*/
|
||||
static int ida_jacobian(realtype t, realtype c_j, N_Vector y, N_Vector yp,
|
||||
N_Vector r, SUNMatrix Jac, void *f_data,
|
||||
N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
|
||||
{
|
||||
Cantera::ResidData* d = (Cantera::ResidData*) f_data;
|
||||
Cantera::ResidJacEval* f = d->m_func;
|
||||
Cantera::IDA_Solver* s = d->m_solver;
|
||||
double delta_t = s->getCurrentStepFromIDA();
|
||||
double** cols;
|
||||
if (SUNMatGetID(Jac) == SUNMATRIX_DENSE) {
|
||||
cols = SM_COLS_D(Jac);
|
||||
} else if (SUNMatGetID(Jac) == SUNMATRIX_BAND) {
|
||||
cols = SM_COLS_B(Jac);
|
||||
} else {
|
||||
return 1; // Unknown SUNMatrix type
|
||||
}
|
||||
f->evalJacobianDP(t, delta_t, c_j, NV_DATA_S(y), NV_DATA_S(yp),
|
||||
cols, NV_DATA_S(r));
|
||||
return 0;
|
||||
}
|
||||
}
|
||||
|
||||
namespace Cantera
|
||||
{
|
||||
|
||||
IDA_Solver::IDA_Solver(ResidJacEval& f) :
|
||||
DAE_Solver(f)
|
||||
{
|
||||
}
|
||||
|
||||
IDA_Solver::~IDA_Solver()
|
||||
{
|
||||
if (m_ida_mem) {
|
||||
IDAFree(&m_ida_mem);
|
||||
}
|
||||
if (m_y) {
|
||||
N_VDestroy_Serial(m_y);
|
||||
}
|
||||
if (m_ydot) {
|
||||
N_VDestroy_Serial(m_ydot);
|
||||
}
|
||||
if (m_abstol) {
|
||||
N_VDestroy_Serial(m_abstol);
|
||||
}
|
||||
if (m_constraints) {
|
||||
N_VDestroy_Serial(m_constraints);
|
||||
}
|
||||
}
|
||||
|
||||
doublereal IDA_Solver::solution(int k) const
|
||||
{
|
||||
return NV_Ith_S(m_y,k);
|
||||
}
|
||||
|
||||
const doublereal* IDA_Solver::solutionVector() const
|
||||
{
|
||||
return NV_DATA_S(m_y);
|
||||
}
|
||||
|
||||
doublereal IDA_Solver::derivative(int k) const
|
||||
{
|
||||
return NV_Ith_S(m_ydot,k);
|
||||
}
|
||||
|
||||
const doublereal* IDA_Solver::derivativeVector() const
|
||||
{
|
||||
return NV_DATA_S(m_ydot);
|
||||
}
|
||||
|
||||
void IDA_Solver::setTolerances(double reltol, double* abstol)
|
||||
{
|
||||
m_itol = IDA_SV;
|
||||
if (!m_abstol) {
|
||||
#if CT_SUNDIALS_VERSION >= 60
|
||||
m_abstol = N_VNew_Serial(m_neq, m_sundials_ctx.get());
|
||||
#else
|
||||
m_abstol = N_VNew_Serial(m_neq);
|
||||
#endif
|
||||
}
|
||||
for (int i = 0; i < m_neq; i++) {
|
||||
NV_Ith_S(m_abstol, i) = abstol[i];
|
||||
}
|
||||
m_reltol = reltol;
|
||||
if (m_ida_mem) {
|
||||
int flag = IDASVtolerances(m_ida_mem, m_reltol, m_abstol);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::setTolerances",
|
||||
"Memory allocation failed.");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void IDA_Solver::setTolerances(doublereal reltol, doublereal abstol)
|
||||
{
|
||||
m_itol = IDA_SS;
|
||||
m_reltol = reltol;
|
||||
m_abstols = abstol;
|
||||
if (m_ida_mem) {
|
||||
int flag = IDASStolerances(m_ida_mem, m_reltol, m_abstols);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::setTolerances",
|
||||
"Memory allocation failed.");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void IDA_Solver::setLinearSolverType(int solverType)
|
||||
{
|
||||
m_type = solverType;
|
||||
}
|
||||
|
||||
void IDA_Solver::setDenseLinearSolver()
|
||||
{
|
||||
setLinearSolverType(0);
|
||||
}
|
||||
|
||||
void IDA_Solver::setBandedLinearSolver(int m_upper, int m_lower)
|
||||
{
|
||||
m_type = 2;
|
||||
m_upper = m_mupper;
|
||||
m_mlower = m_lower;
|
||||
}
|
||||
|
||||
void IDA_Solver::setMaxOrder(int n)
|
||||
{
|
||||
m_maxord = n;
|
||||
}
|
||||
|
||||
void IDA_Solver::setMaxNumSteps(int n)
|
||||
{
|
||||
m_maxsteps = n;
|
||||
}
|
||||
|
||||
void IDA_Solver::setInitialStepSize(doublereal h0)
|
||||
{
|
||||
m_h0 = h0;
|
||||
}
|
||||
|
||||
void IDA_Solver::setStopTime(doublereal tstop)
|
||||
{
|
||||
m_tstop = tstop;
|
||||
}
|
||||
|
||||
doublereal IDA_Solver::getCurrentStepFromIDA()
|
||||
{
|
||||
doublereal hcur;
|
||||
IDAGetCurrentStep(m_ida_mem, &hcur);
|
||||
return hcur;
|
||||
}
|
||||
|
||||
void IDA_Solver::setJacobianType(int formJac)
|
||||
{
|
||||
m_formJac = formJac;
|
||||
if (m_ida_mem && m_formJac == 1) {
|
||||
#if CT_SUNDIALS_VERSION >= 60
|
||||
int flag = IDASetJacFn(m_ida_mem, ida_jacobian);
|
||||
#else
|
||||
int flag = IDADlsSetJacFn(m_ida_mem, ida_jacobian);
|
||||
#endif
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::setJacobianType",
|
||||
"IDADlsSetDenseJacFn failed.");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void IDA_Solver::setMaxErrTestFailures(int maxErrTestFails)
|
||||
{
|
||||
m_maxErrTestFails = maxErrTestFails;
|
||||
}
|
||||
|
||||
void IDA_Solver::setMaxNonlinIterations(int n)
|
||||
{
|
||||
m_maxNonlinIters = n;
|
||||
}
|
||||
|
||||
void IDA_Solver::setMaxNonlinConvFailures(int n)
|
||||
{
|
||||
m_maxNonlinConvFails = n;
|
||||
}
|
||||
|
||||
void IDA_Solver::inclAlgebraicInErrorTest(bool yesno)
|
||||
{
|
||||
if (yesno) {
|
||||
m_setSuppressAlg = 0;
|
||||
} else {
|
||||
m_setSuppressAlg = 1;
|
||||
}
|
||||
}
|
||||
|
||||
void IDA_Solver::init(doublereal t0)
|
||||
{
|
||||
m_t0 = t0;
|
||||
m_told = t0;
|
||||
m_told_old = t0;
|
||||
m_tcurrent = t0;
|
||||
if (m_y) {
|
||||
N_VDestroy_Serial(m_y);
|
||||
}
|
||||
if (m_ydot) {
|
||||
N_VDestroy_Serial(m_ydot);
|
||||
}
|
||||
if (m_id) {
|
||||
N_VDestroy_Serial(m_id);
|
||||
}
|
||||
if (m_constraints) {
|
||||
N_VDestroy_Serial(m_constraints);
|
||||
}
|
||||
|
||||
#if CT_SUNDIALS_VERSION >= 60
|
||||
m_y = N_VNew_Serial(m_neq, m_sundials_ctx.get());
|
||||
m_ydot = N_VNew_Serial(m_neq, m_sundials_ctx.get());
|
||||
m_constraints = N_VNew_Serial(m_neq, m_sundials_ctx.get());
|
||||
#else
|
||||
m_y = N_VNew_Serial(m_neq);
|
||||
m_ydot = N_VNew_Serial(m_neq);
|
||||
m_constraints = N_VNew_Serial(m_neq);
|
||||
#endif
|
||||
|
||||
for (int i=0; i<m_neq; i++) {
|
||||
NV_Ith_S(m_y, i) = 0.0;
|
||||
NV_Ith_S(m_ydot, i) = 0.0;
|
||||
NV_Ith_S(m_constraints, i) = 0.0;
|
||||
}
|
||||
|
||||
// get the initial conditions
|
||||
m_resid.getInitialConditions(m_t0, NV_DATA_S(m_y), NV_DATA_S(m_ydot));
|
||||
|
||||
if (m_ida_mem) {
|
||||
IDAFree(&m_ida_mem);
|
||||
}
|
||||
|
||||
/* Call IDACreate */
|
||||
#if CT_SUNDIALS_VERSION >= 60
|
||||
m_ida_mem = IDACreate(m_sundials_ctx.get());
|
||||
#else
|
||||
m_ida_mem = IDACreate();
|
||||
#endif
|
||||
|
||||
int flag = 0;
|
||||
if (m_itol == IDA_SV) {
|
||||
flag = IDAInit(m_ida_mem, ida_resid, m_t0, m_y, m_ydot);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
if (flag == IDA_MEM_FAIL) {
|
||||
throw CanteraError("IDA_Solver::init",
|
||||
"Memory allocation failed.");
|
||||
} else if (flag == IDA_ILL_INPUT) {
|
||||
throw CanteraError("IDA_Solver::init",
|
||||
"Illegal value for IDAInit input argument.");
|
||||
} else {
|
||||
throw CanteraError("IDA_Solver::init", "IDAInit failed.");
|
||||
}
|
||||
}
|
||||
flag = IDASVtolerances(m_ida_mem, m_reltol, m_abstol);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::init", "Memory allocation failed.");
|
||||
}
|
||||
} else {
|
||||
flag = IDAInit(m_ida_mem, ida_resid, m_t0, m_y, m_ydot);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
if (flag == IDA_MEM_FAIL) {
|
||||
throw CanteraError("IDA_Solver::init",
|
||||
"Memory allocation failed.");
|
||||
} else if (flag == IDA_ILL_INPUT) {
|
||||
throw CanteraError("IDA_Solver::init",
|
||||
"Illegal value for IDAInit input argument.");
|
||||
} else {
|
||||
throw CanteraError("IDA_Solver::init", "IDAInit failed.");
|
||||
}
|
||||
}
|
||||
flag = IDASStolerances(m_ida_mem, m_reltol, m_abstols);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::init", "Memory allocation failed.");
|
||||
}
|
||||
}
|
||||
|
||||
// set the linear solver type
|
||||
if (m_type == 1 || m_type == 0) {
|
||||
long int N = m_neq;
|
||||
int flag;
|
||||
SUNLinSolFree((SUNLinearSolver) m_linsol);
|
||||
SUNMatDestroy((SUNMatrix) m_linsol_matrix);
|
||||
#if CT_SUNDIALS_VERSION >= 60
|
||||
m_linsol_matrix = SUNDenseMatrix(N, N, m_sundials_ctx.get());
|
||||
#else
|
||||
m_linsol_matrix = SUNDenseMatrix(N, N);
|
||||
#endif
|
||||
if (m_linsol_matrix == nullptr) {
|
||||
throw CanteraError("IDA_Solver::init",
|
||||
"Unable to create SUNDenseMatrix of size {0} x {0}", N);
|
||||
}
|
||||
#if CT_SUNDIALS_VERSION >= 60
|
||||
#if CT_SUNDIALS_USE_LAPACK
|
||||
m_linsol = SUNLinSol_LapackDense(m_y, (SUNMatrix) m_linsol_matrix,
|
||||
m_sundials_ctx.get());
|
||||
#else
|
||||
m_linsol = SUNLinSol_Dense(m_y, (SUNMatrix) m_linsol_matrix,
|
||||
m_sundials_ctx.get());
|
||||
#endif
|
||||
flag = IDASetLinearSolver(m_ida_mem, (SUNLinearSolver) m_linsol,
|
||||
(SUNMatrix) m_linsol_matrix);
|
||||
#else
|
||||
#if CT_SUNDIALS_USE_LAPACK
|
||||
m_linsol = SUNLapackDense(m_y, (SUNMatrix) m_linsol_matrix);
|
||||
#else
|
||||
m_linsol = SUNDenseLinearSolver(m_y, (SUNMatrix) m_linsol_matrix);
|
||||
#endif
|
||||
flag = IDADlsSetLinearSolver(m_ida_mem, (SUNLinearSolver) m_linsol,
|
||||
(SUNMatrix) m_linsol_matrix);
|
||||
#endif
|
||||
|
||||
if (flag) {
|
||||
throw CanteraError("IDA_Solver::init", "IDADense failed");
|
||||
}
|
||||
} else if (m_type == 2) {
|
||||
long int N = m_neq;
|
||||
long int nu = m_mupper;
|
||||
long int nl = m_mlower;
|
||||
SUNLinSolFree((SUNLinearSolver) m_linsol);
|
||||
SUNMatDestroy((SUNMatrix) m_linsol_matrix);
|
||||
#if CT_SUNDIALS_VERSION >= 60
|
||||
m_linsol_matrix = SUNBandMatrix(N, nu, nl, m_sundials_ctx.get());
|
||||
#elif CT_SUNDIALS_VERSION >= 40
|
||||
m_linsol_matrix = SUNBandMatrix(N, nu, nl);
|
||||
#else
|
||||
m_linsol_matrix = SUNBandMatrix(N, nu, nl, nu+nl);
|
||||
#endif
|
||||
if (m_linsol_matrix == nullptr) {
|
||||
throw CanteraError("IDA_Solver::init",
|
||||
"Unable to create SUNBandMatrix of size {} with bandwidths "
|
||||
"{} and {}", N, nu, nl);
|
||||
}
|
||||
#if CT_SUNDIALS_VERSION >= 60
|
||||
#if CT_SUNDIALS_USE_LAPACK
|
||||
m_linsol = SUNLinSol_LapackBand(m_y, (SUNMatrix) m_linsol_matrix,
|
||||
m_sundials_ctx.get());
|
||||
#else
|
||||
m_linsol = SUNLinSol_Band(m_y, (SUNMatrix) m_linsol_matrix,
|
||||
m_sundials_ctx.get());
|
||||
#endif
|
||||
IDASetLinearSolver(m_ida_mem, (SUNLinearSolver) m_linsol,
|
||||
(SUNMatrix) m_linsol_matrix);
|
||||
#else
|
||||
#if CT_SUNDIALS_USE_LAPACK
|
||||
m_linsol = SUNLapackBand(m_y, (SUNMatrix) m_linsol_matrix);
|
||||
#else
|
||||
m_linsol = SUNBandLinearSolver(m_y, (SUNMatrix) m_linsol_matrix);
|
||||
#endif
|
||||
IDADlsSetLinearSolver(m_ida_mem, (SUNLinearSolver) m_linsol,
|
||||
(SUNMatrix) m_linsol_matrix);
|
||||
#endif
|
||||
} else {
|
||||
throw CanteraError("IDA_Solver::init",
|
||||
"unsupported linear solver type");
|
||||
}
|
||||
|
||||
if (m_formJac == 1) {
|
||||
#if CT_SUNDIALS_VERSION >= 60
|
||||
flag = IDASetJacFn(m_ida_mem, ida_jacobian);
|
||||
#else
|
||||
flag = IDADlsSetJacFn(m_ida_mem, ida_jacobian);
|
||||
#endif
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::init",
|
||||
"IDADlsSetDenseJacFn failed.");
|
||||
}
|
||||
}
|
||||
|
||||
// pass a pointer to func in m_data
|
||||
m_fdata = make_unique<ResidData>(&m_resid, this, m_resid.nparams());
|
||||
flag = IDASetUserData(m_ida_mem, m_fdata.get());
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::init", "IDASetUserData failed.");
|
||||
}
|
||||
|
||||
// set options
|
||||
if (m_maxord > 0) {
|
||||
flag = IDASetMaxOrd(m_ida_mem, m_maxord);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::init", "IDASetMaxOrd failed.");
|
||||
}
|
||||
}
|
||||
if (m_maxsteps > 0) {
|
||||
flag = IDASetMaxNumSteps(m_ida_mem, m_maxsteps);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::init", "IDASetMaxNumSteps failed.");
|
||||
}
|
||||
}
|
||||
if (m_h0 > 0.0) {
|
||||
flag = IDASetInitStep(m_ida_mem, m_h0);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::init", "IDASetInitStep failed.");
|
||||
}
|
||||
}
|
||||
if (m_tstop > 0.0) {
|
||||
flag = IDASetStopTime(m_ida_mem, m_tstop);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::init", "IDASetStopTime failed.");
|
||||
}
|
||||
}
|
||||
if (m_maxErrTestFails >= 0) {
|
||||
flag = IDASetMaxErrTestFails(m_ida_mem, m_maxErrTestFails);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::init",
|
||||
"IDASetMaxErrTestFails failed.");
|
||||
}
|
||||
}
|
||||
if (m_maxNonlinIters > 0) {
|
||||
flag = IDASetMaxNonlinIters(m_ida_mem, m_maxNonlinIters);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::init",
|
||||
"IDASetmaxNonlinIters failed.");
|
||||
}
|
||||
}
|
||||
if (m_maxNonlinConvFails >= 0) {
|
||||
flag = IDASetMaxConvFails(m_ida_mem, m_maxNonlinConvFails);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::init",
|
||||
"IDASetMaxConvFails failed.");
|
||||
}
|
||||
}
|
||||
if (m_setSuppressAlg != 0) {
|
||||
flag = IDASetSuppressAlg(m_ida_mem, m_setSuppressAlg);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::init", "IDASetSuppressAlg failed.");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void IDA_Solver::correctInitial_Y_given_Yp(doublereal* y, doublereal* yp, doublereal tout)
|
||||
{
|
||||
doublereal tout1 = tout;
|
||||
if (tout == 0.0) {
|
||||
double h0 = 1.0E-5;
|
||||
if (m_h0 > 0.0) {
|
||||
h0 = m_h0;
|
||||
}
|
||||
tout1 = m_t0 + h0;
|
||||
}
|
||||
|
||||
int flag = IDACalcIC(m_ida_mem, IDA_Y_INIT, tout1);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::correctInitial_Y_given_Yp",
|
||||
"IDACalcIC failed: error = {}", flag);
|
||||
}
|
||||
|
||||
flag = IDAGetConsistentIC(m_ida_mem, m_y, m_ydot);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::correctInitial_Y_given_Yp",
|
||||
"IDAGetSolution failed: error = {}", flag);
|
||||
}
|
||||
for (int i = 0; i < m_neq; i++) {
|
||||
y[i] = NV_Ith_S(m_y, i);
|
||||
yp[i] = NV_Ith_S(m_ydot, i);
|
||||
}
|
||||
}
|
||||
|
||||
void IDA_Solver::correctInitial_YaYp_given_Yd(doublereal* y, doublereal* yp, doublereal tout)
|
||||
{
|
||||
int icopt = IDA_YA_YDP_INIT;
|
||||
doublereal tout1 = tout;
|
||||
if (tout == 0.0) {
|
||||
double h0 = 1.0E-5;
|
||||
if (m_h0 > 0.0) {
|
||||
h0 = m_h0;
|
||||
}
|
||||
tout1 = m_t0 + h0;
|
||||
}
|
||||
|
||||
int flag = IDACalcIC(m_ida_mem, icopt, tout1);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::correctInitial_YaYp_given_Yd",
|
||||
"IDACalcIC failed: error = {}", flag);
|
||||
}
|
||||
|
||||
flag = IDAGetConsistentIC(m_ida_mem, m_y, m_ydot);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::correctInitial_YaYp_given_Yd",
|
||||
"IDAGetSolution failed: error = {}", flag);
|
||||
}
|
||||
for (int i = 0; i < m_neq; i++) {
|
||||
y[i] = NV_Ith_S(m_y, i);
|
||||
yp[i] = NV_Ith_S(m_ydot, i);
|
||||
}
|
||||
}
|
||||
|
||||
int IDA_Solver::solve(double tout)
|
||||
{
|
||||
double tretn = tout - 1000;
|
||||
int flag;
|
||||
flag = IDASetStopTime(m_ida_mem, tout);
|
||||
if (flag != IDA_SUCCESS) {
|
||||
throw CanteraError("IDA_Solver::solve", "IDA error encountered.");
|
||||
}
|
||||
while (tretn < tout) {
|
||||
if (tout <= m_tcurrent) {
|
||||
throw CanteraError("IDA_Solver::solve", "tout <= tcurrent");
|
||||
}
|
||||
m_told_old = m_told;
|
||||
m_told = m_tcurrent;
|
||||
flag = IDASolve(m_ida_mem, tout, &tretn, m_y, m_ydot, IDA_ONE_STEP);
|
||||
if (flag < 0) {
|
||||
throw CanteraError("IDA_Solver::solve", "IDA error encountered.");
|
||||
} else if (flag == IDA_TSTOP_RETURN) {
|
||||
// we've reached our goal, and have actually integrated past it
|
||||
} else if (flag == IDA_ROOT_RETURN) {
|
||||
// not sure what to do with this yet
|
||||
} else if (flag == IDA_WARNING) {
|
||||
throw CanteraError("IDA_Solver::solve", "IDA Warning encountered.");
|
||||
}
|
||||
m_tcurrent = tretn;
|
||||
m_deltat = m_tcurrent - m_told;
|
||||
};
|
||||
|
||||
if (flag != IDA_SUCCESS && flag != IDA_TSTOP_RETURN) {
|
||||
throw CanteraError("IDA_Solver::solve", "IDA error encountered.");
|
||||
}
|
||||
return flag;
|
||||
}
|
||||
|
||||
double IDA_Solver::step(double tout)
|
||||
{
|
||||
double t;
|
||||
int flag;
|
||||
if (tout <= m_tcurrent) {
|
||||
throw CanteraError("IDA_Solver::step", "tout <= tcurrent");
|
||||
}
|
||||
m_told_old = m_told;
|
||||
m_told = m_tcurrent;
|
||||
flag = IDASolve(m_ida_mem, tout, &t, m_y, m_ydot, IDA_ONE_STEP);
|
||||
if (flag < 0) {
|
||||
throw CanteraError("IDA_Solver::step", "IDA error encountered.");
|
||||
} else if (flag == IDA_TSTOP_RETURN) {
|
||||
// we've reached our goal, and have actually integrated past it
|
||||
} else if (flag == IDA_ROOT_RETURN) {
|
||||
// not sure what to do with this yet
|
||||
} else if (flag == IDA_WARNING) {
|
||||
throw CanteraError("IDA_Solver::step", "IDA Warning encountered.");
|
||||
}
|
||||
m_tcurrent = t;
|
||||
m_deltat = m_tcurrent - m_told;
|
||||
return t;
|
||||
}
|
||||
|
||||
doublereal IDA_Solver::getOutputParameter(int flag) const
|
||||
{
|
||||
long int lenrw, leniw;
|
||||
switch (flag) {
|
||||
case REAL_WORKSPACE_SIZE:
|
||||
flag = IDAGetWorkSpace(m_ida_mem, &lenrw, &leniw);
|
||||
return doublereal(lenrw);
|
||||
break;
|
||||
}
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
}
|
||||
@ -1,4 +1,4 @@
|
||||
//! @file ODE_integrators.cpp
|
||||
//! @file Integrators.cpp
|
||||
|
||||
// This file is part of Cantera. See License.txt in the top-level directory or
|
||||
// at https://cantera.org/license.txt for license and copyright information.
|
||||
@ -6,6 +6,7 @@
|
||||
#include "cantera/base/ct_defs.h"
|
||||
#include "cantera/numerics/Integrator.h"
|
||||
#include "cantera/numerics/CVodesIntegrator.h"
|
||||
#include "cantera/numerics/IDAIntegrator.h"
|
||||
|
||||
namespace Cantera
|
||||
{
|
||||
@ -14,9 +15,11 @@ Integrator* newIntegrator(const std::string& itype)
|
||||
{
|
||||
if (itype == "CVODE") {
|
||||
return new CVodesIntegrator();
|
||||
} else if (itype == "IDA") {
|
||||
return new IDAIntegrator();
|
||||
} else {
|
||||
throw CanteraError("newIntegrator",
|
||||
"unknown ODE integrator: "+itype);
|
||||
"unknown integrator: "+itype);
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
@ -1,4 +1,5 @@
|
||||
//! @file FlowReactor.cpp A steady-state plug flow reactor
|
||||
//! @file FlowReactor.cpp A steady-state, ideal-gas, adiabatic,
|
||||
//! constant-area (cylindrical), frictionless plug flow reactor
|
||||
|
||||
// This file is part of Cantera. See License.txt in the top-level directory or
|
||||
// at https://cantera.org/license.txt for license and copyright information.
|
||||
@ -7,86 +8,351 @@
|
||||
#include "cantera/base/global.h"
|
||||
#include "cantera/kinetics/Kinetics.h"
|
||||
#include "cantera/thermo/ThermoPhase.h"
|
||||
#include "cantera/zeroD/ReactorSurface.h"
|
||||
#include "cantera/thermo/SurfPhase.h"
|
||||
#include "cantera/numerics/DenseMatrix.h"
|
||||
#include "cantera/kinetics/InterfaceKinetics.h"
|
||||
|
||||
using namespace std;
|
||||
|
||||
namespace Cantera
|
||||
{
|
||||
|
||||
void FlowReactor::getState(double* y)
|
||||
FlowReactor::FlowReactor(double area, double sa_to_vol,
|
||||
double ss_atol, double ss_rtol,
|
||||
int max_ss_steps,
|
||||
int max_ss_error_fails) :
|
||||
m_u(0.0),
|
||||
m_T(0.0),
|
||||
m_area(area),
|
||||
m_sa_to_vol(sa_to_vol),
|
||||
m_sdot_temp(0),
|
||||
m_hk(0),
|
||||
m_ss_rtol(ss_rtol),
|
||||
m_ss_atol(ss_atol),
|
||||
m_max_ss_steps(max_ss_steps),
|
||||
m_max_ss_error_fails(max_ss_error_fails)
|
||||
{
|
||||
}
|
||||
|
||||
void FlowReactor::getState(double* y, double* ydot)
|
||||
{
|
||||
if (m_thermo == 0) {
|
||||
throw CanteraError("FlowReactor::getState",
|
||||
"Error: reactor is empty.");
|
||||
}
|
||||
m_thermo->restoreState(m_state);
|
||||
m_thermo->getMassFractions(y+2);
|
||||
y[0] = 0.0; // distance
|
||||
m_thermo->getMassFractions(y+m_non_spec_eq);
|
||||
const vector_fp& mw = m_thermo->molecularWeights();
|
||||
|
||||
// set the first component to the initial density
|
||||
y[0] = m_rho;
|
||||
|
||||
// set the second component to the initial speed
|
||||
y[1] = m_speed0;
|
||||
y[1] = m_u;
|
||||
|
||||
// set the fourth component to the initial pressure
|
||||
y[2] = m_P;
|
||||
|
||||
// set the third component to the initial temperature
|
||||
y[3] = m_T;
|
||||
|
||||
if (m_chem) {
|
||||
m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
|
||||
}
|
||||
|
||||
// need to advance the reactor surface to steady state to get the initial
|
||||
// coverages
|
||||
for (auto &m_surf : m_surfaces)
|
||||
{
|
||||
auto *kin = static_cast<InterfaceKinetics*>(m_surf->kinetics());
|
||||
kin->advanceCoverages(100.0, m_ss_rtol, m_ss_atol, 0, m_max_ss_steps,
|
||||
m_max_ss_error_fails);
|
||||
auto *surf = static_cast<SurfPhase*>(&kin->thermo(kin->reactionPhaseIndex()));
|
||||
vector_fp cov(surf->nSpecies());
|
||||
surf->getCoverages(&cov[0]);
|
||||
m_surf->setCoverages(&cov[0]);
|
||||
m_surf->syncCoverages();
|
||||
}
|
||||
|
||||
// set the initial coverages
|
||||
getSurfaceInitialConditions(y + m_non_spec_eq + m_nsp);
|
||||
|
||||
// reset ydot vector
|
||||
std::fill(ydot, ydot + m_nv, 0.0);
|
||||
// calculate initial d(coverage)/dt values
|
||||
evalSurfaces(0, ydot + m_non_spec_eq + m_nsp);
|
||||
|
||||
// next, we must solve for the initial derivative values
|
||||
// a . ydot = b
|
||||
// ydot -> [rho', u', P', T', Yk']
|
||||
// a -> coefficients of [rho', u', P', T', Yk'] in the governing eq's, given
|
||||
// the initial state
|
||||
// b -> rhs constant of each conservation equation
|
||||
//
|
||||
// note that the species coverages are included in the algebraic constraints,
|
||||
// hence are not included ehre
|
||||
DenseMatrix a;
|
||||
a.resize(m_non_spec_eq + m_nsp, m_non_spec_eq + m_nsp, 0.0);
|
||||
|
||||
// first row is the ideal gas equation
|
||||
a(0, 0) = - GasConstant * m_T / m_thermo->meanMolecularWeight();
|
||||
a(0, 2) = 1;
|
||||
a(0, 3) = - m_rho * GasConstant / m_thermo->meanMolecularWeight();
|
||||
for (size_t i = 0; i < m_nsp; ++i)
|
||||
{
|
||||
a(0, m_non_spec_eq + i) = - m_rho * m_T / mw[i];
|
||||
}
|
||||
|
||||
|
||||
// initialize the second row from mass conservation equation
|
||||
// Kee 16.48
|
||||
a(1, 0) = m_u; // rho' * u
|
||||
a(1, 1) = m_rho; // u' * rho
|
||||
|
||||
// initialize the third row from momentum conservation equation
|
||||
// Kee 16.62
|
||||
// -> rho * u * u' + P' = -u * (P / A_c) * sum(sdot_k * W_k)
|
||||
a(2, 1) = m_rho * m_u; // rho * u * u'
|
||||
a(2, 2) = 1; // 1 * P'
|
||||
|
||||
// initialize the fourth row from conservation of energy
|
||||
// Kee 16.58, adiabatic
|
||||
const doublereal cp_mass = m_thermo->cp_mass();
|
||||
a(3, 3) = m_rho * m_u * cp_mass; // rho * u * cp * T'
|
||||
|
||||
// initialize the next rows from the mass-fraction equations
|
||||
// Kee 16.51
|
||||
for (size_t i = 0; i < m_nsp; ++i)
|
||||
{
|
||||
a(m_non_spec_eq + i, m_non_spec_eq + i) = m_rho * m_u; // rho * u * Yk'
|
||||
}
|
||||
|
||||
// now set the RHS vector
|
||||
|
||||
// get (perim / Ac) * sum(sk' * wk), used multiple places
|
||||
double h_sk_wk = 0;
|
||||
const doublereal hydraulic = surfaceAreaToVolumeRatio();
|
||||
for (size_t i = 0; i < m_nsp; ++i)
|
||||
{
|
||||
h_sk_wk += hydraulic * m_sdot[i] * mw[i];
|
||||
}
|
||||
|
||||
// RHS of ideal gas eq. is zero
|
||||
ydot[0] = 0;
|
||||
|
||||
// mass continutity, Kee 16.48
|
||||
// (perim / Ac) * sum(sk' * wk)
|
||||
ydot[1] = h_sk_wk;
|
||||
|
||||
// momentum conservation, Kee 16.62, no friction
|
||||
// - u * (perim / Ac) * sum(sk' * wk)
|
||||
ydot[2] = -m_u * h_sk_wk;
|
||||
|
||||
// energy conservation, Kee 16.58, adiabatic
|
||||
// -sum(wk' * Wk * hk) - (perim / Ac) * sum(sk' * Wk * hk)
|
||||
// Note: the partial molar enthalpies are in molar units, while Kee uses mass
|
||||
// units, where:
|
||||
// h_mass = h_mole / Wk
|
||||
// hence:
|
||||
// h_mass * Wk = h_mol
|
||||
m_thermo->getPartialMolarEnthalpies(&m_hk[0]);
|
||||
for (size_t i = 0; i < m_nsp; ++i)
|
||||
{
|
||||
ydot[3] -= m_hk[i] * (m_wdot[i] + hydraulic * m_sdot[i]);
|
||||
}
|
||||
|
||||
// mass-fraction equations Kee 16.51
|
||||
// - Yk * (perim / Ac) * sum(sk' * Wk) + wk' * Wk + (perim / Ac) * sk' * Wk
|
||||
for (size_t i = 0; i < m_nsp; ++i)
|
||||
{
|
||||
ydot[m_non_spec_eq + i] = -y[m_non_spec_eq + i] * h_sk_wk +
|
||||
mw[i] * (m_wdot[i] + hydraulic * m_sdot[i]);
|
||||
}
|
||||
|
||||
// and solve
|
||||
solve(a, ydot, 1, 0);
|
||||
}
|
||||
|
||||
|
||||
void FlowReactor::initialize(doublereal t0)
|
||||
{
|
||||
Reactor::initialize(t0);
|
||||
m_thermo->restoreState(m_state);
|
||||
m_nv = m_nsp + 2;
|
||||
// initialize state
|
||||
m_T = m_thermo->temperature();
|
||||
m_rho = m_thermo->density();
|
||||
m_P = m_thermo->pressure();
|
||||
m_T = m_thermo->temperature();
|
||||
// resize temporary arrays
|
||||
m_wdot.resize(m_nsp);
|
||||
m_hk.resize(m_nsp);
|
||||
// set number of variables to the number of non-species equations
|
||||
// i.e., density, velocity, pressure and temperature
|
||||
// plus the number of species in the gas phase
|
||||
m_nv = m_non_spec_eq + m_nsp;
|
||||
if (m_surfaces.size())
|
||||
{
|
||||
size_t n_surf_species = 0;
|
||||
size_t n_total_species = 0;
|
||||
for (auto const &m_surf : m_surfaces)
|
||||
{
|
||||
// finally, add the number of species the surface for the site-fraction
|
||||
// equations
|
||||
n_surf_species += m_surf->thermo()->nSpecies();
|
||||
n_total_species += m_surf->kinetics()->nTotalSpecies();
|
||||
}
|
||||
m_nv += n_surf_species;
|
||||
m_sdot_temp.resize(n_total_species);
|
||||
}
|
||||
}
|
||||
|
||||
void FlowReactor::syncState()
|
||||
{
|
||||
ReactorBase::syncState();
|
||||
m_rho = m_thermo->density();
|
||||
m_P = m_thermo->pressure();
|
||||
m_T = m_thermo->temperature();
|
||||
}
|
||||
|
||||
void FlowReactor::updateState(doublereal* y)
|
||||
{
|
||||
// Set the mass fractions and density of the mixture.
|
||||
m_dist = y[0];
|
||||
m_speed = y[1];
|
||||
doublereal* mss = y + 2;
|
||||
m_thermo->setMassFractions(mss);
|
||||
doublereal rho = m_rho0 * m_speed0/m_speed;
|
||||
m_rho = y[0];
|
||||
m_u = y[1];
|
||||
m_P = y[2];
|
||||
m_T = y[3];
|
||||
doublereal* mss = y + m_non_spec_eq;
|
||||
m_thermo->setMassFractions_NoNorm(mss);
|
||||
m_thermo->setState_TP(m_T, m_P);
|
||||
|
||||
// assumes frictionless
|
||||
doublereal pmom = m_P0 - rho*m_speed*m_speed;
|
||||
// update surface
|
||||
updateSurfaceState(y + m_nsp + m_non_spec_eq);
|
||||
|
||||
doublereal hmom;
|
||||
// assumes adiabatic
|
||||
if (m_energy) {
|
||||
hmom = m_h0 - 0.5*m_speed*m_speed;
|
||||
m_thermo->setState_HP(hmom, pmom);
|
||||
} else {
|
||||
m_thermo->setState_TP(m_T, pmom);
|
||||
}
|
||||
m_thermo->saveState(m_state);
|
||||
}
|
||||
|
||||
void FlowReactor::setMassFlowRate(double mdot)
|
||||
{
|
||||
m_rho0 = m_thermo->density();
|
||||
m_speed = mdot/m_rho0;
|
||||
m_speed0 = m_speed;
|
||||
m_u = mdot/m_rho0;
|
||||
m_u0 = m_u;
|
||||
m_T = m_thermo->temperature();
|
||||
m_P0 = m_thermo->pressure() + m_rho0*m_speed*m_speed;
|
||||
m_h0 = m_thermo->enthalpy_mass() + 0.5*m_speed*m_speed;
|
||||
m_P0 = m_thermo->pressure() + m_rho0*m_u*m_u;
|
||||
m_h0 = m_thermo->enthalpy_mass() + 0.5*m_u*m_u;
|
||||
}
|
||||
|
||||
void FlowReactor::eval(double time, double* LHS, double* RHS)
|
||||
void FlowReactor::updateSurfaceState(double* y)
|
||||
{
|
||||
size_t loc = 0;
|
||||
for (auto& S : m_surfaces) {
|
||||
// set the coverages without normalization to avoid over-constraining
|
||||
// the system.
|
||||
// note: the ReactorSurface class doesn't normalize when calling setCoverages
|
||||
S->setCoverages(y+loc);
|
||||
S->syncCoverages();
|
||||
loc += S->thermo()->nSpecies();
|
||||
}
|
||||
}
|
||||
|
||||
void FlowReactor::eval(double time, double* y,
|
||||
double* ydot, double* params,
|
||||
double* residual)
|
||||
{
|
||||
m_thermo->restoreState(m_state);
|
||||
|
||||
// distance equation
|
||||
RHS[0] = m_speed;
|
||||
|
||||
// speed equation. Set m_fctr to a large value, so that rho*u is held fixed
|
||||
RHS[1] = m_fctr * (m_speed0 - m_thermo->density() * m_speed / m_rho0);
|
||||
|
||||
// species equations //
|
||||
evalSurfaces(t, &m_sdot_temp[0]);
|
||||
const vector_fp& mw = m_thermo->molecularWeights();
|
||||
|
||||
if (m_chem) {
|
||||
m_kin->getNetProductionRates(RHS + 2); // "omega dot"
|
||||
} else {
|
||||
fill(RHS + 2, RHS + 2 + m_nsp, 0.0);
|
||||
double sk_wk = 0;
|
||||
for (size_t i = 0; i < m_nsp; ++i)
|
||||
{
|
||||
sk_wk = m_sdot[i] * mw[i];
|
||||
}
|
||||
double rrho = 1.0 / m_thermo->density();
|
||||
for (size_t n = 0; n < m_nsp; n++) {
|
||||
RHS[n+2] *= mw[n] * rrho;
|
||||
m_thermo->getPartialMolarEnthalpies(&m_hk[0]);
|
||||
// get net production
|
||||
if (m_chem) {
|
||||
m_kin->getNetProductionRates(&m_wdot[0]);
|
||||
}
|
||||
|
||||
// set dphi/dz variables
|
||||
doublereal m_drhodz = ydot[0];
|
||||
doublereal m_dudz = ydot[1];
|
||||
doublereal m_dPdz = ydot[2];
|
||||
doublereal m_dTdz = ydot[3];
|
||||
|
||||
// use equation of state for density residual
|
||||
const double R_specific = GasConstant / m_thermo->meanMolecularWeight(); // specific gas constant
|
||||
// P' - rho' * (R/M) * T - rho * (R/M) * T'
|
||||
residual[0] = m_dPdz - m_drhodz * R_specific * m_T - m_rho * R_specific * m_dTdz;
|
||||
// next, subtract off rho * T * sum(Yi' / Wi)
|
||||
for (size_t i = 0; i < m_nsp; ++i)
|
||||
{
|
||||
residual[0] -= m_rho * m_T * ydot[m_non_spec_eq + i] / mw[i];
|
||||
}
|
||||
|
||||
//! use mass continuity for velocity residual
|
||||
//! Kee.'s Chemically Reacting Flow, Eq. 16.48
|
||||
const doublereal hydraulic = surfaceAreaToVolumeRatio();
|
||||
residual[1] = m_u * m_drhodz + m_rho * m_dudz - sk_wk * hydraulic;
|
||||
|
||||
//! Use conservation of momentum for pressure residual
|
||||
//! Kee.'s Chemically Reacting Flow, Eq. 16.62
|
||||
//! [NOTE: friction is neglected]
|
||||
residual[2] = m_rho * m_u * m_dudz + m_u * hydraulic * sk_wk + m_dPdz;
|
||||
|
||||
//! use conservation of energy for temperature residual
|
||||
//! Kee.'s Chemically Reacting Flow, Eq. 16.58
|
||||
//! [NOTE: adiabatic]
|
||||
// Also, as above:
|
||||
// h_mass = h_mole / Wk
|
||||
// hence:
|
||||
// h_mass * Wk = h_mol
|
||||
const doublereal cp_mass = m_thermo->cp_mass();
|
||||
residual[3] = m_rho * m_u * cp_mass * m_dTdz;
|
||||
for (size_t i = 0; i < m_nsp; ++i)
|
||||
{
|
||||
residual[3] += m_hk[i] * (m_wdot[i] + hydraulic * m_sdot[i]);
|
||||
}
|
||||
|
||||
//! species conservation equations
|
||||
//! Kee.'s Chemically Reacting Flow, Eq. 16.51
|
||||
for (size_t i = 0; i < m_nsp; ++i)
|
||||
{
|
||||
residual[i + m_non_spec_eq] = m_rho * m_u * ydot[i + m_non_spec_eq] +
|
||||
y[i + m_non_spec_eq] * hydraulic * sk_wk -
|
||||
mw[i] * (m_wdot[i] + hydraulic * m_sdot[i]);
|
||||
}
|
||||
|
||||
// surface algebraic constraints
|
||||
if (m_surfaces.size())
|
||||
{
|
||||
size_t loc = m_non_spec_eq + m_nsp; // offset into residual vector
|
||||
for (auto &m_surf : m_surfaces)
|
||||
{
|
||||
Kinetics* kin = m_surf->kinetics();
|
||||
SurfPhase* surf = m_surf->thermo();
|
||||
size_t nk = surf->nSpecies();
|
||||
kin->getNetProductionRates(&m_sdot_temp[0]);
|
||||
size_t ns = kin->surfacePhaseIndex();
|
||||
size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
|
||||
doublereal sum = y[loc];
|
||||
for (size_t i = 1; i < nk; ++i)
|
||||
{
|
||||
//! net surface production rate residuals
|
||||
//! Kee.'s Chemically Reacting Flow, Eq. 16.63
|
||||
residual[loc + i] = m_sdot_temp[surfloc + i];
|
||||
//! next, evaluate the algebraic constraint to explicitly conserve
|
||||
//! surface site fractions.
|
||||
//! Kee.'s Chemically Reacting Flow, Eq. 16.64
|
||||
sum += y[loc + i];
|
||||
}
|
||||
// set residual of the first species sum - 1 to ensure
|
||||
// surface site conservation.
|
||||
// Note: the first species is selected to be consistent with
|
||||
// Reactor::evalSurfaces
|
||||
residual[loc] = sum - 1.0;
|
||||
loc += nk;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -95,11 +361,15 @@ size_t FlowReactor::componentIndex(const string& nm) const
|
||||
// check for a gas species name
|
||||
size_t k = m_thermo->speciesIndex(nm);
|
||||
if (k != npos) {
|
||||
return k + 2;
|
||||
} else if (nm == "X" || nm == "distance") {
|
||||
return k + m_non_spec_eq;
|
||||
} else if (nm == "rho" || nm == "density") {
|
||||
return 0;
|
||||
} else if (nm == "U" || nm == "velocity") {
|
||||
return 1;
|
||||
} else if (nm == "P" || nm == "pressure") {
|
||||
return 2;
|
||||
} else if (nm == "T" || nm == "temperature") {
|
||||
return 3;
|
||||
} else {
|
||||
return npos;
|
||||
}
|
||||
|
||||
@ -9,24 +9,18 @@
|
||||
#include "cantera/base/utilities.h"
|
||||
#include "cantera/base/Array.h"
|
||||
#include "cantera/numerics/Integrator.h"
|
||||
#include "cantera/zeroD/FlowReactor.h"
|
||||
|
||||
#include <cstdio>
|
||||
|
||||
namespace Cantera
|
||||
{
|
||||
|
||||
ReactorNet::ReactorNet() :
|
||||
m_integ(newIntegrator("CVODE"))
|
||||
m_is_dae(false),
|
||||
m_nmax(0)
|
||||
{
|
||||
suppressErrors(true);
|
||||
|
||||
// use backward differencing, with a full Jacobian computed
|
||||
// numerically, and use a Newton linear iterator
|
||||
m_integ->setMethod(BDF_Method);
|
||||
m_integ->setLinearSolverType("DENSE");
|
||||
}
|
||||
|
||||
ReactorNet::~ReactorNet()
|
||||
{
|
||||
// Defined in .cpp to limit dependence on Integrator.h
|
||||
}
|
||||
|
||||
void ReactorNet::setInitialTime(double time)
|
||||
@ -69,9 +63,22 @@ void ReactorNet::setSensitivityTolerances(double rtol, double atol)
|
||||
m_init = false;
|
||||
}
|
||||
|
||||
void ReactorNet::setMaxSteps(int nmax)
|
||||
{
|
||||
m_nmax = nmax;
|
||||
m_init= false;
|
||||
}
|
||||
|
||||
void ReactorNet::initialize()
|
||||
{
|
||||
m_nv = 0;
|
||||
m_integ.reset(newIntegrator(m_is_dae ? "IDA" : "CVODE"));
|
||||
// use backward differencing, with a full Jacobian computed
|
||||
// numerically, and use a Newton linear iterator
|
||||
m_integ->setMethod(BDF_Method);
|
||||
m_integ->setProblemType(DENSE + NOJAC);
|
||||
m_integ->setIterator(Newton_Iter);
|
||||
|
||||
debuglog("Initializing reactor network.\n", m_verbose);
|
||||
if (m_reactors.empty()) {
|
||||
throw CanteraError("ReactorNet::initialize",
|
||||
@ -104,11 +111,13 @@ void ReactorNet::initialize()
|
||||
m_integ->setSensitivityTolerances(m_rtolsens, m_atolsens);
|
||||
m_integ->setMaxStepSize(m_maxstep);
|
||||
m_integ->setMaxErrTestFails(m_maxErrTestFails);
|
||||
if (m_nmax)
|
||||
m_integ->setMaxSteps(m_nmax);
|
||||
m_integ->initialize(m_time, *this);
|
||||
if (m_verbose) {
|
||||
writelog("Number of equations: {:d}\n", neq());
|
||||
writelog("Maximum time step: {:14.6g}\n", m_maxstep);
|
||||
}
|
||||
m_integ->initialize(m_time, *this);
|
||||
if (m_integ->preconditionerSide() != PreconditionerSide::NO_PRECONDITION) {
|
||||
checkPreconditionerSupported();
|
||||
}
|
||||
@ -226,7 +235,7 @@ double ReactorNet::step()
|
||||
} else if (!m_integrator_init) {
|
||||
reinitialize();
|
||||
}
|
||||
m_time = m_integ->step(m_time + 1.0);
|
||||
m_time = m_integ->step(m_time + 1.0);;
|
||||
updateState(m_integ->solution());
|
||||
return m_time;
|
||||
}
|
||||
@ -259,6 +268,19 @@ int ReactorNet::lastOrder()
|
||||
void ReactorNet::addReactor(Reactor& r)
|
||||
{
|
||||
r.setNetwork(this);
|
||||
try
|
||||
{
|
||||
dynamic_cast<FlowReactor&>(r);
|
||||
m_is_dae = true;
|
||||
}
|
||||
catch(bad_cast)
|
||||
{
|
||||
if (m_is_dae)
|
||||
{
|
||||
throw CanteraError("ReactorNet::addReactor",
|
||||
"Cannot mix Reactors and FlowReactor's");
|
||||
}
|
||||
}
|
||||
m_reactors.push_back(&r);
|
||||
}
|
||||
|
||||
@ -286,6 +308,24 @@ void ReactorNet::eval(doublereal t, doublereal* y,
|
||||
checkFinite("ydot", ydot, m_nv);
|
||||
}
|
||||
|
||||
void ReactorNet::eval(doublereal t, doublereal* y,
|
||||
doublereal* ydot, doublereal* p,
|
||||
doublereal* residual)
|
||||
{
|
||||
updateState(y);
|
||||
for (size_t n = 0; n < m_reactors.size(); n++) {
|
||||
m_reactors[n]->evalEqs(t, y + m_start[n], ydot + m_start[n], p, residual);
|
||||
}
|
||||
checkFinite("ydot", ydot, m_nv);
|
||||
}
|
||||
|
||||
void ReactorNet::getConstraints(doublereal* constraints)
|
||||
{
|
||||
for (size_t n = 0; n < m_reactors.size(); n++) {
|
||||
m_reactors[n]->getConstraints(constraints + m_start[n]);
|
||||
}
|
||||
}
|
||||
|
||||
double ReactorNet::sensitivity(size_t k, size_t p)
|
||||
{
|
||||
if (!m_init) {
|
||||
@ -295,7 +335,8 @@ double ReactorNet::sensitivity(size_t k, size_t p)
|
||||
throw IndexError("ReactorNet::sensitivity",
|
||||
"m_sens_params", p, m_sens_params.size()-1);
|
||||
}
|
||||
double denom = m_integ->solution(k);
|
||||
double denom;
|
||||
denom = m_integ->solution(k);
|
||||
if (denom == 0.0) {
|
||||
denom = SmallNumber;
|
||||
}
|
||||
@ -333,13 +374,6 @@ void ReactorNet::updateState(doublereal* y)
|
||||
}
|
||||
}
|
||||
|
||||
void ReactorNet::getState(double* y)
|
||||
{
|
||||
for (size_t n = 0; n < m_reactors.size(); n++) {
|
||||
m_reactors[n]->getState(y + m_start[n]);
|
||||
}
|
||||
}
|
||||
|
||||
void ReactorNet::getDerivative(int k, double* dky)
|
||||
{
|
||||
double* cvode_dky = m_integ->derivative(m_time, k);
|
||||
@ -374,6 +408,11 @@ bool ReactorNet::getAdvanceLimits(double *limits)
|
||||
has_limit |= m_reactors[n]->getAdvanceLimits(limits + m_start[n]);
|
||||
}
|
||||
return has_limit;
|
||||
void ReactorNet::getState(double* y, double* ydot)
|
||||
{
|
||||
for (size_t n = 0; n < m_reactors.size(); n++) {
|
||||
m_reactors[n]->getState(y + m_start[n], ydot + m_start[n]);
|
||||
}
|
||||
}
|
||||
|
||||
size_t ReactorNet::globalComponentIndex(const string& component, size_t reactor)
|
||||
|
||||
@ -1378,13 +1378,18 @@ class TestFlowReactor(utilities.CanteraTest):
|
||||
|
||||
t = 0
|
||||
v0 = r.speed
|
||||
times = [0]
|
||||
velocities = [r.speed]
|
||||
self.assertNear(v0, 10 / r.density)
|
||||
while t < 10.0:
|
||||
t = net.step()
|
||||
times.append(t)
|
||||
velocities.append(r.speed)
|
||||
|
||||
self.assertNear(v0, r.speed)
|
||||
self.assertNear(r.distance, v0 * t)
|
||||
self.assertNear(np.trapz(velocities, x=times), v0 * t, 1e-3)
|
||||
|
||||
@unittest.expectedFailure
|
||||
def test_reacting(self):
|
||||
g = ct.Solution(yaml=self.gas_def)
|
||||
g.TPX = 1400, 20*101325, 'CO:1.0, H2O:1.0'
|
||||
@ -1394,20 +1399,23 @@ class TestFlowReactor(utilities.CanteraTest):
|
||||
|
||||
net = ct.ReactorNet()
|
||||
net.add_reactor(r)
|
||||
net.atol = 1e-18
|
||||
net.rtol = 1e-9
|
||||
net.atol = 1e-10
|
||||
net.rtol = 1e-5
|
||||
net.max_err_test_fails = 10
|
||||
net.max_steps = int(1e6)
|
||||
|
||||
t = 0
|
||||
times = [0]
|
||||
velocities = [r.speed]
|
||||
v0 = r.speed
|
||||
self.assertNear(r.speed, 10 / r.density)
|
||||
while t < 1.0:
|
||||
t1 = net.time
|
||||
x1 = r.distance
|
||||
|
||||
t = net.step()
|
||||
times.append(t)
|
||||
velocities.append(r.speed)
|
||||
|
||||
v = (r.distance - x1) / (net.time - t1)
|
||||
self.assertNear(r.speed, v, 1e-3)
|
||||
self.assertNear(r.speed, v0, 1e-3)
|
||||
self.assertNear(np.trapz(velocities, x=times), v0 * t, 1e-3)
|
||||
|
||||
|
||||
class TestSurfaceKinetics(utilities.CanteraTest):
|
||||
|
||||
Loading…
Reference in New Issue
Block a user