**3.1 Explicit block**

6 Will-be-set-by-IN-TECH

Explicit Hydrodynamics Block Based on Second order R−K method

> n+1 , (ρu ) n+1 , E\* )

*<sup>∂</sup><sup>r</sup>* <sup>=</sup> 0, (20)

*<sup>∂</sup><sup>r</sup>* ), (21)

k

Input T Return ( ρ

ρ ρu E

For k =1,...,kmax

T is available k

n+1= cv ρn+1T

Calculate T

Tk+1= Tk+ <sup>δ</sup> <sup>T</sup>

δ k

k

E ρ u

hydrodynamics equations can be written as

ρ (n+1)

Newton Iteration

Res =

E

diffusion equation becomes

augmented by Equation (21).

where *V* = <sup>4</sup>

n+1 t

End

t n (n)

Call Hydrodynamics Block with T Form Non−Linear Residual

> k T n+1

k ρ T

i i = 1,2,...,N

n+1 <sup>2</sup> (u n+1 <sup>+</sup> <sup>ρ</sup>

n+1 cv − E\* Δ t

, ,

<sup>ρ</sup> n T

−( RHS( ) + RHS( ))/2

k+1+ <sup>ρ</sup>n+1(u n+1 ) /2

i i = 1,2,...,N

Fig. 1. Flowchart of the second order self-consistent IMEX algorithm

*∂***U** *<sup>∂</sup><sup>t</sup>* <sup>+</sup>

> *∂E <sup>∂</sup><sup>t</sup>* <sup>=</sup> *<sup>∂</sup>*

)/2

n ρ

k

2

accounting for the effects of radiation transport (diffusion equation). For instance, the pure

where **U** = (*ρ*, *ρu*, *E*)*T*, **F**(**U**)=(*ρu*, *ρu*2, *u*(*E* + *p*))*T*, and **G**(**U**)=(0, *p*, 0)*T*. Then the

*<sup>∂</sup><sup>V</sup>* (*A<sup>κ</sup>*

and *A* = 4*πr*<sup>2</sup> is the associated cross-sectional area. Notice that the total energy density, *E*, obtained by Equation (20) just represents the hydrodynamics component and it must be

Our algorithm consists of an explicit and an implicit block. The explicit block solves Equation (20) and the implicit block solves Equation (21). We will briefly describe these algorithm blocks in the following subsections. However, we note again that the explicit block is embedded within the implicit block as part of a nonlinear function evaluation as it is depicted in Fig. 1. This is done to obtain a nonlinearly converged algorithm that leads to second order calculations. We also note that similar discretizations, but without converging nonlinearities, can lead to order reduction in time convergence (Bates et al., 2001). Before we go into details of the individual algorithm blocks, we would like to present a flow diagram that illustrates the

*∂***G**

*∂T*

<sup>3</sup>*πr*<sup>3</sup> is the generalized volume coordinate in one-dimensional spherical geometry,

*∂*(*A***F**) *<sup>∂</sup><sup>V</sup>* <sup>+</sup> Our explicit time discretization is based on a second order TVD Runge-Kutta method (Gottlieb & Shu, 1998; Gottlieb et al., 2001; Shu & Osher, 1988; 1989). The main reason why we choose this methodology is that it preserves the strong stability properties of the explicit Euler method. This is important because it is well known that solutions to the conservation laws usually involve discontinuities (e.g, shock or contact discontinuities) and (Gottlieb & Shu, 1998; Gottlieb et al., 2001) suggest that a time integration method which has the strong stability preserving property leads to non-oscillatory calculations (especially at shock or contact discontinuities).

A second order two-step TVD Runge-Kutta method for (20) can be cast as Step-1 :

$$\begin{aligned} \rho^1 &= \rho^n - \Delta t \frac{1}{r^2} \frac{\partial}{\partial r} (r^2 \rho u)^n, \\ (\rho u)^1 &= (\rho u)^n - \Delta t [\frac{1}{r^2} \frac{\partial}{\partial r} (r^2 \rho u^2) + \frac{\partial p}{\partial r}]^n, \\ E^1 &= E^n - \Delta t \{\frac{1}{r^2} \frac{\partial}{\partial r} [r^2 u (E+p)]\}^n. \end{aligned} \tag{22}$$

Step-2 :

$$\begin{split} \rho^{n+1} &= \frac{\rho^{n} + \rho^{1}}{2} - \frac{\Delta t}{2} \frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} \rho u)^{1}, \\ (\rho u)^{n+1} &= \frac{(\rho u)^{n} + (\rho u)^{1}}{2} - \frac{\Delta t}{2} \{\frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} \rho u^{2})^{1} + \frac{\partial}{\partial r} (\rho^{1} R \mathbf{T}^{\mathbf{n}+1})\}, \\ E^{\*} &= \frac{E^{n} + E^{1}}{2} - \frac{\Delta t}{2} \{\frac{1}{r^{2}} \frac{\partial}{\partial r} [r^{2} u^{1} (c\_{v} \rho^{1} \mathbf{T}^{\mathbf{n}+1} + \frac{1}{2} \rho^{1} (u^{1})^{2} + \rho^{1} R \mathbf{T}^{\mathbf{n}+1})\}). \end{split} \tag{23}$$

