A Preconditioned Iterative Solution Scheme for Nonlinear Parabolic Systems Arising in Air Pollution Modeling

A preconditioned iterative solution method is presented for nonlinear parabolic transport systems. The ingredients are implicit Euler discretization in time and finite element discretization in space, then an outer-inner (outer damped inexact Newton method with inner preconditioned conjugate gradient) iteration, further, as a main part, preconditioning via an `-tuple of independent elliptic operators. Numerical results show that the suggested method works properly for a test problem in air pollution modeling.


Introduction
Nonlinear time-dependent reaction-convection-diffusion (transport) systems arise in various situations in applied mathematics and mathematical modeling, often leading to large-scale, computationally challenging problems. Parabolic transport systems are often encountered in the form x, u 1 , . . . , u ) = g i , i = 1, . . . , , (1.1) which contains nonlinear coupling in the reaction terms. Such problems frequently arise in environmental modeling, for instance in the study of the transport of air pollutants, where u i are concentrations of chemical species. These systems may consist of a huge number of equations: for instance, in [14] a model with more than 70 chemical reaction is analysed in details, but larger models can also occur. A further specific feature of air pollution models is that chemical reactions are described by very stiff systems of ODEs, which makes the system ill-conditioned. There are various techniques for the numerical handling of such types of problems, see e.g. the monographs [9,13,14,15]. Due to the compound structure of these systems, a widespread approach is to handle certain parts of the problem separately. For example, operator splitting is often used to reduce the original problem to successive simpler ones, or implicit-explicit (IMEX) time stepping is widely applied in order to enable to treat the less stiff part of the problem explicitly. The cost of such simplifications can be certain extra error, as is mostly the case with operator splitting, or a reduced stability, which arises for IMEX methods compared to fully implicit time stepping.
Any of the usual solution methods falls into one of two types: one discretizes either first in time or first in space. In our paper we discretize first in time. After time-discretization a coupled elliptic system has then to be solved on each time level. Since the system of ODEs that stems from the chemical part is generally stiff, the usage of implicit time-discretization schemes seems inevitable, in fact we apply fully implicit time stepping for greater stability. Hence a nonlinearity appears in the elliptic problems. Following [1], we propose an outer-inner (outer damped inexact Newton method with inner preconditioned conjugate gradient) iteration for solving the finite element discretization of the nonlinear elliptic problems. The main part of this method is preconditioning using the discretization of an -tuple of independent elliptic operators as preconditioner. This implies that the preconditioning matrix has a block-diagonal structure, and the auxiliary problems can be solved with a cost proportional to that of a single PDE, in contrast to solving the linearized PDE systems. Here we note that there are various preconditioning strategies for such problems, some of them similar to ours, for instance, approximate matrix factorization or block Gauss-Seidel iterations. Our approach uses equivalent operators following [2], which results in block diagonal preconditioners and lead to mesh independent convergence rates for the conjugate gradient method.
We note that the applicability of the proposed method is more general than the mentioned air pollution models. In fact, this approach can be used in a similar way for other time-dependent parabolic reaction-convection-diffusion systems arising in the context of chemical or biological interactions. For instance, models describing chemotaxis-growth systems, when nonlinear advection part of the model starts to be a leading mechanism of dynamics, are described in [3,12], and various other reaction-convection-diffusion systems are given in [5].
The main goal of the paper is the application and numerical illustration of the above-mentioned theoretical results to a real-life parabolic transport system in air pollution modeling. We consider a simplified two dimensional air-pollution model described in [7], where the nonlinear chemical part is taken from [8]. The results of [1] are adapted, involving the time-step parameter, to the solution of the elliptic systems arising after the implicit time-discretization. The numerical results show that the suggested method performs satisfactorily for this realistic model problem.

A two-dimensional model problem
In [7] the following two-dimensional environmental model is proposed: x, y, u 1 , . . . , u ) = 0, i = 1, . . . , , (2.1) x, y) are concentrations of chemical species (pollutants) and the diffusion coefficient K is some positive constant. This system describes the horizontal submodule of the process. The assumption of constant diffusion coefficient K := K x = K y is not a severe restriction, since (in contrast to the vertical coefficient) the horizontal diffusion coefficients are often chosen in this way. It is also assumed that Dirichlet boundary conditions are imposed, where the boundary values are usually interpolated from numerical results of a large scale global model. Following [7], the advection part in Eq. (2.1) has the form The spatial domain is a square with side length 500 km (a 1 = a 2 = 0 =: a, b 1 = b 2 = 500 =: b) and the time span is a whole day. The time is measured in minutes, hence the length of the time interval is 1440 min. This problem can be considered as a special case of the Unified Danish Eulerian Model (cf. [4,14]), and is proposed in [7] as a proper model for testing numerical methods for air pollution problems. In particular, the original number of species = 56 is simplified to = 10 using the most characteristic pollutants in the process, and an explicit circular wind field is involved.

The chemical part of the model
The main part of the model Eq. (2.1) comes from chemistry. It consists of a source term for each pollutant and a nonlinear reaction term: This part causes both the nonlinearity and the coupling. The emission E i may depend on t, but not all chemical species are emitted, some of them are created by the chemical reactions during the transport in the atmosphere. For this experiment we have chosen zero emission for all pollutants (a so-called puff), hence the reaction part Q i equals to R i . The chemistry of the simplified model is defined in such a way that the number of chemical components can be kept on a minimal level. It is noted in [8] that the obtained chemical module is not claimed to reconstruct the complete real process, i.e. the computed concentrations do not necessarily reflect factual values of chemical species in the atmosphere. As mentioned before, this simplified model reflects the essential behaviour and can serve as a proper test problem for numerical computations. The chemical scheme used in the model can be found in [8, Tab. 1], but we repeat them in Table 1 for the sake of completeness.

FEM discretization and iterative solution
First we apply implicit time discretization to Eq. (2.1), and the nonhomogeneous boundary condition can be transformed into a homogeneous one by Hence we now assume that the boundary condition is homogeneous. Then our system has the form on each time level, which is a nonlinear elliptic system. Its more concrete form will be recalled later in (4.2), including the form of the functions f i ; the presence of the factor 1/τ implies the coercivity of the problem for small enough τ and thus the uniform ellipticity of the linearized systems. For brevity, we write (3.1) as using the vector notation. For any u ∈ H 1 0 (Ω) let We apply the finite element method (FEM) for the discretization of (3.3) using the N -dimensional FEM subspace V h = span{ϕ 1 , . . . , ϕ N } ⊂ H 1 0 (Ω) and seek the FEM solution u h ∈ V h satisfying We note that such a standard Galerkin method will prove to be appropriate, since the diffusion term is not negligible in the model, therefore no oscillations will appear in the numerical tests even without using some stabilization technique.

Outer iteration: Newton's method
The operator F h : V h → V h and the function g h ∈ V h are defined by the identities via the Riesz representation theorem, thus the problem can be written as a nonlinear algebraic system We apply the damped inexact Newton method (DIN) for the iterative solution of problem (3.4). The construction of the DIN method and the related convergence results are well-known, for completeness we briefly summarize them as follows.
Algorithm 1 [DIN]. Let u 0 ∈ V h be arbitrary. The sequence (u n ) ⊂ V h is constructed as follows: • Denoting the residual by Under suitable smoothness, growth and coercivity conditions the following theorem holds, for the proof see [1], or [6, Thm. 5.12] for a more general case.
with some 0 < γ ≤ 1, then the convergence is locally of order 1+γ, that is the convergence is linear for n 0 steps 2L (here and in the definition of σ n the constant L comes from the Lipschitz continuity of F ), and further on (as σ n ≡ 1) u n − u h H 1 0 ≤ d 1 q (1+γ) n−n 0 with some d 1 > 0, 0 < q < 1.
We note that the order 1+γ of convergence depends on the order of accuracy in the inner solution, i.e. on the estimate δ n ≤ const · F h (u n ) − g h γ H 1 0 that appears in the condition of the theorem. Here γ = 1 is the optimal choice since 1 + γ can equal at most 2.

Inner iteration: preconditioned CG for the linearized problems
In each step the construction of u n requires the approximate solution of the linearized problem which is equivalent to the FEM solution in V h of the linear elliptic system The theory of equivalent operators (cf. [2]) can be applied to the auxiliary linear problem (3.6) which can be solved by a conjugate gradient (CG) method using a suitable preconditioner. Let η i ∈ L ∞ (Ω), η i ≥ 0 be suitable functions and for p i ∂Ω = 0, and define the -tuple of independent elliptic operators The terms η i p i approximate the lower order term in the ith component of the system: since the system comes from a discretized parabolic problem, a natural choice is the constant coefficient η i ≡ 1/τ , see in the tests in Subsection 4.2. The preconditioner for the discrete system (3.7) is defined as the stiffness matrix S h of S in H 1 0 (Ω) . Then we apply the CGN algorithm (see, e.g., [2]) for the preconditioned nonsymmetric system Since the operators S i are decoupled, in each Newton step the linearized system (3.6) is preconditioned by an -tuple of independent symmetric elliptic operators, which means that the preconditioning matrix S h has a blockdiagonal structure. This enables parallel computation of the solution of the auxiliary problems in the CGN, which was demonstrated for a linear elliptic test system in [10]. Moreover, combining the convergence results for the CGN (cf. [2,11]) and the DIN Algorithm 1, the combined iteration provides mesh independent convergence, with superlinear convergence rate for both the inner and outer iterations.

Numerical treatment of the model problem
where Lu := −K ∆u + b · ∇u stands for the linear elliptic part of the problem.

Time discretization
We assume as before that the nonhomogeneous boundary condition is already eliminated, i.e. γ = 0. Using the backward Euler method for time discretization, for any j ∈ N and steplength τ > 0 we get which is a nonlinear elliptic problem having the form (3.2) with u j+1 The presence of the factor 1/τ implies the coercivity of the problem for small enough τ and thus the uniform ellipticity of the linearized systems; moreover, the lower order part is thus dominating in the problem if τ is small. In the numerical experiments, based on [7,8], the spatial domain is the square Ω = [0, 500] 2 with side length 500 km, the length of the time interval [T 0 , T 1 ] is 1440 min, the time discretization parameter τ = 5 min and the number of equations is = 10. The initial conditions (on time level j = 0) are the constant functions u 0 = 10 3 10 3 10 3 5 · 10 3 5 · 10 3 10 2 10 −2 10 −2 10 −3 10 −11 , measured in mol/km 3 , and the boundary functions are chosen to be periodic: γ i has the form γ i (t) = const i ·(sin(t/C) + 2), where C is a constant and the constants const i are chosen in such a way that the compatibility of the boundary and initial data is ensured. We use a daily average value for the solar zenith angle θ. Finally, the diffusion coefficient is set to be K = 1.8 km 2 /min (as suggested in [14]).

Numerical experiments and results
Now we solve system (4.1) with = 10 equations, with initial and boundary conditions described above, using the implicit Euler method with time discretization parameter τ = 5. On each time level the obtained nonlinear elliptic problem (4.2) is discretized by using an equidistant, piecewise linear finite element grid with spatial mesh parameter h with N h = 500, where N is the number of subintervals in both spatial directions.
We use Algorithm 1 as outer iteration. Since the Lipschitz constant L is unknown, an adaptive strategy is used for choosing σ n and δ n . The linearized problem in each Newton-step is solved by the CGN algorithm, where the preconditioner comes from the stiffness matrix of (3.9) which consists of decoupled, hence independent symmetric Helmholtz-type preconditioners S i p i := −K ∆p i + 1 τ p i (i = 1, . . . , 10).
Here the chosen preconditioner does not depend on the time level, the same matrix S h is used throughout the whole time interval. In each time level we used r h,n H 1 0 ≤ 10 −3 r h,0 H 1 0 as stopping criterion for the Newton method. Inside each Newton-step the auxiliary linear equations (3.5) were solved up to the relative precision δ n .
Numerical experiments were carried out in the cases N = 4, 8, 16, 32, where N denotes the number of subintervals in each spatial directions. In Table 4 the spatial convergence is illustrated, where the numbers reflect a kind of worst case scenario. This means that -considering all the four mesh parameters -the maximum number of Newton iterations was 10 during the time integration. We have calculated on each time level jτ the values of the relative errors Q  The ith row in Table 4 shows the maximal values of the quotients (4.3) and the maximum number of PCGs needed in the ith DIN step. When the value of Q reaches the relative tolerance 10 −3 , the Newton iteration terminates. Fig. 1 shows the number of Newton iterations for all mesh parameters as the time elapses. The apparently periodically fluctuating behaviour can be explained by the imposed time-periodic boundary condition.
In each fixed time level jτ and in each Newton iteration step we solve the auxiliary linear problem (3.5) approximately by a preconditioned CG method, until the relative error is reached, where δ n depends on the current Newton steplength σ n and has an approximate magnitude 10 −2 to 10 −4 . In Table 5  PCG iterations for all the four meshes. As it can be seen from Table 5, the finest mesh has required the most PCG iterations, although less Newton steps were needed in that case. This behaviour can also be seen in Table 4. Some experiments were run with different time stepsizes τ = 20, 10, 5, 2.5, 1.25 for the meshes with parameter h = 1/16. The differences of the solutions (in the common points of the time interval) are shown in Fig. 2, when the time parameter is halved. It can be seen that these differences tends to zero as τ → 0, which shows the convergence of the proposed method in time. The speed with which the graphs of these difference functions tend to zero supports numerically the expected first order convergence in time. It can also be noticed that the error can be quite large, hence the choice of small values of τ is recommended for the time integration.  Fig. 3 the distribution of some characteristic pollutants are shown at the end of the time interval. It can be seen that the constant initial values have   been left relatively intact in the middle of the domain, but they have been stretched near the boundary by the sinusoidal boundary condition.
As a conclusion of the numerical experiments, we may observe that the suggested method performs satisfactorily for this realistic model problem, in particular, reasonable convergence has been achieved without the need of some splitting of the reaction-convection-diffusion operator.