You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: docs/src/Analysis.md
+40-7Lines changed: 40 additions & 7 deletions
Display the source diff
Display the rich diff
Original file line number
Diff line number
Diff line change
@@ -10,11 +10,19 @@ properties of the solution.
10
10
The Simulation object can be defined:
11
11
12
12
```
13
-
bsol = Simulation(sol,domain)
13
+
bsol = Simulation(sol,domain,interfaces,p)
14
14
```
15
15
16
-
where `sol` is the ODESolution object output by the DifferentialEquations package and `domain` is the
17
-
domain `sol` corresponds to.
16
+
where `sol` is the ODESolution object output by the DifferentialEquations package, `domain` is the
17
+
domain `sol` corresponds to, interfaces is the array of interface objects and p is the parameter vector.
18
+
19
+
## The SystemSimulation Object
20
+
21
+
When a system involves multiple domains a single Simulation object is insufficient. For these systems you need to construct a SystemSimulation object. This works in much the same way except that the domains should be listed as a tuple and the domain and interface ordering should be the same as that used when constructing the Reactor object. In theory the SystemSimulation object should be able to be used in place of a Simulation object in most places (If this is not the case and should be the case please make an issue!).
@@ -52,9 +60,21 @@ species for reactions this is the index of the reaction.
52
60
53
61
The function `rates` can be used to calculate the rates of all reactions at specific time points.
54
62
`rates(bsol,t)` will give the array of reaction rates at time `t`
55
-
while `rates(bsol;ts)` will give a matrix of reaction rates at all times in `ts`.
63
+
while `rates(bsol;ts=ts)` will give a matrix of reaction rates at all times in `ts`.
56
64
Note that `ts` defaults to `bsol.sol.t`.
57
65
66
+
### Adjoint Sensitivities
67
+
68
+
Sensitivity values to a target species or thermodynamic variable can be computed from a `Simulation` or `SystemSimulation` object using the `getadjointsensitivities(bsol::Q,target::String,solver::W;sensalg::W2=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(false)),abstol::Float64=1e-6,reltol::Float64=1e-3,normalize=true,kwargs...)` function. This computes the sensitivity with respect to the target at the final time point of bsol. It uses `solver`, `sensalg`, `abstol`, and `reltol` for the adjoint solve and by default will give the normalized sensitivity values (note these are molar sensitivities, concentration sensitivities can't be computed from a single adjoint pass). This is usually much faster than forward sensitivities.
69
+
70
+
### Transitory Sensitivities
71
+
72
+
Transitory sensitivity values can be computed using several different algorithms. `transitorysensitivitiesfullexact(sim::Simulation,t;tau=NaN,
abstol=1e-16,reltol=1e-6)` gives you the exact full matrix of transitory sensitivities using the forward sensitivity algorithm, while `transitorysensitivitiesfulltrapezoidal(sim,t;tau=NaN,normalized=true)` gives the approximate full matrix of transitory sensitivities using the trapezoidal method. `transitorysensitivitiesparamexact` and `transitorysensitivitiesparamtrapezoidal` are available for computing a single column of the matrix (with respect to a single parameter)). Lastly `transitorysensitivitiesadjointexact(sim::Simulation,t,name;tau=NaN,
abstol=1e-16,reltol=1e-6)` is available for computing a single row of the matrix (sensitivity to a specified species with respect to all parameters). The adjoint algorithm is jacobian free if `tau` is specified and the solver is jacobian free.
77
+
58
78
### Other Useful Properties
59
79
60
80
Please let us know on our Github issues page if we're missing any important
@@ -71,16 +91,20 @@ Mole fractions can be plotting using the `plotmolefractions` function
71
91
`tol` at the points in `bsol.sol.t`.
72
92
`plotmolefractions(bsol,spcnames)` plots all the species with names in `spcnames` at the points in `bsol.sol.t`.
73
93
74
-
### Plotting Concentration Sensitivity
94
+
### Plotting Forward Sensitivities
75
95
76
-
Concentration sensitivities can be plotted using the `plotmaxthermosensitivity` and `plotmaxratesensitivity` functions.
96
+
Sensitivities (normalized molar sensitivities) can be plotted using the `plotmaxthermoforwardsensitivity` and `plotmaxrateforwardsensitivity` functions.
`spcname` corresponds to the species sensitivities are being calculated for, `N` is the maximum number of
80
100
sensitive species/reactions plotted (0 corresponds to all of them), sensitive species/reactions with sensitivities
81
101
less than `tol` are not included in the plot. Note that the thermo sensitivities are given in mol/kcal while the rate
82
102
sensitivities are fully non-dimensionalized (as displayed on the plots).
83
103
104
+
### Plotting Adjoint Sensitivities
105
+
106
+
Adjoint sensitivities can be plotted using the `plotthermoadjointsensitivities(bsol::Y,name::X,dps::Z;N=0,tol=0.01)` and `plotrateadjointsensitivities(bsol::Y,name::X,dps::Z;N=0,tol=0.01)` functions where `dps` is the normalized adjoint sensitivity values.
107
+
84
108
### Plotting ROPs
85
109
86
110
ROPs can be plotted with the `plotrops` function.
@@ -98,6 +122,15 @@ than `tol` * the largest absolute rate at every time point is excluded from the
98
122
99
123
The analogous functions `plotradicalrops(bsol,t;N=0,tol=0.01)` and `plotradicalrops(bsol;rxnrates=Array{Float64,1}(),ts=Array{Float64,1}(),tol=0.05)` are available for plotting the rops for the sum of all radicals.
100
124
125
+
### Plotting Transitory Sensitivities
126
+
127
+
Transitory sensitivities and be combusted and plotted using `plotrxntransitorysensitivities(bsol,name,t;dSdt=nothing,tau=nothing,tol=1e-3,N=0,rxntol=1e-6)`
128
+
and `plotthermotransitorysensitivities(bsol,name,t;dSdt=nothing,tau=nothing,tol=1e-3,N=0)` where dSdt contains the transitory sensitivity values otherwise these values will be computed automatically using the trapezoidal method using the input `tau` or automatically selecting `tau` and the `rxntol` value. At most `N` reactions are included in the plot and the plot will not include reactions with absolute transitory sensitivities less than `tol` times the largest absolute value.
129
+
130
+
### Plotting Time Scales
131
+
132
+
The timescale distribution of a simulation at a point can be plotted using `plottimescales(sim,t;taumax=1e6,taumin=1e-18,taures=10.0^0.5,usediag=true)` or `plottimescales(Jy;taumax=1e6,taumin=1e-18,taures=10.0^0.5,usediag=true)` where `taumax`, `taumin` and `taures` control the bins. If `usediag=true` it will simply determine the timescales using the diagonal of the Jacobian otherwise it will compute and use the eigenvalues. In general we've observed no significant differences in distributions generated using the diagonal values vs the eigenvalues although there may be significant differences when the mechanism is small.
133
+
101
134
### Other Plots
102
135
103
136
While we are trying to build up a library of plotting functions that make mechanism analysis easier
Automatic mechanism analysis can be run for a single species at a single time point using the function `analyzespc(sim,spcname,t;N=10,tol=1e-3,branchthreshold=0.9,
)`. This returns an array of ReactionAnalysis objects corresponding to each reaction found in the analysis (based on `N` and `tol` choices).
10
+
11
+
## ReactionAnalysis Object
12
+
13
+
The ReactionAnalysis object has seven attributes, `branchings` the array of potentially important Branching objects, `paths` the array of potentially important ReactionPath objects, `radprodlossfract` the fraction of production (+) or loss (-) of radicals that this reaction accounts for, `spcind` the index of the target species, `spcname` the name of the target species, `rxnind` the index of the reaction and `sens` the transitory sensitivity value. Can be dumpped to a string report by `getrxnanalysisstring(sim,ra;branchingcutoff=1e-2,radbranchfract=0.01)` or simply printed with `printrxnanalysis(sim,ra;branchingcutoff=1e-2,radbranchfract=0.01)` where the `branchingcutoff` is the fraction of the branching at which reactions will no longer show up as part of a branching and `radbranchfract` is the fraction of radical production/loss above which the value is displayed as important.
14
+
15
+
## Branching Object
16
+
17
+
The Branching object has three attributes, `spcind` the index of the target species, `rxninds` the array of the indices of the reactions in order of branching fraction, `branchingratios` the fraction of the branching accounted for by each reaction. Can be used to generate a flux diagram with `getfluxdiagram(bsol,t,b::Branching; branchtol=1.0e-2, kwargs...)` where `branchtol` is the fraction of the branching at which the product of a reaction will not be included in the flux diagram. Note this flux diagram includes all reaction between species part of the important branching reactions and not just the branching reactions themselves.
18
+
19
+
## ReactionPath Object
20
+
21
+
The ReactionPath object has six important attributes, `forward` indices whether the path was generated following the flux forward from the target species or backwards from the target species, `spcsinds` is the array of the species indices followed in order along the flux path, `rxninds` is the array of reactions that connect these `spcsinds`, `spcind` is the index of the target species, `branchfracts` is the branching fraction of the reaction in `rxninds` for the species in `spcsinds`, `branchfract` is the fraction of the flux following from the start of the path to the end. Can be used to generate a flux diagram with `getfluxdiagram(bsol,t,rp::ReactionPath; radius=0, kwargs...)`. Note this flux diagram includes all reactions between species that are parts of the ReactionPath and not just the reaction path itself.
In many cases you will have connections between multiple domains in your system or between one domain and something outside the system. In these cases you need to define interfaces.
The Reactor object exists primarily to automatically sets up the ODEProblem object for you to solve.
65
-
It takes the domain, the initial condition vector returned when constructing the domainand a time interval.
80
+
For single domain reactors it takes the domain, the initial condition vector returned when constructing the domain, a time interval and an array of any interface objects and returns a Reactor object. If your system has multiple domains you need to pass the domains, initialconditions time interval, array of interfaces and parameters with the domains, associated initial conditions and associated parameters as tuples in consistent order from domains to interfaces.
RMS purposefully exposes the solver interface provide users with all the options available from
75
95
Julia's DifferentialEquations package. The ODEProblem object is a field of the Reactor
76
-
object `react.ode` and can be solved as the user desires.
96
+
object `react.ode` and can be solved as the user desires. A recommended solver choice is stored in `react.recommendedsolver`. User can also specify their own choice of solver.
97
+
98
+
Forward sensitivity analysis can also be requested on the Reactor object by setting `forwardsensitivities=true`. Note that adjoint sensitivity analysis is usually much faster and can be done as a postprocessing analysis after the simulation is complete without a need to set `forwardsensitivities=true` during the simulation (this is discussed in the Analysis section).
99
+
100
+
Example:
101
+
102
+
```
103
+
sol = solve(react.ode,react.recommendedsolver,abstol=1e-20,reltol=1e-12)
104
+
```
77
105
78
-
Example:
79
106
```
80
-
sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-12)
107
+
sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-12;forwardsensitivities=true)
81
108
```
82
109
83
110
In general CVODE_BDF tends to work well on these problems.
0 commit comments