Chemistry Solvers in CMAQ

1

• Chemistry Solvers in CMAQ
• 1 SMVGEAR
• 2 QSSA
• 3 MEBI

Daewon W. Byun
UH Air Quality Modeling and Monitoring Center
2
Chemical Reaction Kinetics

• The rates of chemical reaction determine whether
  a species is formed or destroyed by gas-phase
  chemistry.
• AQ gas mechanisms treat all reactions as if they
  are elementary
• We use the laws of reaction kinetics directly to
  develop mathematical expressions.
• We will study
• (1) the rate expressions
• (2) the forms of the rate constants

3
Chemical Reaction Rates
4
Reactions involving third body
Several important ter-molecular reactions involve
O2 and/or N2 which mediate those reactions by
absorbing energy from exothermic bi-molecular
reactions. The third body M N2 O2. Since
their concentrations are relatively stable,
mechanism developers convert second- or
third-order reactions that include these species
to a reaction one order lower by multiplying the
higher-order rate coefficient by the
concentration of M, O2, or N2.
5
Rate Constant Formula
6
Rate Constant Formula
7
Reaction Equations
The system of non-linear, ordinary differential
equations (ODEs) for N species can be expressed
as follows
(8-14a)
with the initial conditions
(8-14b)
Nonlinear Stiff equation
8
Stiff Nonlinear Equations

• It is stiff because of the widely varying time
  scales of the chemical reactions and complex
  interactions among species.
• A stiff system can be described mathematically as
  one in which all the eigenvalues of the Jacobian
  matrix of the system are negative, and the ratio
  of the absolute values of the largest-tosmallest
  real parts of the eigenvalues is much greater
  than one.

• For atmospheric chemistry problems, the ratio is
  often greater than 1010, making the system very
  stiff

9
Solver 1 SMVGEAR
Sparse-Matrix, Vectorized Gear (M. Jacobson
Turco, 1994). .
Based on a numerical solver with the algorithm
developed by Gear (1971b). The technique is
implicit in method, does not amplify errors from
one step to another and incorporates automatic
step size and error control. The solver is often
been used to evaluate other faster solution
methods for accuracy.
10
GEAR Solver
The method is implicit since concentrations at
the desired time step n depend on values of the
first derivatives contained in f(cn,tn) that are
functions of the concentrations at the same time.
The order of the method corresponds to the
number of concentrations at previous time steps
that are incorporated in the summation on the
right hand side of Equation 8-15.
11
GEAR Solver
12
GEAR Solver
13
GEAR Solver
14
SMVGEAR Solver

• Jacobson and Turco (1994) have modified the Gear
  algorithm to incorporate additional computational
  efficiencies that can achieve speedups on the
  order of 100 on vector computers.
• About half of the improvement is attributed to
  enhanced vectorization, and half to improved
  matrix operations.
• Because of the improved matrix operations,
  SMVGEAR also runs faster than traditional Gear
  solvers on non-vector machines, but the greatest
  benefit will be obtained with vector machines.
• In SMVGEAR, the modeling domain is divided into
  blocks of cells, and the Gear algorithm is
  applied to each cell within a block
  simultaneously. With this structure, the length
  of the innermost loops for most calculations is
  equal to the number of cells in a block (500
  cells).

15
SMVGEAR Solver

• Much of the computational intensity associated
  with the Gear method arises from the matrix
  operations that are needed to perform the Newton
  iterations in the corrector step. Two techniques
  
• First, all known cases of multiplication by zero
  in matrix multiplication and decomposition are
  eliminated. This is particularly beneficial for
  atmospheric chemistry problems since the Jacobian
  matrices are almost always sparse (i.e., they
  contain a large number of zero entries).

16
SMVGEAR Solver

• However, decomposition techniques that are
  applied to these matrices often result in
  substantial fill-in, thereby reducing the
  benefits of employing sparse-matrix techniques.
• To maintain maximum sparsity after the
  decomposition operation, SMVGEAR orders the
  species in the Jacobian such that those with the
  fewest partial derivative terms are located in
  the top rows, and those with the most are in the
  bottom rows.
• At the very beginning of the program, the
  ordering is done and a symbolic decomposition is
  performed to identify multiplies by zero.
• Since the computations associated with the matrix
  operations are determined entirely by the
  structure of the chemical mechanism, do this only
  once and the results can then be applied to every
  cell uniformly.

