A local character based method for solving linear systems of radiation diffusion problems

1 Introduction

In the simulation of astrophysics and inertial confinement fusion (ICF), it is needed to solve radiation diffusion equations. This is a strongly nonlinear system, and the discretized equations are solved by some nonlinear iterative methods, such as Newton kind of methods, Picard method, and source iteration method for multi-group radiation diffusion equations [2, 13, 17, 25]. In each nonlinear iteration, some linearized equations should be solved. In numerical simulations of ICF, usually $O (10^{4})$ time steps are needed to finish the whole simulation, and a series of linear equations $A x = b$ should be solved. The solution time for linear equations dominates the real numerical simulation cost [1, 23, 25]. It is important to design an efficient method for solving linear equations in ICF and astrophysics simulations.

Usually, the preconditioned Krylov methods [19] are used to solve the linear systems in ICF simulations. For a preconditioned Krylov method, there are four basic procedures: the construction and application of a preconditioner, the choice of an initial iterate, the computation of the next iterate and the check for stopping criterion [9, 22]. Each procedure has its own effect on the iterative method. In past years, much research work has been focused on preconditioning. Among them, AMG preconditioner is the most widely used one in the simulation of radiation diffusion equations [2, 17, 23, 25].

For a series of linear systems, several techniques for choosing an efficient initial iterate have been proposed in recent years. For linear systems with the same coefficient matrix and different right-hand sides, a method for constructing an efficient initial iterate was given in [5]. Later, Chan further considered the choice of initial iterate for a series of equations with different coefficient matrices in [4]. Fischer proposed another method, in which an approximate solution based upon previous solutions was used [7]. In [22], a model reduction method based on proper orthogonal decomposition (POD) was used to guess a better initial iterate. The projection methods were also used in [3, 6, 10, 14]. In [9, 20], a better initial guess was given by the extrapolation of previous solutions. An et al. proposed a kind of method for choosing a nonlinear initial iterate for solving two dimension three-temperature (2D 3T) heat conduction equations in [1]. By solving the nonlinear system on a sub-domain where the temperature varied greatly, the method in [1] can improve the computational efficiency effectively. In [8, 12, 24], some subsystems are constructed to accelerate the convergence of nonlinear problems.

In the simulation of ICF, with time step advancing, the physical variables, such as density, pressure, temperature, etc., vary greatly in some local domain and mildly in the rest domain. This is reflected in the radiation diffusion equations, where the unknown is the temperature or radiation energy. From the perspective of solving linear systems from one nonlinear iteration to the next iteration, or from one time step to the next, the solution only varies greatly in some local domain while mildly in the rest. In ICF simulations, the initial iterate for solving linear equations is usually set as the previous nonlinear iterate or previous time step solution. Therefore, the solution on some local domain where the solution varies greatly has a strong influence on the convergence of an iterative method, while the solution on the rest domain has little influence on the method.

In this paper, we propose a local character-based method for solving the linear systems arising from the simulation of radiation diffusion equations. In this method, a local domain is first constructed, corresponding to a subset of the global degrees-of-freedom, on which the solution is expected to vary greatly. Then, the subsystem on the local domain is solved. At last, the whole system is solved by using the solution on the local domain. Two methods are given to construct the local domain. Specifically, the contributions of this paper are summarized as follows:

A local character-based method for solving linear systems of radiation diffusion equations is proposed.
Two methods are given to construct the local domain. One is based on the spatial gradient and the other is based on the residual.
Some numerical experiments are conducted to verify the effectiveness of the proposed method. Compared with the AMG preconditioned Krylov method, the local character-based method can reduce the solution cost in most cases, and 40% cost can be saved at most.

The rest of the paper is organized as follows. The detail of the local character-based method is presented in section 2, and also two methods for constructing the local domain are presented in this section. In section 3, some numerical results are given and also the choice of the parameter for the local character-based method is analyzed in this section. Finally, section 4 is the summary and some remarks of the paper.

2 The local character-based method

In this section, the local character-based method is introduced. The basic algorithm is presented in section 2.1. Two methods for constructing the local domain are introduced in the rest. The first one, which is based on the spatial gradient, will be presented in section 2.2. The second one, which is based on the residual, will be presented in section 2.3.

2.1 The local character-based algorithm

Consider a linear system resulting from implicit time propagation of discretized time dependent partial differential equations (PDEs), or the implicit linearization of a nonlinear system of PDEs:

A x = b, A \in R^{N \times N}, x, b \in R^{N},

where

A = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} a_{11} & a_{12} & \dots & a_{1 N} a_{21} & a_{22} & \dots & a_{2 N} ⋮ & ⋮ & ⋱ & ⋮ a_{N 1} & a_{N 2} & \dots & a_{N N} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠, x = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} x_{1} x_{2} ⋮ x_{N} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎠, b = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} b_{1} b_{2} ⋮ b_{N} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ .

Eq (2.1) is usually obtained from some real applications, such as the application of radiation diffusion problem, etc. For solving Eq (2.1), an iterative method, particularly a Krylov method is often used. Because the linear system is obtained in the nonlinear iteration on some time steps, an effective initial iterate $x^{(0)}$ , which is the solution of the last nonlinear iterate or the solution at last time step, is used. In the simulation process of ICF, the temperature and radiation energy only varies greatly in some local domain and varies mildly in the rest domain. Therefore, if an iterative method is used to solve Eq (2.1), only some components of the iterate will change greatly from the initial iterate $x^{(0)}$ to the converged solution.

For purpose of clarity, in the following we introduce some notations. Let $Ω = {1, 2, \dots, N}$ be the index set of all components (degrees-of-freedom) of the linear system (2.1). The linear system considered in this paper is discretized from a partial differential equations, each component corresponds to one spatial point (for example, a node of the grid) in the computing domain of PDEs. Therefore, a component will also be called a point in the following discussion. $Ω$ can also be regarded as an algebraic “domain”. Assume that, from the initial iterate $x^{(0)}$ to the converged iterate, the solution varies greatly at $K$ components. For convenience, the index set of these $K$ components is denoted as $Ω\scriptsize local={i1,i2,…,iK}$ . Each component in $Ω\scriptsize local$ is also called a bad point. $Ω\scriptsize local$ is called a local domain in the following. Let $¯Ω\scriptsize local=Ω∖Ω% \scriptsize local$ . Based on the partition