We used the following equation of state relations in (22)- (23);

$$p = \rho \text{RTE} = c\_v \rho T + \frac{1}{2} \rho u^2 \,, \tag{24}$$

where *cv* = *<sup>R</sup> <sup>γ</sup>*−<sup>1</sup> is the fluid specific heat with *<sup>R</sup>* being the universal gas constant. This explicit algorithm block interacts with the implicit block through the highlighted *Tn*+<sup>1</sup> terms in Equation (23). We can observe that the implicit equation (21) is practically solved for *T* by using the energy relation. Therefore, the explicit block is continuously impacted by the

where *<sup>∂</sup><sup>p</sup>*

where

where

where

**3.2 Implicit block**

(*cvρn*+1*Tn*+<sup>1</sup> + <sup>1</sup>

*∂T*

*∂ <sup>∂</sup><sup>V</sup>* (*A<sup>κ</sup>*

*∂ρ* <sup>=</sup> *RT* in this study. **<sup>U</sup>***<sup>R</sup>*

edge from the right and left side, i.e,

*<sup>i</sup>*+1/2 and **<sup>U</sup>***<sup>L</sup>*

Non-Linear Heat Conduction and Stiff Source Terms (Radiation Hydrodynamics)

*<sup>i</sup>*+1/2 <sup>=</sup> **<sup>U</sup>***i*+<sup>1</sup> <sup>−</sup> <sup>Δ</sup>*<sup>r</sup>*

<sup>301</sup> An IMEX Method for the Euler Equations That Posses Strong

⎧ ⎨ ⎩

*<sup>a</sup>* <sup>=</sup> **<sup>U</sup>***i*+<sup>1</sup> <sup>−</sup> **<sup>U</sup>***<sup>i</sup>*

*<sup>b</sup>* <sup>=</sup> **<sup>U</sup>***<sup>i</sup>* <sup>−</sup> **<sup>U</sup>***i*−<sup>1</sup>

⎛ ⎝

2 *∂*

Notice that this implicit discretization resembles to the *Crank-Nicolson method* (Strikwerda, 1989; Thomas, 1998). We solve Equation (33) iteratively for *Tn*<sup>+</sup>1. The nonlinear solver needed by Equation (33) is based on the Jacobian-Free Newton Krylov method which is described in the next subsection. When the Newton method converges all the nonlinearities in this

> ⎛ ⎝

*ρn*+<sup>1</sup> (*ρu*)*n*+<sup>1</sup> *En*+<sup>1</sup>

⎞ ⎠ .

*ρn*+<sup>1</sup> (*ρu*)*n*+<sup>1</sup> *E*∗

*<sup>∂</sup><sup>V</sup>* (*Aκn*+<sup>1</sup> *<sup>∂</sup>Tn*+<sup>1</sup>

*<sup>∂</sup><sup>r</sup>* )*<sup>i</sup>* <sup>+</sup>

1 2 *∂*

<sup>−</sup> *Ai*−1/2*κi*−1/2(*Ti* <sup>−</sup> *Ti*−1)/Δ*<sup>r</sup>* Δ*Vi*

*<sup>∂</sup><sup>V</sup>* (*Aκ<sup>n</sup> <sup>∂</sup>T<sup>n</sup>*

*<sup>∂</sup><sup>r</sup>* )*i*, (33)

. (34)

⎞ ⎠ .

**<sup>U</sup>***<sup>n</sup>* <sup>→</sup> **<sup>U</sup>**<sup>∗</sup> <sup>=</sup>

Δ*r*

0 if *ab* ≤ 0,

*<sup>i</sup>*+1/2 = **U***<sup>i</sup>* +

<sup>2</sup> **<sup>U</sup>***r*,*i*<sup>+</sup>1,

*a* if |*a*| < |*b*| and *ab* > 0, *b* if |*b*| < |*a*| and *ab* > 0,

**U***<sup>R</sup>*

**U***<sup>L</sup>*

**U***r*,*<sup>i</sup>* = *minmod*(*a*, *b*) =

The explicit block produces the following solution vector

This information is used to discretize Equation (21) as follows,

<sup>2</sup> *<sup>ρ</sup>n*+1(*un*+1)<sup>2</sup> <sup>−</sup> *<sup>E</sup>*∗)*<sup>i</sup>* <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> <sup>1</sup>

*<sup>∂</sup><sup>r</sup>* )*<sup>i</sup>* <sup>=</sup> *Ai*<sup>+</sup>1/2*κi*+1/2(*Ti*+<sup>1</sup> <sup>−</sup> *Ti*)/Δ*<sup>r</sup>* Δ*Vi*