17
QSSA Solver
Quasi-Steady-State Approximation. The QSSA
solver is a low order, explicit solver that
exhibits good stability for stiff systems.
Although less accurate than the Gear solver, it
is still a reasonably accurate, fast solver that
is especially suitable for large scale grid
models. The solver developed for the CMAQ is a
predictor/corrector version based on the one
developed by Lamb and used in the Regional
Oxidant Model (Lamb, 1983, and Young et al.,
1993).
18
QSSA Solver
19
QSSA Solver
20
QSSA Solver
In the first stage, an optimal chemistry time
step interval Dt is determined at time tn based
on a tolerance parameter, l , such that
The algorithm attempts to weight the time step
determination against species whose
concentrations are very small compared to the
primary oxidants and some of the fastest reacting
radicals.
21
QSSA Solver
22
QSSA Solver
23
MEBI Solver
The MEBI (Modified Euler Backward Iterative)
solver developed by Hertel et al. 1993 uses
functional iteration to obtain a solution to the
implicit Euler backward approximation. To
speed up convergence to a solution, however,
groups of species that are strongly coupled are
isolated from the rest of the species in the
mechanism and analytical solutions are derived
for their concentrations. The analytical
results are then used with the backward Euler
approximation to obtain estimates of the
concentrations for the remainder of the mechanism
species. Estimates of the group and Euler
backward species concentrations are constantly
updated using the results from the previous
iteration until convergence is achieved.
24
MEBI Solver
The accuracy of the solution obtained with the
MEBI is controlled by the size of the chemistry
time step as well as by convergence tolerances
for Newton-Raphson iterations and convergence
tolerances for each species for the Euler
iterations. The tolerances used in the CMAQ
system were selected to provide relatively
accurate solutions without greatly sacrificing
computational efficiency. The time step has a
default value of 2.5 minutes, and convergence
tolerances are set for each species individually
on the basis of the work done by Huang and Chang
2000 and subsequent testing and comparison with
other solvers.  
25
MEBI Solver

• Hertel et al. 1993 applied the technique to a
  version of the CB4 mechanism and found it to be
  more efficient than a version of the QSSA that
  employed species lumping.
• In the paper, four different groups of strongly
  coupled species were defined
• NO, NO2, O3, and O(3P)
• OH, HO2, HONO, and HNO4 and
• NO3 and N2O5, and
• PAN and C2O3.
• Analytical solutions were obtained from the
  production and loss rates for each species in the
  group independently of the other species in the
  mechanism.

26
MEBI Solver
Since the species in the first two groups
interact with a number of other species in the
mechanism and they interact with each other, the
analytical solutions can become quite complex.
Instead of relying on analytical solutions for
these two groups, Huang and Chang 2000 used
Newton-Rapson iteration to obtain numerical
solutions. Although this approach requires
matrix inversions, they are not particularly time
consuming since the size of the matrices is small
(4 by 4). Analytical solutions are retained for
the last two groups.
27
MEBI Solver
The CMAQ implementation of the MEBI solver is
based on the approach of Huang and Chang 2000
as applied to the CMAQ version of the
CB4/RADM/SAPRC99 mechanisms linked to aerosols
and aqueous chemistry. With this approach, one
step towards generalization has been taken since
the Euler backward iterations and the numerical
solutions to the first two groups of coupled
species can be generalized. It would be possible
to completely generalize the method by obtaining
numerical solutions for the other two groups as
well. The MEBI solver has been implemented in
CMAQ in a non-generalized fashion to avoid
potential misapplication to other mechanisms.
The feasibility of using this approach with the
other CMAQ mechanisms will be investigated in the
future.   One additional modification was made
to the version of the MEBI solver incorporated in
CMAQ.
28
MEBI Solver
Preliminary testing of the original solver
revealed that, in some rare instances, the
backward Euler iteration did not converge within
the allotted number of iterations. The
convergence problem was associated with the
negative stoichiometric coefficients that are
used for a few products in the CB4 mechanism.
When initial concentrations and production rates
of these species were very low, the negative
product coefficients can cause the concentrations
of these species to go negative. In the
standard backward Euler iterations, the
concentrations of these species were constantly
decreasing, but at a very slow rate. The Euler
iterative algorithm was modified to set the
concentration of these species to zero whenever
the loss rate exceeded the initial concentration
plus the amount produced over the time step.
This resulted in significant increases in the
convergence rate in these instances. Since these
situations happened very rarely, however, the
change had virtually no effect on over all
computer time.