Ω=Ω\scriptsize local∪¯Ω% \scriptsize local,

the linear system (2.1) can be partitioned as

(\begin{matrix} B & E F & C \end{matrix}) (\begin{matrix} x_{B} x_{C} \end{matrix}) = (\begin{matrix} b_{B} b_{C} \end{matrix}),

(1)

where

x_{B} = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} x_{i_{1}} x_{i_{2}} ⋮ x_{i_{K}} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠, b_{B} = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} b_{i_{1}} b_{i_{2}} ⋮ b_{i_{K}} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠,

and

(\begin{matrix} B & E F & C \end{matrix}) = Q^{T} A Q, (\begin{matrix} x_{B} x_{C} \end{matrix}) = Q^{T} x, (\begin{matrix} b_{B} b_{C} \end{matrix}) = Q^{T} b .

Here $Q \in R^{n \times n}$ is a permutation matrix obtained by exchanging the rows $l$ and $i_{l}$ ( $l = 1, 2, \dots, K$

) of the identity matrix

I \in R^{n \times n}

In Eq (1), $x_{B}$ is expected to change greatly from the initial iterate to the converged solution, while $x_{C}$ is expected to change mildly from the initial iterate to the converged solution. In the local character-based method, first the subsystem

B x_{B} = b_{B} - E x_{C}^{(0)}

(2)

will be solved, and the solution ${¯ x}_{B}$ will be obtained. Then

~ x = Q (\begin{matrix} {¯ x}_{B} x_{C}^{(0)} \end{matrix})

(3)

will be used as an initial iterate to solve the whole system (2.1) iteratively. The specific local character-based algorithm is described by Algorithm 1.

1:Input:

A

x^{(0)}

and

b

2:Construct the subset

Ω\scriptsize local

3:Construct the permutation matrix

Q

4:Construct the sub-matrix

B

C

E

and

F

. Partition the right hand side into

b_{B}

and

b_{C}

, and partition the initial iterate into

x_{B}^{(0)}

and

x_{C}^{(0)}

5:Solve the subsystem (2), obtain the solution

{¯ x}_{B}

6:Assemble

~ x

defined in (3) by using

{¯ x}_{B}

and

x_{C}^{(0)}

7:Apply several Gauss-Seidel iterations to Equations (2.1) so that

~ x

is updated.

8:if

~ x

satisfies the convergence criterion then

9: Let

x = ~ x

10:else

11: Solve the whole system (2.1) by using the initial iterate

~ x

, and obtain the solution

x

12:end if

13:Output:

x

Algorithm 1 Local character-based algorithm

In Algorithm 1, there are two key ingredients: the first one is the construction of the subset $Ω\scriptsize local$ in Line 2; the second one is the solution of the local subsystem in Line 5. The construction of $Ω\scriptsize local$ is the base of the algorithm. By $Ω\scriptsize local$ , the whole system can be partitioned into two subsystems, and particularly the local subsystem is defined. By solving the local subsystem, a solution on $Ω\scriptsize local$ is obtained. After assembling the whole solution $~ x$ as the initial iterate, the whole system will be solved by an iterative method. Because the solution on $¯Ω\scriptsize local$ is expected to change mildly, $x_{C}^{(0)}$ is expected to be very near to the final solution. Therefore $~ x$ will be very near to the solution. Actually, in some of our numerical cases, the residual corresponding to $~ x$ is so small that the convergence criteria is satisfied, and $~ x$ is a final solution.

The construction of the local domain, or the subset $Ω\scriptsize local$ , plays a key role in Algorithm 1. In the following, two methods will be introduced to construct the local domain.

2.2 Gradient-based local domain construction

The linear system (2.1) is discretized from a partial differential equations (such as radiation diffusion equations). For the solution of a partial differential equations, its gradient can be defined. For example, if $T^{n} (x, y)$ is a solution of a two dimensional radiation diffusion equation at time step $n$ , then the gradient of $T^{n} (x, y)$ is defined as $\nabla T^{n} (x, y) = (\partial T^{n} / \partial x, \partial T^{n} / \partial y)$ . The norm of the gradient, $∥ \nabla T^{n} (x, y) ∥$ , can be used to predict the variation of the solution from time step $n$ to time step $n + 1$ .

For the linear system (2.1), the similar norm of the gradient can be defined, which can be used to construct the local domain $Ω\scriptsize local$ . For a point $i$ , if the norm of the gradient of the initial iterate is large at this point, then it is likely that the solution may change greatly around this point, which will be considered as an element for the local domain $Ω\scriptsize local$ . Specifically, for a given row $i$ of the matrix $A$ , the non-zero entries of this row reflect the adjacency relationship of the point $i$ to other points. Therefore, the gradient $g_{i}$ at point $i$ (component $i$ ) is defined by

g_{i} = \sum \begin{matrix} 1 \leq j \leq N a_{i j} \neq 0 \end{matrix} ∣ ∣ x_{i}^{(0)} - x_{j}^{(0)} ∣ ∣ .

(4)

Note that this definition is not the exact gradient of a field (such as temperature, density, etc.) in two-dimensional or three-dimensional space, but just a mimicry. More precisely, this is very similar to the $l_{1}$ norm of the gradient at each point $i$ . For example, if $A$ is obtained by the five-point difference scheme, then the summation in (4) includes four terms, as shown in Fig. 1, in which $p$ and $q$ are used as indices for cells. By definition in (4), the gradient of cell $(p, q)$ is

	$g_{p, q}$	$=$	$\| g_{p - \frac{1}{2}, q} \| + \| g_{p + \frac{1}{2}, q} \| + \| g_{p, q - \frac{1}{2}} \| + \| g_{p, q + \frac{1}{2}} \|$
		$=$	$\| x_{p - 1, q} - x_{p, q} \| + \| x_{p + 1, q} - x_{p, q} \| + \| x_{p, q - 1} - x_{p, q} \| + \| x_{p, q + 1} - x_{p, q} \| .$