discretization, we obtain the following fully updated solution vector,

**<sup>U</sup>**<sup>∗</sup> <sup>→</sup> **<sup>U</sup>***n*+<sup>1</sup> <sup>=</sup>

*<sup>i</sup>*+1/2 are the interpolated values at (*<sup>i</sup>* <sup>+</sup> 1/2)*th* cell

<sup>2</sup> **<sup>U</sup>***r*,*i*, (29)

<sup>Δ</sup>*<sup>r</sup>* , (31)

<sup>Δ</sup>*<sup>r</sup>* . (32)

(30)

n

Fig. 2. Computational Conventions.

implicit *Tn*+<sup>1</sup> solutions at each non-linear Newton iteration. This provides the tight nonlinear coupling between the two algorithm blocks. Notice that the *kth* nonlinear Newton iteration of the implicit block corresponds to *<sup>T</sup>n*+<sup>1</sup> <sup>←</sup> *<sup>T</sup><sup>k</sup>* and *<sup>k</sup>* <sup>→</sup> (*<sup>n</sup>* <sup>+</sup> <sup>1</sup>) upon the convergence of the Newton method (refer to Fig. 1). Also, the ∗ values in Equation (23) are predicted intermediate values and later they are corrected by the implicit block which is given in the next subsection. One observation about this algorithm block is that some calculations are redundant related to Equation (22). In other words, Equation (22) can be computed only once at the beginning of each Newton iteration, because the non-linear iterations do not impact (22). This can lead overall less number of function evaluations.

Now we shall describe how we evaluate the numerical fluxes needed by Equations (22) and (23). For simplicity, we consider (20) to describe our fluxing procedure. Basically, it is based on the Local Lax Friedrichs (LLF) method (we refer to (LeVeque, 1998; Thomas, 1999) for the details of the LLF method and for more information in regards to the explicit discretizations of conservation laws). For instance, if we consider the following simple discretization for Equation (20),

$$\mathbf{U}\_{i}^{1} = \mathbf{U}\_{i}^{n} - \frac{\Delta t}{\Delta V\_{i}} (A\_{i+1/2} F\_{i+1/2}^{\rm n} - A\_{i-1/2} F\_{i-1/2}^{\rm n}) - \frac{\Delta t}{\Delta r} (\mathbf{G}\_{i+1/2}^{\rm n} - \mathbf{G}\_{i+1/2}^{\rm n}),\tag{25}$$

where <sup>Δ</sup>*Vi* = *<sup>V</sup>*(*ri*<sup>+</sup>1/2) − *<sup>V</sup>*(*ri*−1/2), *Ai*±1/2 = *<sup>A</sup>*(*ri*±1/2), and indices *<sup>i</sup>* and *<sup>i</sup>* + 1/2 represent cell center and cell edge values respectively (refer to Fig. 2), then the *Local Lax Friedrichs method* defines *Fi*<sup>+</sup>1/2 and *Gi*<sup>+</sup>1/2 as

$$F\_{i+1/2} = \frac{F(\mathbf{U}\_{i+1/2}^R) + F(\mathbf{U}\_{i+1/2}^L)}{2} - a\_{i+1/2} \frac{\mathbf{U}\_{i+1/2}^R - \mathbf{U}\_{i+1/2}^L}{2},\tag{26}$$

$$G\_{i+1/2} = \frac{G(\mathbf{U}\_{i+1/2}^R) + G(\mathbf{U}\_{i+1/2}^L)}{2},\tag{27}$$

where *<sup>α</sup>* <sup>=</sup> *max*{|*λ<sup>L</sup>* <sup>1</sup> <sup>|</sup>, <sup>|</sup>*λ<sup>R</sup>* <sup>1</sup> <sup>|</sup>, <sup>|</sup>*λ<sup>L</sup>* <sup>2</sup> <sup>|</sup>, <sup>|</sup>*λ<sup>R</sup>* <sup>2</sup> <sup>|</sup>, <sup>|</sup>*λ<sup>L</sup>* <sup>3</sup> <sup>|</sup>, <sup>|</sup>*λ<sup>R</sup>* <sup>3</sup> |} in which *λ*<sup>1</sup> = *u* − *c*, *λ*<sup>2</sup> = *u*, *λ*<sup>3</sup> = *u* + *c*, and *c* is the sound speed. The sound speed is defined by