Figure 1: Gradient for five-point scheme.

Based on the definition of the gradient (4), the local domain can be constructed. For point $i$ , if

g_{i} > α max 1 \leq j \leq N g_{j},

then $i$ will be a member of the local domain index set $Ω\scriptsize local$ . Here $α$ is a prescribed parameter. The specific algorithm is shown in Algorithm 2.

1:Input:

A

x^{(0)}

α

2:Initialize:

g_{max} = 0

Ω\scriptsize local={}

K = 0

3:for

i = 1 : N

4: Compute

g_{i}

by equation (4).

5: if

g_{i} > g_{max}

then

g_{max} = g_{i}

7: end if

8:end for

9:for

i = 1 : N

10: if

g_{i} > α g_{max}

then

11:

Ω\scriptsize local={Ω\scriptsize local,i}

12:

K = K + 1

13: end if

14:end for

15:Output:

Ω\scriptsize local

K

Algorithm 2 Gradient-based local domain construction

2.3 Residual-based local domain construction

The second method for constructing the local domain is based on the residual of the linear system. The method consists of two steps. First find the components with large residual, and the bad points set, or the local domain, $Ω\scriptsize local$ , will be constructed. Second consider the influence of the bad points on other degrees-of-freedom, and the set of bad points then extended to include other degrees-of-freedom on which the solution is also likely to change.

2.3.1 Finding bad points

To iteratively solve the linear system (2.1), a commonly used convergence criteria is given by

∥ b - A x^{(k)} ∥_{2} \leq ϵ ∥ b ∥_{2},

(5)

where $ϵ$ is a prescribed tolerance and $x^{(k)}$ is an iterative solution. Since

∥ b - A x^{(k)} ∥_{2} \leq \sqrt{N} ∥ b - A x^{(k)} ∥_{\infty},

the condition (5) will be satisfied if

∥ b - A x^{(k)} ∥_{\infty} \leq \frac{ϵ ∥ b ∥_{2}}{\sqrt{N}},

max i ∣ ∣ ∣ ∣ b_{i} - N \sum j = 1 a_{i, j} x_{j}^{(k)} ∣ ∣ ∣ ∣ \leq \frac{ϵ ∥ b ∥_{2}}{\sqrt{N}},

which is equivalent to

∣ ∣ ∣ ∣ b_{i} - N \sum j = 1 a_{i, j} x_{j}^{(k)} ∣ ∣ ∣ ∣ \leq \frac{ϵ ∥ b ∥_{2}}{\sqrt{N}}, i = 1, 2, \dots, N .

Now, by the initial iterate $x^{(0)}$ and the convergence criteria tolerance $ϵ$ , a component $i$ is defined as a bad point if

∣ ∣ ∣ ∣ b_{i} - N \sum j = 1 a_{i, j} x_{j}^{(0)} ∣ ∣ ∣ ∣ > \frac{ϵ ∥ b ∥_{2}}{\sqrt{N}} .

Specifically, the initial local domain $Ω\scriptsize local$ is defined by

Ω\scriptsize local={i:i=1,2,…,N;∣∣ ∣∣bi−N∑j=1ai,jx(0)j∣∣ ∣∣>ϵ∥b∥2√N}.

(6)

2.3.2 Expanding the local domain

When the local subsystem of the bad points is solved, the solution at the bad points is changed. This then has further influence on points adjacent to the bad points. Therefore, the set of the bad points will further be extended. That is, the subset of the bad points (local domain) obtained by (6) will further be expanded.

As discussed in Subsection 2.1, assume that ${¯ x}_{B}$ is the solution of the local domain subsystem, and let $~ x$ defined by (3). The $j$ -th residual corresponding to $x^{(0)}$ and $~ x$ are given, respectively, as

r_{j}^{(0)} = b_{j} - N \sum i = 1 a_{j, i} x_{i}^{(0)}, % and {~ r}_{j} = b_{j} - N \sum i = 1 a_{j, i} {~ x}_{i} .

By the residuals, the influence of the local domain to the point $j$ can be given.

Definition 2.1

Assume that $Ω\scriptsize local$ is a local domain. For one point $j∈¯Ω\scriptsize local$ , let $r_{j}^{(0)}$ be the residual at point $j$ with the initial iterate, and let ${~ r}_{j}$ be the residual at point $j$ after solving the local domain subsystem. If $∣ ∣ r_{j}^{(0)} ∣ ∣ \leq \frac{ϵ ∥ b ∥_{2}}{\sqrt{N}}$ and $∣ ∣ {~ r}_{j} ∣ ∣ > \frac{ϵ ∥ b ∥_{2}}{\sqrt{N}}$ , then the local domain $Ω\scriptsize local$ has influence on the point $j$ ; otherwise, the local domain $Ω\scriptsize local$ has no influence on the point $j$ .

It is easy to see that

${~ r}_{j}$	$=$	$b_{j} - N \sum i = 1 a_{j, i} {~ x}_{i}$
	$=$	$b_{j} - K \sum l = 1 a_{j, i_{l}} {~ x}_{i_{l}} - N \sum l = K + 1 a_{j, i_{l}} {~ x}_{i_{l}}$
	$=$	$b_{j} - K \sum l = 1 a_{j, i_{l}} {¯ x}_{l} - N \sum l = K + 1 a_{j, i_{l}} x_{i_{l}}^{(0)}$
	$=$	$b_{j} - K \sum l = 1 a_{j, i_{l}} ({¯ x}_{l} - x_{i_{l}}^{(0)}) - N \sum l = 1 a_{j, i_{l}} x_{i_{l}}^{(0)}$
	$=$	$r_{j}^{0} - K \sum l = 1 a_{j, i_{l}} ({¯ x}_{l} - x_{i_{l}}^{(0)}) .$

Assume that

∣ ∣ {¯ x}_{l} - x_{i_{l}}^{(0)} ∣ ∣ \leq ∣ ∣ x_{i_{l}}^{(0)} ∣ ∣, l = 1, 2, \dots, K,

which means the relative variation of the solution at all the bad points is less than one (This is satisfied in our numerical tests). Then it is easy to see that

∣ ∣ {~ r}_{j} ∣ ∣ \leq ∣ ∣ r_{j}^{(0)} ∣ ∣ + K \sum l = 1 ∣ ∣ a_{j, i_{l}} x_{i_{l}}^{(0)} ∣ ∣, j = 1, 2, \dots, N .

(7)

This inequality shows that if

∣ ∣ r_{j}^{(0)} ∣ ∣ + K \sum l = 1 ∣ ∣ a_{j, i_{l}} x_{i_{l}}^{(0)} ∣ ∣ \leq \frac{ϵ ∥ b ∥_{2}}{\sqrt{N}},

then the local domain will have not influence on point $j$ (by Definition 2.1). Otherwise, the local domain will be expected to have influence on point $j$ . In this case, the local domain will be extended. Specifically, assume that $j∈¯Ω\scriptsize local$ is a neighbor of $Ω\scriptsize local$ . If

∣ ∣ r_{j}^{(0)} ∣ ∣ + K \sum l = 1 ∣ ∣ a_{j, i_{l}} x_{i_{l}}^{(0)} ∣ ∣ > \frac{ϵ ∥ b ∥_{2}}{\sqrt{N}},

then $j$ will be incorporated into the local domain subset $Ω\scriptsize local$ .

Note that the extension process can be implemented several times until the size of the subset $Ω\scriptsize local$ no longer increases. In practice, a few times of expanding is enough. The specific process is given in Algorithm 3.

In Algorithm 3, two steps are used to construct the local domain $Ω\scriptsize local$ . The first step is to construct $Ω\scriptsize local$ by the component residuals, and those components with large residuals are incorporated into $Ω\scriptsize local$ . The second step is to extend $Ω\scriptsize local$ by at most $E_{max}$ iterations. In each expanding iteration, some neighbors of $Ω\scriptsize local$ are further incorporated into $Ω\scriptsize local$ .

1:Input:

N

A

b

x^{(0)}

ϵ

E_{max}

2:Initialize:

Ω\scriptsize local={}

N = {}

K = 0

3:for

i = 1 : N

4: Compute the initial residual

r_{i}^{(0)} = b_{i} - \sum_{j = 1}^{N} a_{i, j} x_{j}^{(0)}

5: if

| r_{i}^{(0)} | > \frac{ϵ ∥ b ∥_{2}}{\sqrt{N}}

then

Ω\scriptsize local={Ω\scriptsize local,i}

K = K + 1

8: end if

9:end for

10:for

m = 1 : E_{max}

11: Let

N

be the neighbors of

Ω\scriptsize local

12: Let

N\scriptsize local={}

K_{0} = 0

13: for all

j \in N

14: if

∣ ∣ r_{j}^{(0)} ∣ ∣ + \sum_{l = 1}^{K} ∣ ∣ a_{j, i_{l}} x_{i_{l}}^{(0)} ∣ ∣ > \frac{ϵ ∥ b ∥_{2}}{\sqrt{N}}

then

15:

N\scriptsize local={N\scriptsize local,j}

16:

K_{0} = K_{0} + 1

17: end if

18: end for

19: if

K_{0} > 0

then

20: .

21:

K = K + K_{0}

22: else

23: break;

24: end if

25:end for

26:Output:

Ω\scriptsize local

K

Algorithm 3 Residual based local domain construction

Example 1

The following example shows the specific implementation of Algorithm 3. Consider a $9 \times 9$ linear system with the coefficient matrix given by

A = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} 1 & - \frac{1}{2} - \frac{1}{2} & 1 & - \frac{1}{3} - \frac{1}{3} & 1 & - \frac{1}{4} - \frac{1}{4} & 1 & - \frac{1}{5} - \frac{1}{5} & 1 & - \frac{1}{6} - \frac{1}{6} & 1 & - \frac{1}{7} - \frac{1}{7} & 1 & - \frac{1}{8} - \frac{1}{8} & 1 & - \frac{1}{9} - \frac{1}{9} & 1 \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

The solution and the initial iterate are set, respectively, as

x = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} 1 e - 1 1 e - 2 1 e - 3 1 e - 4 1 e - 5 1 e - 6 1 e - 7 1 e - 8 1 e - 9 \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠, x^{(0)} = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} 1 1 e - 1 1.001 e - 3 1.001 e - 4 1.001 e - 5 1.001 e - 6 1.001 e - 7 1.001 e - 8 1.001 e - 9 \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .

The right hand side is set to $b = A x$ .

Consider the application of Algorithm 3 to this linear system. The convergence criterion is set to $ϵ = 1 e - 5$ , and the parameter for maximal number of extension is set to $E_{max} = 6$ . It is easy to check that $ϵ ∥ b ∥_{2} / \sqrt{N} = 3.44 e - 7$ . The initial residual is

r^{(0)} = b - A x^{(0)} = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} - 8.55 e - 1 3.66 e - 1 3.00 e - 2 1.52 e - 7 1.02 e - 8 6.81 e - 10 4.41 e - 11 2.61 e - 12 1.11 e - 13 \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

In the first step of Algorithm 3, points 1, 2 and 3 are marked as the bad points, i.e., $Ω_{l o c a l} = {1, 2, 3}$ . In the second step of Algorithm 3, $Ω_{l o c a l}$ is further extended. The loop of extension process is implemented 4 times. After each loop, the main results are listed in the following table. In the fourth iteration, $N = {7}$ and $| r_{7}^{(0)} | + | a_{76} \cdot x_{6}^{(0)} | = 1.430 e - 7 < ϵ ∥ b ∥_{2} / \sqrt{N} = 3.44 e - 7$ , therefore, $K_{0} = 0$ and no point will be incorporated into $Ω\scriptsize local$ . The loop will break.