$$c = \sqrt{\frac{\partial p}{\partial \rho}}\,\,\,\,\,\,\tag{28}$$

where *<sup>∂</sup><sup>p</sup> ∂ρ* <sup>=</sup> *RT* in this study. **<sup>U</sup>***<sup>R</sup> <sup>i</sup>*+1/2 and **<sup>U</sup>***<sup>L</sup> <sup>i</sup>*+1/2 are the interpolated values at (*<sup>i</sup>* <sup>+</sup> 1/2)*th* cell edge from the right and left side, i.e,

$$\mathbf{U}\_{i+1/2}^{R} = \mathbf{U}\_{i+1} - \frac{\Delta r}{2} \mathbf{U}\_{r, i+1 \prime}$$

$$\mathbf{U}\_{i+1/2}^{L} = \mathbf{U}\_{i} + \frac{\Delta r}{2} \mathbf{U}\_{r, i\prime} \tag{29}$$

where

8 Will-be-set-by-IN-TECH

Cell Center Cell Edge

implicit *Tn*+<sup>1</sup> solutions at each non-linear Newton iteration. This provides the tight nonlinear coupling between the two algorithm blocks. Notice that the *kth* nonlinear Newton iteration of the implicit block corresponds to *<sup>T</sup>n*+<sup>1</sup> <sup>←</sup> *<sup>T</sup><sup>k</sup>* and *<sup>k</sup>* <sup>→</sup> (*<sup>n</sup>* <sup>+</sup> <sup>1</sup>) upon the convergence of the Newton method (refer to Fig. 1). Also, the ∗ values in Equation (23) are predicted intermediate values and later they are corrected by the implicit block which is given in the next subsection. One observation about this algorithm block is that some calculations are redundant related to Equation (22). In other words, Equation (22) can be computed only once at the beginning of each Newton iteration, because the non-linear iterations do not impact (22). This can lead

Now we shall describe how we evaluate the numerical fluxes needed by Equations (22) and (23). For simplicity, we consider (20) to describe our fluxing procedure. Basically, it is based on the Local Lax Friedrichs (LLF) method (we refer to (LeVeque, 1998; Thomas, 1999) for the details of the LLF method and for more information in regards to the explicit discretizations of conservation laws). For instance, if we consider the following simple discretization for

*<sup>i</sup>*+1/2 <sup>−</sup> *Ai*−1/2*F<sup>n</sup>*

where <sup>Δ</sup>*Vi* = *<sup>V</sup>*(*ri*<sup>+</sup>1/2) − *<sup>V</sup>*(*ri*−1/2), *Ai*±1/2 = *<sup>A</sup>*(*ri*±1/2), and indices *<sup>i</sup>* and *<sup>i</sup>* + 1/2 represent cell center and cell edge values respectively (refer to Fig. 2), then the *Local Lax Friedrichs method*

> *<sup>i</sup>*+1/2) <sup>2</sup> <sup>−</sup> *<sup>α</sup>i*+1/2

> > *<sup>i</sup>*+1/2) + *<sup>G</sup>*(**U***<sup>L</sup>*

*<sup>i</sup>*−1/2) <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>*

Δ*r* (*G<sup>n</sup>*

**U***<sup>R</sup>*

*<sup>i</sup>*+1/2)

*<sup>i</sup>*+1/2 <sup>−</sup> **<sup>U</sup>***<sup>L</sup>*

*<sup>i</sup>*+1/2 <sup>−</sup> *<sup>G</sup><sup>n</sup>*

*i*+1/2

<sup>2</sup> , (27)

*∂ρ* , (28)

<sup>3</sup> |} in which *λ*<sup>1</sup> = *u* − *c*, *λ*<sup>2</sup> = *u*, *λ*<sup>3</sup> = *u* + *c*, and

<sup>2</sup> , (26)

*<sup>i</sup>*+1/2), (25)

i+1/2

r = R **<sup>0</sup>**

t

n

t

n+1

A Computational Cell

: Represents a Cell Centered Quantity at time level n

: Represents a Cell Edge Quantity at time level n

r <sup>i</sup> <sup>r</sup>

r = 0

i

u u

Equation (20),

**U**1 *<sup>i</sup>* <sup>=</sup> **<sup>U</sup>***<sup>n</sup>*

defines *Fi*<sup>+</sup>1/2 and *Gi*<sup>+</sup>1/2 as