$m$	$N$	$N\scriptsize local$	$Ω\scriptsize local$	$K$	$\| r_{j}^{(0)} \| + \sum_{l = 1}^{K} \| a_{j, i_{l}} x_{i_{l}}^{(0)} \|$
1	${4}$	${4}$	${1, 2, 3, 4}$	4	$\| r_{4}^{(0)} \| + \| a_{43} \cdot x_{3}^{(0)} \| = 2.504 e - 4$
2	${5}$	${5}$	${1, 2, 3, 4, 5}$	5	$\| r_{5}^{(0)} \| + \| a_{54} \cdot x_{4}^{(0)} \| = 2.003 e - 5$
3	${6}$	${6}$	${1, 2, 3, 4, 5, 6}$	6	$\| r_{6}^{(0)} \| + \| a_{65} \cdot x_{5}^{(0)} \| = 1.669 e - 6$
4	${7}$	$\emptyset$	${1, 2, 3, 4, 5, 6}$	6	$\| r_{7}^{(0)} \| + \| a_{76} \cdot x_{6}^{(0)} \| = 1.430 e - 7$

Table 1: The extension process of

Ω\scriptsize local

2.4 Parallelization of the local character-based methods

For parallel computing, assume that $p$ processors are used. For simplicity, assume that $N = ¯ N \times p$ , where $¯ N$ is an integer. The matrix $A$ is partitioned by rows, that is

A = {[A_{1 -}^{T} A_{2 -}^{T} \dots A_{p -}^{T}]}^{T},

where $A_{l -}$ is a $¯ N \times N$ matrix, $l = 1, 2, \dots, p$ . The solution and right-hand side are partitioned correspondingly

x = {[x_{1 -}^{T} x_{2 -}^{T} \dots x_{p -}^{T}]}^{T}, b = {[b_{1 -}^{T} b_{2 -}^{T} \dots b_{p -}^{T}]}^{T},

where $x_{l -}$ and $b_{l -}$ are $¯ N$ vectors, $l = 1, 2, \dots, p$ . For convenience, let

I_{l -} = {(l - 1) ¯ N, \dots, l ¯ N}

be the component index set of the $l$ -th processor, $l = 1, 2, \dots, p$ .

The key work in the parallel version of the algorithm is to construct the local domain $Ω\scriptsize local$ . When the local domain is constructed, the local domain subsystem can be obtained easily. For solving the local domain subsystem, any parallel solution method can be used. In the following, Algorithm 4 and Algorithm 5 are given for the gradient based and residual based local domain construction in parallel computing case.

1:Input:

A_{l -}

x_{l -}^{(0)}

I_{l -}

α

2:Initialize:

g_{max} = 0

Ω\scriptsize local={}

K = 0

3:for

i \in I_{l -}

4: Communicate all adjacency components of

x_{i}^{(0)}

from other processors.

5: Compute

g_{i}

by equation (4).

6: if

g_{i} > g_{max}

then

g_{max} = g_{i}

8: end if

9:end for

10:Reduction

g_{max}

by max.

11:for

i \in I_{l -}

12: if

g_{i} > α g_{max}

then

13:

Ω\scriptsize local={Ω\scriptsize local,i}

14:

K = K + 1

15: end if

16:end for

17:Reduction

Ω\scriptsize local

and

K

by sum.

18:Output:

Ω\scriptsize local

K

Algorithm 4 Gradient-based local domain construction in parallel

Compared to Algorithm 2, some communications should be implemented for parallel case, which is shown in Line 4, 10, and 17 in Algorithm 4.

1:Input:

N

¯ N

A_{l -}

b_{l -}

x_{l -}^{(0)}

I_{l -}

ϵ

E_{max}

2:Initialize:

Ω\scriptsize local={}

N = {}

K = 0

3:for

i \in I_{l -}

4: Communicate all adjacency components of

x_{i}^{(0)}

from other processors.

5: Compute the initial residual

r_{i}^{(0)} = b_{i} - \sum_{j = 1}^{N} a_{i, j} x_{j}^{(0)}

6: if

| r_{i}^{(0)} | > \frac{ϵ ∥ b ∥_{2}}{\sqrt{N}}

then

Ω\scriptsize local={Ω\scriptsize local,i}

K = K + 1

9: end if

10:end for

11:Reduction

Ω\scriptsize local

and

K

by sum.

12:for

m = 1 : E_{max}

13: Let

N={the neighbors of Ω\scriptsize local}∩Il–

14: Let

N\scriptsize local={}

K_{0} = 0

15: for all

j \in N

16: Communicate all adjacency components of

x_{j}^{(0)}

from other processors.

17: if

∣ ∣ r_{j}^{(0)} ∣ ∣ + \sum_{l = 1}^{K} ∣ ∣ a_{j, i_{l}} x_{i_{l}}^{(0)} ∣ ∣ > \frac{ϵ ∥ b ∥_{2}}{\sqrt{N}}

then

18:

N\scriptsize local={N\scriptsize local,j}

19:

K_{0} = K_{0} + 1

20: end if

21: end for

22: if

K_{0} > 0

then

23: .

24:

K = K + K_{0}

25: Reduction

Ω\scriptsize local

and

K

by sum.

26: else

27: break;

28: end if

29:end for

30:Output:

Ω\scriptsize local

K

Algorithm 5 Residual based local domain construction in parallel

Compared to Algorithm 3, it is needed to have some communication in Line 4, 11, 16, and 25 in Algorithm 5.

It should be noted that for both Algorithm 4 and Algorithm 5, some global reductions are needed to be implemented. For example, the computation of the right hand side norm, etc.

3 Numerical Results

In this section, one test model and two suites of real application linear systems are used to test the effectiveness of the local character-based method. The test model is a two-dimensional nonlinear heat conduction equation. The two suites of real application linear systems are respectively the multi-group radiation diffusion equations and the three temperature energy equations.

All experiments are carried out on an E5-2620 CPU, which is clocked at 2.10 GHz and has a 16 GB memory. E5-2620 has a total of 12 cores on two sockets.

The local character-based method is compared with the typical solution method. In the test, the compared method is a BoomerAMG preconditioned GMRES method. For BoomerAMG, Falgout coarsening (a combination of CLJP and the classical Ruge-Stuben coarsening) and “classical” interpolation

[18] are used, the maximal number of levels is set to 8, and one symmetric Gauss-Seidel iteration is used for both pre-relaxation and post-relaxation. The maximal GMRES iteration number is set to 80, and the Krylov dimension is 40.

In the local character-based method, the local domain subsystem is solved by the BoomerAMG preconditioned GMRES method with completely the same parameters for solving the global system. One Gauss-Seidel iteration is applied to smooth $~ x$ . In the local character-based methods, the parameters $α$ in Algorithm 2 and $E_{max}$ in Algorithm 3 are predefined by users, and some numerical analysis is given for different $α$ and $E_{max}$ .

For the local character-based method, after solving the sub-domain system, the assembled solution $~ x$ in (3) satisfies the convergence criteria in most of the tests. Therefore, no further iteration is implemented for the global system in most cases.

For the purpose of presenting the results, the following notations are used to represent the specific methods:

$METHOD 0$ : BoomerAMG preconditioned GMRES method.
$METHOD 1$ : Local character-based method in which the local domain is constructed by Algorithm 2.
$METHOD 2$ : Local character-based method in which the local domain is constructed by Algorithm 3.

Furthermore, the following notations will be used to report the results.

$N$ : the scale of the global linear system.
$K$ : the scale of the sub-domain linear system in the local character-based method.
$η = \frac{K}{N} \times 100 %$ : percentage of the scale of the sub-domain linear system to the scale of the global system.
$CPU\scriptsize glb-meth\footnotesizetotal$ : CPU time for solving the linear system by BoomerAMG preconditioned GMRES method.
$CPU\scriptsize loc-meth\footnotesizetotal$ : CPU time for solving the linear system by the local character-based method.
$CPU\scriptsize loc-constr$ : CPU time for constructing the local domain in the local character-based method.
$CPU\scriptsize loc-solv$ : CPU time for solving the sub-domain system in the local character-based method.
$S_{CPU}$ : the speedup of CPU time for solving the linear system by the local character-based method to BoomerAMG preconditioned GMRES method.

3.1 Test of two-dimensional heat conduction equation

Consider the following two-dimensional (2-D) heat conduction equations:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} \frac{\partial T}{\partial t} = \nabla \cdot (κ (T) \nabla T), & 0 < x < 1, 0 < y < 1, 0 < t < \infty T (0, x, y) = T_{0} (x, y), & 0 < x < 1, 0 < y < 1, t = 0 T (t, 0, y) = T_{l}, & x = 0, 0 < t < \infty T (t, 1, y) = T_{r}, & x = 1, 0 < t < \infty \frac{\partial T}{\partial y} = 0, & y = 0, y = 1, 0 < t < \infty \end{matrix}

In this test, $T_{l}$ , $T_{r}$ , $κ (T)$ and $T_{0} (x, y)$ are set respectively as

T_{l} = 1, T_{r} = 1 \times 10^{- 4}, κ (T) = T^{3.5}, T_{0} (x, y) = e^{- 100 x} \times T_{l} + T_{r} .

By using the backward Euler method to discretize in time, one obtains

\frac{T^{n + 1} - T^{n}}{Δ t} = {(\nabla \cdot (κ (T) \nabla T))}^{n + 1}

(8)

The mesh is $99 \times 99$ , and a five-point finite difference method is used to discretize in space, and a nonlinear system is obtained.

In this test, the time step size is $Δ t = 10^{- 2}$ . The simulation is implemented from physical time 0 to 1, and the total number of time steps is 100. The Picard method is used to solve the system, and the linear system in Picard iteration will be solved by $METHOD 0$ , $METHOD 1$ , or $METHOD 2$ . The parameter $α$ in $METHOD 1$ is set to $10^{- 4}$ , and the $E_{max}$ in $METHOD 2$ is set to 1. For clarity, the time advancing and nonlinear iteration process is given in Algorithm 6.

The $n$ loop in Algorithm 6 is for time advancing, and the $s$ loop is for nonlinear iteration. Eq. (9) is the Picard linearized system. The convergence tolerance for the linear system is $ϵ = 10^{- 10}$ , and the stopping criterion for the nonlinear Picard iteration is ${∥ ∥ T^{n + 1, s + 1} - T^{n + 1, s} ∥ ∥}_{2}^{n + 1, s + 1} < 10^{- 8}$ . In the computation, the previous nonlinear iterate is used as the initial guess for solving the current linear system.

1:Input: The number of mesh for

x

and

y

directions,

N_{x}

N_{y}

, time step size

Δ t

and the initial

T^{0}

2:for

n = 0, 1, \dots

T^{n + 1, 0} = T^{n}

4: for

s = 0, 1, \dots

5: Obtain

T^{n + 1, s + 1}

by solving the following linear equations

\frac{T^{n + 1, s + 1} - T^{n}}{Δ t} = \nabla \cdot (κ (T^{n + 1, s}) \nabla T^{n + 1, s + 1})

(9)

6: if (nonlinear iteration converge) then

7: break;

8: end if

T^{n + 1, s} = T^{n + 1, s + 1}

10: end for

11:end for

Algorithm 6 Solution for 2-D heat conduction equation

First a comparison between the solutions obtained by the three different methods is conducted. Fig. 2 shows the difference in solution obtained by different methods as a function of time, given by