where *<sup>α</sup>* <sup>=</sup> *max*{|*λ<sup>L</sup>*

n n

i+1/2

Fig. 2. Computational Conventions.

overall less number of function evaluations.

*<sup>i</sup>* <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>* Δ*Vi*

*Fi*<sup>+</sup>1/2 <sup>=</sup> *<sup>F</sup>*(**U***<sup>R</sup>*

<sup>1</sup> <sup>|</sup>, <sup>|</sup>*λ<sup>R</sup>*

<sup>1</sup> <sup>|</sup>, <sup>|</sup>*λ<sup>L</sup>*

*c* is the sound speed. The sound speed is defined by

<sup>2</sup> <sup>|</sup>, <sup>|</sup>*λ<sup>R</sup>*

(*Ai*+1/2*F<sup>n</sup>*

*<sup>i</sup>*+1/2) + *<sup>F</sup>*(**U***<sup>L</sup>*

*Gi*<sup>+</sup>1/2 <sup>=</sup> *<sup>G</sup>*(**U***<sup>R</sup>*

<sup>2</sup> <sup>|</sup>, <sup>|</sup>*λ<sup>L</sup>*

<sup>3</sup> <sup>|</sup>, <sup>|</sup>*λ<sup>R</sup>*

*c* = *∂p*

$$\mathbf{U}\_{r,i} = \min \text{mod}(a,b) = \begin{cases} a & \text{if } |a| < |b| \text{ and } ab > 0, \\ b & \text{if } |b| < |a| \text{ and } ab > 0, \\ 0 & \text{if } ab \le 0, \end{cases} \tag{30}$$

where

$$a = \frac{\mathbf{U}\_{i+1} - \mathbf{U}\_i}{\Delta r},\tag{31}$$

$$b = \frac{\mathbf{U}\_i - \mathbf{U}\_{i-1}}{\Delta r}.\tag{32}$$

#### **3.2 Implicit block**

The explicit block produces the following solution vector

$$\mathbf{U}^n \to \mathbf{U}^\* = \begin{pmatrix} \rho^{n+1} \\ (\rho u)^{n+1} \\ E^\* \end{pmatrix}.$$

This information is used to discretize Equation (21) as follows,

$$\frac{(c\_v \rho^{n+1} T^{n+1} + \frac{1}{2} \rho^{n+1} (u^{n+1})^2 - E^\*)\_{\mathrm{i}}}{\Delta t} = \frac{1}{2} \frac{\partial}{\partial V} (A \boldsymbol{\kappa}^{n+1} \frac{\partial T^{n+1}}{\partial r})\_{\mathrm{i}} + \frac{1}{2} \frac{\partial}{\partial V} (A \boldsymbol{\kappa}^{\mathrm{n}} \frac{\partial T^{\mathrm{n}}}{\partial r})\_{\mathrm{i}} \tag{33}$$

where

$$\frac{\partial}{\partial V}(A\boldsymbol{\kappa}\frac{\partial T}{\partial r})\_{\boldsymbol{i}} = \frac{A\_{\boldsymbol{i}+1/2}\boldsymbol{\kappa}\_{\boldsymbol{i}+1/2}(T\_{\boldsymbol{i}+1}-T\_{\boldsymbol{i}})/\Delta r}{\Delta V\_{\boldsymbol{i}}} - \frac{A\_{\boldsymbol{i}-1/2}\boldsymbol{\kappa}\_{\boldsymbol{i}-1/2}(T\_{\boldsymbol{i}}-T\_{\boldsymbol{i}-1})/\Delta r}{\Delta V\_{\boldsymbol{i}}}.\tag{34}$$

Notice that this implicit discretization resembles to the *Crank-Nicolson method* (Strikwerda, 1989; Thomas, 1998). We solve Equation (33) iteratively for *Tn*<sup>+</sup>1. The nonlinear solver needed by Equation (33) is based on the Jacobian-Free Newton Krylov method which is described in the next subsection. When the Newton method converges all the nonlinearities in this discretization, we obtain the following fully updated solution vector,

$$\mathbf{U}^\* \to \mathbf{U}^{n+1} = \begin{pmatrix} \rho^{n+1} \\ (\rho u)^{n+1} \\ E^{n+1} \end{pmatrix}.$$

where *�* = <sup>1</sup>

refer to Fig. 1). **Evaluating** *F*(*Tk*) :

*<sup>F</sup>*(*Tk*) = [*cvρn*+1*T<sup>k</sup>*<sup>+</sup> <sup>1</sup>

**3.4 Time step control**

approach.

*<sup>n</sup>*�*v*�<sup>2</sup> <sup>∑</sup>*<sup>n</sup>*

(typically 10−<sup>6</sup> for 64-bit double precision).

**Given** *T<sup>k</sup>* where *k* represents the current Newton iteration.

rather than using the entire system of the governing equations

where the unknown *υ<sup>f</sup>* represents the front velocity. This gives

*υn*

numerical approximation is used to calculate *υ<sup>f</sup>*

*∂E <sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>υ</sup><sup>f</sup>*

*<sup>f</sup>* <sup>=</sup> <sup>∑</sup>(|*E<sup>n</sup>*

Δ*t*

<sup>∑</sup>(|*E<sup>n</sup>*

Then the new time step is determined by the Courant-Friedrichs-Lewy (CFL) condition

**Form** *F*(*Tk*) based on the Crank-Nicolson method,

<sup>2</sup> *<sup>ρ</sup>n*+<sup>1</sup>(*un*+<sup>1</sup>)<sup>2</sup>−*E*<sup>∗</sup> ] <sup>Δ</sup>*<sup>t</sup>* <sup>−</sup> <sup>1</sup>

**Call** Hydrodynamics block with (*ρn*, *un*, *En*, *Tk*) to compute *ρn*<sup>+</sup>1, *un*<sup>+</sup>1, *E*∗.

Non-Linear Heat Conduction and Stiff Source Terms (Radiation Hydrodynamics)

2 *∂ <sup>∂</sup><sup>V</sup>* (*Aκ<sup>k</sup> <sup>∂</sup>T<sup>k</sup>*

*<sup>i</sup>*=<sup>1</sup> *b*|*ui*| + *b*, *n* is the dimension of the linear system and *b* is a constant

whose magnitude is within a few orders of magnitude of the square root of machine roundoff

<sup>303</sup> An IMEX Method for the Euler Equations That Posses Strong

Here, we briefly describe how to form the IMEX function *F*(*T*). We refer *F*(*T*) as the IMEX function, since it uses both explicit (hydrodynamics) and implicit (diffusion) information. Notice that for a method that uses all implicit information, *F*(*T*) would correspond to a regular nonlinear residual function. The following pseudo code describes how to form *F*(*T*) (we also

> *<sup>∂</sup><sup>r</sup>* ) <sup>−</sup> <sup>1</sup> 2 *∂ <sup>∂</sup><sup>V</sup>* (*Aκ<sup>n</sup> <sup>∂</sup>T<sup>n</sup>*

It is important to note that we are not iterating between the implicit and explicit blocks. Instead we are executing the explicit block inside of a nonlinear function evaluation defined by *F*(*Tk*). The unique properties of JFNK allow us to perform a Newton iteration on this IMEX function, and thus JFNK is a required component of this nonlinearly converged IMEX

In this section, we describe two procedures to determine the computational time steps that are used in our test calculations. The first one was originally proposed by (Rider & Knoll, 1999). The idea is to estimate the dominant wave propagation speed in the problem. In one dimension this involves calculating the ratio of temporal to spatial derivatives of the dependent variables. In principle, it is sufficient to consider the following hyperbolic equation

*∂E*

*<sup>υ</sup><sup>f</sup>* <sup>=</sup> <sup>−</sup> *<sup>∂</sup>E*/*∂<sup>t</sup> ∂E*/*∂r*

As noted in Rider & Knoll (1999), to avoid problems from lack of smoothness the following

*<sup>i</sup>* <sup>−</sup> *<sup>E</sup>n*−<sup>1</sup>

*<sup>i</sup>*+<sup>1</sup> <sup>−</sup> *<sup>E</sup><sup>n</sup>*

*<sup>n</sup>*+<sup>1</sup> <sup>=</sup> *<sup>C</sup>*� <sup>Δ</sup>*<sup>r</sup>* � *υn f*

*<sup>i</sup>* |/Δ*t*)

*<sup>i</sup>*−1|/2Δ*r*)

*<sup>∂</sup><sup>r</sup>* ).

*<sup>∂</sup><sup>r</sup>* <sup>=</sup> 0, (41)

. (42)

. (43)

, (44)

#### **3.3 The Jacobian-Free Newton Krylov method and forming the IMEX function**

The Jacobian-Free Newton Krylov method (e.g, refer to (Brown & Saad, 1990; Kelley, 2003; Knoll & Keyes, 2004)) is a combination of the Newton method that solves a system of nonlinear equations and a Krylov subspace method that solves the Newton correction equations. With this method, Newton-like super-linear convergence is achieved in the nonlinear iterations, without the complexity of forming or storing the Jacobian matrix. The effects of the Jacobian matrix are probed only through approximate matrix-vector products required in the Krylov iterations. Below, we provide more details about this technique.

The Newton method solves **F**(*T*) = 0 (e.g, assume Equation (33) is written in this form) iteratively over a sequence of linear system defined by

$$\begin{aligned} \mathbf{J}(T^k)\delta T^k &= -\mathbf{F}(T^k),\\ T^{k+1} &= T^k + \delta T^k, \qquad k = 0, 1, \cdots \end{aligned} \tag{35}$$

where **J**(*Tk*) = *<sup>∂</sup>***<sup>F</sup>** *<sup>∂</sup><sup>T</sup>* is the Jacobian matrix and *<sup>δ</sup>T<sup>k</sup>* is the update vector. The Newton iteration is terminated based on a required drop in the norm of the nonlinear residual, i.e,

$$\|\|\mathbf{F}(T^k)\|\|\_2 < \operatorname{tol}\_{\operatorname{res}} \|\|\mathbf{F}(T^0)\|\|\_2 \tag{36}$$

where *tolres* is a given tolerance. The linear system, Newton correction equation (35), is solved by using the Arnoldi based Generalized Minimal RESidual method (GMRES)(Saad, 2003) which belongs to the general class of the Krylov subspace methods(Reid, 1971). We note that these subspace methods are particularly suitable choice when dealing with non-symmetric linear systems. In GMRES, an initial linear residual, **r**0, is defined for a given initial guess *δT*0,

$$\mathbf{r}\_0 = -\mathbf{F}(T) - \mathbf{J}\delta T\_0. \tag{37}$$

Here we dropped the index *k* convention since the Krylov (GMRES) iteration is performed at a fixed *k*. Let *j* be the Krylov iteration index. The *j th* Krylov iteration minimizes �**J***δTj* + **F**(*T*)�<sup>2</sup> within a subspace of small dimension, relative to *n* (the number of unknowns), in a least-squares sense. *δTj* is drawn from the subspace spanned by the Krylov vectors, {**r**0,**Jr**0,**J**2**r**0, ··· ,**J***j*−1**r**0} , and can be written as

$$
\delta T\_j = \delta T\_0 + \sum\_{i=0}^{j-1} \beta\_i(\mathbf{J})^i \mathbf{r}\_{0\prime} \tag{38}
$$

where the scalar *β<sup>i</sup>* minimizes the residual. The Krylov iteration is terminated based on the following inexact Newton criteria (Dembo, 1982)

$$\|\mathbf{J}\delta T\_{\rangle} + \mathbf{F}(T)\|\_{2} < \gamma \|\mathbf{F}(T)\|\_{2} \tag{39}$$

where the parameter *γ* is set in terms of how tight the linear solver should converge at each Newton iteration (we typically use *γ* = 10−3). One particularly attractive feature of this methodology is that it does not require forming the Jacobian matrix. Instead, only matrix-vector multiplications, **<sup>J</sup>***v*, are needed, where *<sup>v</sup>* ∈ {**r**0,**Jr**0,**J**2**r**0, ···} . This leads to the so-called *Jacobian-Free* implementations in which the action of the Jacobian matrix can be approximated by

$$\mathbf{J}v = \frac{\mathbf{F}(T + \epsilon v) - \mathbf{F}(T)}{\varepsilon},\tag{40}$$

where *�* = <sup>1</sup> *<sup>n</sup>*�*v*�<sup>2</sup> <sup>∑</sup>*<sup>n</sup> <sup>i</sup>*=<sup>1</sup> *b*|*ui*| + *b*, *n* is the dimension of the linear system and *b* is a constant whose magnitude is within a few orders of magnitude of the square root of machine roundoff (typically 10−<sup>6</sup> for 64-bit double precision).

Here, we briefly describe how to form the IMEX function *F*(*T*). We refer *F*(*T*) as the IMEX function, since it uses both explicit (hydrodynamics) and implicit (diffusion) information. Notice that for a method that uses all implicit information, *F*(*T*) would correspond to a regular nonlinear residual function. The following pseudo code describes how to form *F*(*T*) (we also refer to Fig. 1).

**Evaluating** *F*(*Tk*) :

10 Will-be-set-by-IN-TECH

The Jacobian-Free Newton Krylov method (e.g, refer to (Brown & Saad, 1990; Kelley, 2003; Knoll & Keyes, 2004)) is a combination of the Newton method that solves a system of nonlinear equations and a Krylov subspace method that solves the Newton correction equations. With this method, Newton-like super-linear convergence is achieved in the nonlinear iterations, without the complexity of forming or storing the Jacobian matrix. The effects of the Jacobian matrix are probed only through approximate matrix-vector products required in the Krylov iterations. Below, we provide more details about this technique. The Newton method solves **F**(*T*) = 0 (e.g, assume Equation (33) is written in this form)

),

where *tolres* is a given tolerance. The linear system, Newton correction equation (35), is solved by using the Arnoldi based Generalized Minimal RESidual method (GMRES)(Saad, 2003) which belongs to the general class of the Krylov subspace methods(Reid, 1971). We note that these subspace methods are particularly suitable choice when dealing with non-symmetric linear systems. In GMRES, an initial linear residual, **r**0, is defined for a given initial guess *δT*0,

Here we dropped the index *k* convention since the Krylov (GMRES) iteration is performed

�**J***δTj* + **F**(*T*)�<sup>2</sup> within a subspace of small dimension, relative to *n* (the number of unknowns), in a least-squares sense. *δTj* is drawn from the subspace spanned by the Krylov vectors,

> *j*−1 ∑ *i*=0

where the scalar *β<sup>i</sup>* minimizes the residual. The Krylov iteration is terminated based on the

where the parameter *γ* is set in terms of how tight the linear solver should converge at each Newton iteration (we typically use *γ* = 10−3). One particularly attractive feature of this methodology is that it does not require forming the Jacobian matrix. Instead, only matrix-vector multiplications, **<sup>J</sup>***v*, are needed, where *<sup>v</sup>* ∈ {**r**0,**Jr**0,**J**2**r**0, ···} . This leads to the so-called *Jacobian-Free* implementations in which the action of the Jacobian matrix can be

**<sup>J</sup>***<sup>v</sup>* <sup>=</sup> **<sup>F</sup>**(*<sup>T</sup>* <sup>+</sup> *�v*) <sup>−</sup> **<sup>F</sup>**(*T*)

*βi*(**J**)*<sup>i</sup>*

*δTj* = *δT*<sup>0</sup> +

*<sup>∂</sup><sup>T</sup>* is the Jacobian matrix and *<sup>δ</sup>T<sup>k</sup>* is the update vector. The Newton iteration

*<sup>T</sup>k*+<sup>1</sup> <sup>=</sup> *<sup>T</sup><sup>k</sup>* <sup>+</sup> *<sup>δ</sup>Tk*, *<sup>k</sup>* <sup>=</sup> 0, 1, ··· (35)

)�<sup>2</sup> <sup>&</sup>lt; *tolres*�**F**(*T*0)�<sup>2</sup> (36)

**r**<sup>0</sup> = −**F**(*T*) − **J***δT*0. (37)

�**J***δTj* + **F**(*T*)�<sup>2</sup> < *γ*�**F**(*T*)�2, (39)

*th* Krylov iteration minimizes

**r**0, (38)

*�* , (40)

**3.3 The Jacobian-Free Newton Krylov method and forming the IMEX function**

**<sup>J</sup>**(*Tk*)*δT<sup>k</sup>* <sup>=</sup> <sup>−</sup>**F**(*T<sup>k</sup>*

�**F**(*T<sup>k</sup>*

at a fixed *k*. Let *j* be the Krylov iteration index. The *j*

{**r**0,**Jr**0,**J**2**r**0, ··· ,**J***j*−1**r**0} , and can be written as

following inexact Newton criteria (Dembo, 1982)

is terminated based on a required drop in the norm of the nonlinear residual, i.e,

iteratively over a sequence of linear system defined by

where **J**(*Tk*) = *<sup>∂</sup>***<sup>F</sup>**

approximated by

**Given** *T<sup>k</sup>* where *k* represents the current Newton iteration. **Call** Hydrodynamics block with (*ρn*, *un*, *En*, *Tk*) to compute *ρn*<sup>+</sup>1, *un*<sup>+</sup>1, *E*∗. **Form** *F*(*Tk*) based on the Crank-Nicolson method, *<sup>F</sup>*(*Tk*) = [*cvρn*+1*T<sup>k</sup>*<sup>+</sup> <sup>1</sup> <sup>2</sup> *<sup>ρ</sup>n*+<sup>1</sup>(*un*+<sup>1</sup>)<sup>2</sup>−*E*<sup>∗</sup> ] <sup>Δ</sup>*<sup>t</sup>* <sup>−</sup> <sup>1</sup> 2 *∂ <sup>∂</sup><sup>V</sup>* (*Aκ<sup>k</sup> <sup>∂</sup>T<sup>k</sup> <sup>∂</sup><sup>r</sup>* ) <sup>−</sup> <sup>1</sup> 2 *∂ <sup>∂</sup><sup>V</sup>* (*Aκ<sup>n</sup> <sup>∂</sup>T<sup>n</sup> <sup>∂</sup><sup>r</sup>* ).

It is important to note that we are not iterating between the implicit and explicit blocks. Instead we are executing the explicit block inside of a nonlinear function evaluation defined by *F*(*Tk*). The unique properties of JFNK allow us to perform a Newton iteration on this IMEX function, and thus JFNK is a required component of this nonlinearly converged IMEX approach.