Here, $T_{METHOD 0}$ , $T_{METHOD 1}$ , and $T_{METHOD 2}$ are the solutions obtained by $METHOD 0$ , $METHOD 1$ , and $METHOD 2$ , respectively, at each time step. From this figure, one can see that for all 100 time steps, the maximal value of $e_{1}^{0}$ and $e_{2}^{0}$ are, respectively, $7.6 \times 10^{- 9}$ and $7.0 \times 10^{- 9}$ . This shows that the approximation solution obtained by different methods are very close and they can be considered almost the same.

In the test, if $METHOD 0$ , $METHOD 1$ , or $METHOD 2$ is used to solve all the linear systems for all 100 time steps, the total CPU time for solving the linear systems is respectively $39.19 s$ , $28.09 s$ , and $29.76 s$ . Therefore, the speedups of $METHOD 1$ and $METHOD 2$ to $METHOD 0$ are about 1.40 and 1.32, respectively.

Figure 2: The differences of the solution obtained by $METHOD 0$ , $METHOD 1$ , and $METHOD 2$ for two dimensional heat conduction problem.

Fig. 3 shows the average speedups of the two local character-based methods at each time step. At same time, the average percentage of the scale of the sub-domain linear system to the scale of the global system, $η$ , for two local character-based methods is shown in the figure. From this figure, one can see that both methods have the similar performance in the simulation. With time advancing, the speedup deceases, and $η$ increases. At later simulation time, $METHOD 2$ performs a little better than $METHOD 1$ .

Figure 3: The average speedup and average percentage of the scale of the sub-domain linear system to the scale of the global system for two local character-based methods.

Fig. 4 shows the specific speedup and percentage of the scale of the sub-domain linear system to the scale of the global system at each nonlinear iteration for the two local character-based methods. In the figure, the results at four time steps are given. From this figure, one can see that $METHOD 1$ have a relatively stable performance at each time step. At the same time, the percentage of the sub-domain linear system is stable for all nonlinear iterations at each time. With time advancing, the performance of $METHOD 1$ decreases. For $METHOD 2$ , one can see that it has similar performance for the early 6 nonlinear iterations at all time steps. For the later nonlinear iterations, $METHOD 2$ performs unstable. At about three nonlinear iterations, its performance deteriorates.

Figure 4: The percentage of the scale of the subsystem and the speedup of the local character-based methods for solving linear systems of the 2-D heat conduction equations at time steps 25, 50, 75, and 100. The upper is the result of $METHOD 1$ , the lower is the result of $METHOD 2$ .

For 2-D heat conduction equation, the local domains determined by $METHOD 1$ and $METHOD 2$ are almost the same. As examples, Fig. 5 shows the specific local domain in the first nonlinear iteration at time steps 25, 50, 75, and 100. At the same time, the solutions at these two time steps are shown in the figure. One see can that the local domain $Ω1\scriptsize local$ and $Ω2\scriptsize local$ are almost the same, and they just contain the region where the solution varies greatly. Here $Ω1\scriptsize local$ represents the local domain obtained by Algorithm 2 and $Ω2\scriptsize local$ represents the local domain obtained by Algorithm 3.

Figure 5: The local domain for two dimensional heat conduction problem at time steps 25, 50, 75, and 100 (the first nonlinear iteration).

3.2 Scalability of the local character-based methods for the 2-D heat conduction equation

In this section, the scalability test of the local character-based methods are done for the 2-D heat conduction equation. In the test, the scale on each processor is fixed as $512 \times 512$ , and the number of processors varies from 1 to 256. For all tests, the time step size is fixed as $Δ t = 10^{- 2}$ , and the number of time steps is 100. The test results are listed in Table 2. In the table, the following notations are used:

$CPU\footnotesizetotal$ : CPU time for solving all the linear systems in the simulation.
$N\footnotesize{eq}$ : Total number of linear equations in the simulation.
$¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯CPU\footnotesizetotal$ : The average CPU time for solving one linear system.
$S_{CPU}$ : The speedup of $METHOD 1$ and $METHOD 2$ to $METHOD 0$ (based on the average CPU time for solving one linear system).

#cores	Method	$N$	$CPU\footnotesizetotal$	$N\footnotesize{eq}$	$¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯CPU\footnotesizetotal$	$S_{CPU}$
	$METHOD 0$		1087.45	2757	0.39	-
1	$METHOD 1$	$512^{2}$	693.38	2757	0.25	1.56
	$METHOD 2$		760.71	2758	0.28	1.42
	$METHOD 0$		2364.74	3281	0.72	-
4	$METHOD 1$	$1024^{2}$	927.42	3288	0.28	2.55
	$METHOD 2$		1148.46	3319	0.35	2.06
	$METHOD 0$		3783.12	3655	1.04	-
16	$METHOD 1$	$2048^{2}$	1105.13	3770	0.29	3.42
	$METHOD 2$		1608.81	3844	0.42	2.35
	$METHOD 0$		4172.75	3955	1.06	-
64	$METHOD 1$	$4096^{2}$	1100.34	4253	0.26	3.79
	$METHOD 2$		1789.89	4293	0.42	2.33
	$METHOD 0$		4621.16	4052	1.14	-
256	$METHOD 1$	$8192^{2}$	1283.76	4741	0.27	3.60
	$METHOD 2$		2173.07	4676	0.46	2.13

Table 2: Scalability test of

METHOD 0

METHOD 1

, and

METHOD 2

for 2-D heat conduction equation

From Table 2, one can see that $METHOD 1$ and $METHOD 2$ scale better than $METHOD 0$ by comparing the the total CPU time for solving all the linear equations, $CPU\footnotesizetotal$ . By comparing the numbers of linear equations for different methods, we find that more linear equations are solved for $METHOD 1$ and $METHOD 2$ than $METHOD 0$ . By comparing the average CPU time for solving one linear equation, one can see that $METHOD 1$ and $METHOD 2$ perform much better than $METHOD 0$ . Correspondingly, the speedups of $METHOD 1$ and $METHOD 2$ (for solving one linear equation) are from 1.42 to 3.79. With the increase of the scale, the speedup increases for both $METHOD 1$ and $METHOD 2$ .

3.3 Test of multi-group radiation diffusion equations

In this subsection, the local character-based method will be used to solve the multi-group radiation diffusion equations in real applications. The linear equations for the test are extracted from the radiation hydrodynamics code—Lared-integration, which is used for inertial confinement fusion (ICF) study [21].

3.3.1 Multi-group radiation diffusion equations

The multi-group radiation diffusion equations are a set of $G$ coupled partial differential equations, which are given as follows:

\frac{\partial E_{g}}{\partial t} = \nabla \cdot κ_{g} \nabla E_{g} + 4 π σ_{g} B_{g} - c σ_{g} E_{g}, g = 1, \dots, G

(10)

In the equations,

$G$ is the number of energy groups.
$E_{g}$ is the radiation energy density of the $g$ -th group.
$κ_{g}$ is the radiation diffusion coefficient of the $g$ -th group.
$B_{g}$ is the integration of the Planck function of the $g$ -th group.
$σ_{g}$ is the cross section of the $g$ -th group.

In real applications, the multi-group radiation diffusion equations are coupled to electron-ion equations, and they form a strong nonlinear system. For more details about the radiation diffusion model, see [11, 25]. We will focus on the solution for the multi-group radiation diffusion equations in this paper. Equations (10) are a time-dependent nonlinear system. The backward Euler method is used to discretize in time. For solving the nonlinear system after temporal discretization, a nonlinear iteration (based on Picard-Newton linearization) is used, and in each nonlinear iteration, the source iteration method is used to solve the multi-group radiation diffusion equations [11, 16]. In each source iteration, the $G$ radiation diffusion equations are decoupled, and each group radiation diffusion equation can be solved independently. By considering the time advancing, Algorithm 7 can be used to describe the main solution procedure for multi-group radiation diffusion equations.

1:for (Time stepping loop) do

2: for (Nonlinear iteration loop) do

3: for (Source iteration loop) do

4: Solve

G

linearized multi-group radiation diffusion equations.

5: end for

6: end for

7:end for

Algorithm 7 Solution procedure for multi-group radiation diffusion equations

A Lagrangian method is used in the Lared-integration code and the mesh moves with the motion of the fluid. Consequently, the multi-group radiation diffusion equations are discretized on two-dimensional deforming meshes. A nine-point scheme is used to discretize the diffusion equations.

3.3.2 Linear systems for test

In the test, the number of cells is 7040. The number of groups is $G = 64$ , and therefore, 64 equations must be solved in each source iteration.

Linear systems are tested for all groups, in a given source iteration of a given nonlinear iteration, at time steps 1,000 and 5,000. For notation, each linear system is denoted as T $n_{1}$ N $n_{2}$ S $n_{3}$ G $n_{4}$ , which means the linear system of group $n_{4}$ in the $n_{3}$ th source iteration of $n_{2}$ th nonlinear iteration at $n_{1}$ th time step. For example, T1000N0S0G0 represents the linear system of $0$ -th group in $0$ -th source iteration of $0$ -th nonlinear iteration at time step 1,000. All the coefficient matrices of the linear system have the same sparsity pattern, shown in Fig. 6.

Figure 6: The sparsity pattern of the matrix of the multi-group radiation diffusion linear system.

3.3.3 Test results

We test a total of 128 linear systems on two time steps, i.e., T1000N0S0G0-G63 and T5000N0S0G0-G63. In this test, the parameter $α$ in $METHOD 1$ is set to $10^{- 10}$ , and $E_{max}$ in $METHOD 2$ is set to 1. The convergence tolerance $ϵ =$

A local character based method for solving linear systems of radiation diffusion problems

Authors

Virtual element method for the system of time dependent nonlinear convection-diffusion-reaction equation

Convergence Analysis of a Cooperative Diffusion Gauss-Newton Strategy

A preconditioner based on sine transform for two-dimensional Riesz space factional diffusion equations in convex domains

Time-dependent stochastic basis adaptation for uncertainty quantification

A modified Kačanov iteration scheme with application to quasilinear diffusion models

A note on parallel preconditioning for the all-at-once solution of Riesz fractional diffusion equations

Diffusion synthetic acceleration for heterogeneous domains, compatible with voids

1 Introduction

2 The local character-based method

2.1 The local character-based algorithm

2.2 Gradient-based local domain construction

2.3 Residual-based local domain construction

2.3.1 Finding bad points

2.3.2 Expanding the local domain

Definition 2.1

Example 1

2.4 Parallelization of the local character-based methods

3 Numerical Results

3.1 Test of two-dimensional heat conduction equation

3.2 Scalability of the local character-based methods for the 2-D heat conduction equation

3.3 Test of multi-group radiation diffusion equations

3.3.1 Multi-group radiation diffusion equations

3.3.2 Linear systems for test

3.3.3 Test results

A local character based method for solving linear systems of radiation diffusion problems

Authors

Related Research

Virtual element method for the system of time dependent nonlinear convection-diffusion-reaction equation

Convergence Analysis of a Cooperative Diffusion Gauss-Newton Strategy

A preconditioner based on sine transform for two-dimensional Riesz space factional diffusion equations in convex domains

Time-dependent stochastic basis adaptation for uncertainty quantification

A modified Kačanov iteration scheme with application to quasilinear diffusion models

A note on parallel preconditioning for the all-at-once solution of Riesz fractional diffusion equations

Diffusion synthetic acceleration for heterogeneous domains, compatible with voids

This week in AI

1 Introduction

2 The local character-based method

2.1 The local character-based algorithm

2.2 Gradient-based local domain construction

2.3 Residual-based local domain construction

2.3.1 Finding bad points

2.3.2 Expanding the local domain

Definition 2.1

Example 1

2.4 Parallelization of the local character-based methods

3 Numerical Results

3.1 Test of two-dimensional heat conduction equation

3.2 Scalability of the local character-based methods for the 2-D heat conduction equation

3.3 Test of multi-group radiation diffusion equations

3.3.1 Multi-group radiation diffusion equations

3.3.2 Linear systems for test

3.3.3 Test results