Machine learning moment closure models for the radiative transfer equation II: enforcing global hyperbolicity in gradient based closures

This is the second paper in a series in which we develop machine learning (ML) moment closure models for the radiative transfer equation (RTE). In our previous work \cite{huang2021gradient}, we proposed an approach to directly learn the gradient of the unclosed high order moment, which performs much better than learning the moment itself and the conventional $P_N$ closure. However, the ML moment closure model in \cite{huang2021gradient} is not able to guarantee hyperbolicity and long time stability. We propose in this paper a method to enforce the global hyperbolicity of the ML closure model. The main idea is to seek a symmetrizer (a symmetric positive definite matrix) for the closure system, and derive constraints such that the system is globally symmetrizable hyperbolic. It is shown that the new ML closure system inherits the dissipativeness of the RTE and preserves the correct diffusion limit as the Knunsden number goes to zero. Several benchmark tests including the Gaussian source problem and the two-material problem show the good accuracy, long time stability and generalizability of our globally hyperbolic ML closure model.


Introduction
In this paper we introduce an extension to our previous work on ML closures for radiative transfer modeling [20]. The new approach offers long time stability by ensuring mathematical consistency between the closure and the macroscopic model. Further, the numerical results in section 4 demonstrate the plausibility of capturing kinetic effects in a moment system with a handful of moments and an appropriate closure model.
The study of radiative transfer is of vital importance in many fields of science and engineering including astrophysics [35], heat transfer [24], and optical imaging [23]. The fundamental equation describing radiative transfer is an integro-differential equation, termed as the radiative transfer equation (RTE). Nevertheless, in most cases, any mesh based numerical discretization of the RTE equation leads to unacceptable computational costs due to the curse of dimensionality.
Moment methods study the evolution of a finite number of moments of the specific intensity in the RTE. Typically, the evolution of a lower-order moment depends on a higher-order moment, leading to what is known as the moment closure problem. Hence, one has to introduce suitable closure relations that relates the highest order moment included with the lower order moments in order to get a closed system of equations. Many moment closure models have been developed, including the P N model [8]; the variable Eddington factor models [28,33]; the entropy-based M N models [19,2,1]; the positive P N models [18]; the filtered P N models [32,25]; the B 2 models [3]; and the M P N model [12,13,29].
In moment closure problems, hyperbolicity is a critical issue, which is essential for a system of first-order partial differential equations (PDEs) to be well-posed [38]. The pioneering work on the moment closure for the gas kinetic theory by Grad in [15] is the most basic one among the moment models. However, it was discovered in [7] that the equilibrium is on the boundary of the hyperbolicity region for the Grad's 13-moment model in the three-dimensional case. Due to such a deficiency, the application of the moment method is severely restricted. This issue also attracts a lot of attention, with many papers in the literature focusing on the development of globally hyperbolic moment systems [6,12,13,29].
The traditional trade off in introducing a closure relation and solving a moment model instead of a kinetic equation is generic accuracy verses practical computability. However, thanks to the rapid development of machine learning (ML) and data-driven modeling [5,36,16], a new approach to solve the moment closure problem has emerged based on ML [17,37,21,4,30,40,31,20].
This approach offers a path for multi-scale problems that is relatively unique, promising to capture kinetic effects in a moment model with only a handful of moments. For the detailed literature review, we refer readers to [20]. We only remark that most of the works mentioned above are not able to guarantee hyperbolicity or long time stability, except the work in [21]. In [21], based on the conservation-dissipation formalism [43] of irreversible thermodynamics, the authors proposed a stable ML closure model with hyperbolicity and Galilean invariance for the Boltzmann BGK equation. Nevertheless, their model is limited to one extra non-equilibrium variable and it is still not clear how to generalize to arbitrary number of moments.
The work in this paper is a continuation of the previous work in [20], where we proposed to directly learn a closure that relates the gradient of the highest order moment to gradients of lower order moments. This new approach is consistent with the exact closure we derived for the free streaming limit and also provides a natural output normalization [20]. A variety of numerical tests show that the ML closure model in [20] has better accuracy than an ML closure based on learning a relation between the moments, as opposed to a relation between the gradients, and the conventional P N closure [20]. However, it is not able to guarantee hyperbolicity and long time simulations are not always satisfactory. Consequently, the focus of this work is to develop a method to enforce global hyperbolicity of the ML closure model. The main idea is motivated by the observation that the coefficient matrix for the P N closure is a tridiagonal matrix with positive off-diagonal entries. This indicates the existence of a diagonal symmetrizer matrix such that the P N closure is symmetrizable hyperbolic. Motivated by this observation, we propose to only keep the last several components in the ansatz, see equation (2.4) in Section 2.1, and seek a symmetrizer, which is a symmetric positive definite (SPD) matrix, of a block diagonal form. For degrees of freedom no larger than four in the ansatz, relating the gradient of the (N + 1) th moment to the gradient of the N th , (N − 1) th , (N − 2) th , and (N − 3) th moments, we derive constraints to guarantee the hyperbolicity of the ML closure model, see Theorem 2.1 in Section 2.1. Moreover, due to the block diagonal structure of the symmetrizer, we show that the ML closure system also satisfies the structural stability condition in [41] when taking into account the relaxation effects of the source terms, see Theorem 2.4 in Section 2.1. This condition characterises the dissipation feature of the moment closure system and is analogous to the H-theorem for the Boltzmann equation. We also consider our moment closure system under a diffusive scaling and formally show that our model preserves the correct diffusion limit as the Knunsden number goes to zero. The justification of the formal asymptotic analysis could be rigorously verified [34,27], with the aid of the structural stability condition. We also remark that, unlike [21], the work in this paper constructs ML closure models that are hyperbolic with an arbitrary number of moments. We numerically tested that the hyperbolic ML closure model has good accuracy in a variety of numerical examples. Moreover, the non-hyperbolic ML closure model in [20] blows up for long time simulations, while the hyperbolic one remains stable and has good accuracy.
The remainder of this paper is organized as follows. In Section 2, we introduce the hyperbolic ML moment closure model and discuss the diffusion limit of the model. In Section 3, we present the details in the training of the neural networks. The effectiveness of our ML closure model is demonstrated through extensive numerical results in Section 4. Some concluding remarks are given in Section 5.

ML moment closure for radiative transfer equation
In this section, we first review the ML moment closure method for the RTE in slab geometry proposed in [20] and propose our approach to enforce the hyperbolicity of the ML moment closure method. We then formally show that the resulting closure model can capture the correct diffusion limit.

Hyperbolic ML closure model
We consider the time-dependent RTE for a gray medium in slab geometry: Here, f = f (x, v, t) is the specific intensity of radiation. The variable v ∈ [−1, 1] is the cosine of the angle between the photon velocity and the x-axis. σ s = σ s (x) ≥ 0 and σ a = σ a (x) ≥ 0 are the scattering and absorption coefficients.
Denote the k-th order Legendre polynomial by P k = P k (x). Define the k-th order moment by Multiplying by P k (v) on both sides of (2.1) and integrating over v ∈ [−1, 1], we derive the moment equations: The learned gradient approach proposed in [20] is to find a relation between ∂ x m N +1 and the gradients on lower order moments: approximated by a neural network and learned from data. In this way, we obtain a quasi-linear first-order systems. This approach is shown to have uniform accuracy in the optically thick regime, intermediate regime and the optically thin regime. Moreover, the accuracy of this gradient-based model is much better than the approach based on creating a ML closure directly trained to match the moments, as well as the conventional P N closure. However, the learned gradient model exhibits numerical instability due to the loss of hyperbolicity [20]. This severely restricts the application of this model, especially for long time simulations.
The main idea of this work is to enforce the hyperbolicity by demanding that the coefficient matrix of the closure system A is real diagonizable. For this purpose, we seek a SPD matrix A 0 (also called a symmetrizer) such that A 0 A is a symmetric matrix. Namely, the system is symmetrizable hyperbolic, see the details in Appendix A. This places certain constraints on the functions N i .
Plugging (2.4) into the moment closure system, we have: Then, we can write down the coefficient matrix of this moment closure system: We first observe that, for the P N closure, i.e. when N j = 0 for 0 ≤ j ≤ N , A is a tridiagonal matrix with positive off-diagonal entries. In this case, A 0 can be taken as a diagonal matrix: and A 0 A is a tridiagonal symmetric matrix: Inspired by this observation, we propose to first investigate the case when the ansatz is given by only the last (k + 1) components where k is a parameter and k ≤ N . We will find a SPD matrix A 0 of the form With some algebraic calculation (details given in the appendix), when k = 3, we can derive the constraints explicitly as shown in the following theorem.
where a i has been specified in (2.6). If the coefficients a i for i = N − 3, N − 2, N − 1, N satisfy the following constraints: equivalent to the following constraints on N i : Proof. The proof is given in Appendix B.
Corollary 2.2. If we keep only 3 degrees of freedoms by setting a N −3 = 0, the constraint (2.12) reduces to 15) and the equivalent constraint (2.13) reduces to If we further set a N −3 = a N −2 = 0 by taking 2 degrees of freedoms, the constraint (2.12) reduces to a N −1 > 0, (2.17) and the equivalent constraint (2.13) reduces to In this case, the coefficient matrix A is a tridiagonal matrix with positive off-diagonal entries, which is real diagonalizable.
Remark 2.3. In principle, one could generalize the result in Theorem 2.1 to more than 4 degrees of freedom, by following the same lines in the proof given in Appendix B. However, the constraints for hyperbolicity will be a set of implicit inequalities. How to incorporate the implicit constraint into the architecture of the neural network would be an interesting topic to explore.
Besides the hyperbolicity for the convection term, one may also be interested in the relaxation effect of the source term. In this regard, the related property is the structural stability condition proposed in [41], which includes the constraints on the convection term, collision term, and the coupling of both. This stability condition is shown to be critical for the existence of the solutions [34], and satisfied by many moment closure systems [11]. In addition, the structural stability condition to the moment systems is same as H-theorem to the Boltzmann equation, which characterizes the dissipation property of the moment systems. For the convenience of the reader, we review the structural stability condition in Appendix C. It is easy to show that our new ML moment closure system also satisfied the structural stability condition: Theorem 2.4. For k = 3 and N ≥ 3, if the constraint (2.12) (or equivalently (2.13)) is satisfied, then the ML moment closure system with the closure relation (2.9) satisfies the structural stability condition in [41].
Proof. We only prove for the case σ a = 0. The other case can be proven by following the same line.
Denote Q = Q(U ) the source term of the closure system and Q U = ∂Q ∂U . Then we can compute Here, C = diag(0, σ s I N −3 ) and I r denote the identity matrix of order r. Notice that A 0 Q U is a symmetric negative semi-definite matrix. Then the structural stability condition can be easily checked for this moment closure system by taking P to be a scalar matrix in Appendix C.

Diffusion limit
In this part, we show formally that the hyperbolic ML moment closure model has the correct diffusion limit. Consider the RTE under a diffusive scaling: where ε > 0 denotes the Knudsen number, which is the ratio of mean free path of particles to the characteristic length. It is well known that as ε → 0, the kinetic transport equation (2.20) will converge to the diffusion limit-macroscopic model: where m 0 is the zeroth-order moment of f . For the RTE under a diffusive scaling (2.20), the corresponding moment closure equation is: Next, we formally show that the ML moment closure model (2.22) with N ≥ 3 converges to the diffusion limit (2.21) as ε → 0. We expand m k for k = 0, · · · , 2 as: and plug into the first three equations of (2.22). Collecting the O(ε) term of the first equation, we Collecting the O(1) term of the second equation, we have 1 . Combining the above three relations, we eventually reach: This indicates that the leading-order term of m 0 satisfies the diffusion equation (2.21).
To justify the above formal asymptotic analysis, it turns out that the structural stability condition is essential. Indeed, under the structural stability condition, one can rigorously prove that the solution of the ML moment closure model will converge to that of the macroscopic model (2.21) as ε goes to 0 [34,27], if the constraint (2.12) (or equivalently (2.13)) is satisfied and the neural network in (2.22) satisfies some regularity conditions. Although the regularity of the neural network seems not easy to be validated, we will show numerically in Section 4 that our hyperbolic ML closure model can capture the correct diffusion limit.

Training of the neural network
In this section, we provide details of the training of the proposed neural network that enforces the constraint (2.12) (or equivalently (2.13)) for the hyperbolicity for the case when k = 3. We design a fully-connected neural network denoted by M : R N +1 → R 4 with the input being the lower order moments (m 0 , m 1 , · · · , m N ) and the output being M = (M 1 , M 2 , M 3 , M 4 ). Then, we do the following post-processing to the output of the neural network: where h is the function defined in (2.14). Here σ : R → R is a positive function, i.e., σ(x) > 0 for any x ∈ R. We test several positive functions to enforce this constraint, including exponential function, ReLU (Rectified Linear Unit) function, softplus function and square function. Our numerical tests show that the neural network with σ to be the softplus function has the smallest L 2 error in the training data.
In the numerical implementation, we take the number of layers to be 6 and the number of nodes to be 256 in the fully-connected neural network and use the hyperbolic tangent activation function.
Other hyperparameters in the training are the same as those in [20]. The training data comes from numerically solving the RTE using the space-time discontinuous Galerkin (DG) method [9,10] with a range of initial conditions in the form of truncated Fourier series and different scattering and absorption coefficients which are constants over the computational domain. We train the neural network with 100 different initial data. For each initial data, we run the numerical solver up to t = 1. The other parameters are the same as in [20].

Numerical results
In this section, we show the performance of our ML closure model on a variety of benchmark tests, including problems with constant scattering and absorption coefficients, Gaussian source problems and two-material problems. The main focus of the tests is on the comparison of hyperbolic ML closure (termed as "hyperbolic closure") with k = 3, ML closure with learning gradient in [20] (termed as "non-hyperbolic closure") with k = N , and the classical P N closure [8].
In all numerical examples, we take the physical domain to be the unit interval [0, 1] and periodic boundary conditions are imposed. To numerically solve the moment closure system, we apply the fifth-order finite difference WENO scheme [22] with a Lax-Friedrichs flux splitting for the spatial discretization, and the third order strong-stability-preserving Runge-Kutta (RK) scheme [39] for the time discretization. We take the grid number in space to be N x = 256. For the hyperbolic ML closure model, the CFL condition is taken to be ∆t = 0.8∆x/c where c denote the maximum eigenvalues in all the grid points, and the penalty constant in the Lax-Friedrichs numerical flux α LF = c. For the non-hyperbolic ML closure model, we impose larger numerical diffusion by taking the penalty constant in the Lax-Friedrichs numerical flux α LF = 5 and fixing ∆t = 0.1∆x [20].  system (N = 6). It is observed that, at t = 0.5, the two ML moment closures agree well with the RTE, except that some minor oscillations appear in the non-hyperbolic ML closure. At t = 1, the solutions of non-hyperbolic ML moment closure blow up, while the hyperbolic one stays stable and has good agreement with the RTE. As a comparison, the P N closure has large deviations from the exact solution at both t = 0.5 and t = 1.
In Figure 4.2, we display the log-log scatter plots of the relative L 2 error versus the scattering coefficient for N = 6. We observe that, at t = 0.5, both ML closures have better accuracy than the P N closure. Moreover, in the optically thick regime, all the closures perform well. At t = 1, the numerical solutions of the non-hyperbolic ML moment closure blow up at some data points, while the hyperbolic one still has good accuracy.
In    Results obtained by the hyperbolic closure.
closure system and the solution generated by the RTE in the optically thin regime (σ s = σ a = 1).
We observe that the hyperbolic closure generates good predictions in the long time simulation up to t = 10. We also display the maximum eigenvalues in all the grid points during the time evolution in Figure 4.3 (b). It is observed that the eigenvalues are always real numbers, which validates the hyperbolicity feature of the closure system. However, the eigenvalues could be as large as four during the time evolution, which violates the physical characteristic speed of 1 in the RTE. This also results in smaller time step size in the simulation of the closure system. It is unclear how to preserve the physical characteristic speed in the moment closure system in this setting.
Example 4.2 (Gaussian source problem). In this example, we investigate the RTE with the initial condition to be Gaussian distribution in the physical domain, named Gaussian source problem in literature [14,12]: In this test, we take c 1 = 0.5, c 2 = 2.5, x 0 = 0.5 and θ = 0.01.
In Figure 4.4, we present the results obtained using various closure models. We observe good agreement between the two ML closure model and the kinetic model at t = 1, while the P N model has large deviations from the kinetic model. Moreover, the non-hyperbolic closure is more accurate than the hyperbolic closure, perhaps due to the fact that the non-hyperbolic closure has larger k values and thus can approximate the closure better. The solution of the non-hyperbolic ML closure model blows up at t = 2, while the hyperbolic ML model stays stable and is more accurate than the P N closure. Although the deviation of the hyperbolic ML model "looks" large in  respectively. This again illustrates that the hyperbolicity is an essential property in the ML closure model. Example 4.3 (two-material problem). In this example, we consider the two-material problem [26].
The numerical results are shown in Figure 4.5. The gray part is in the optically thin regime and the other part is in the intermediate regime. We observe that the two ML closure models are more accurate than the P N closure at t = 0.5 and the hyperbolic ML closure is slightly better than the non-hyperbolic one. At t = 1, the hyperbolic ML closure still agrees well with the kinetic model.
Large deviations exist for the P N closure in the optically thin regime and severe oscillations occur for the non-hyperbolic ML closure.

Concluding remarks
In this paper, we propose a method to enforce the global hyperbolicity of the ML closure model.
We find a symmetrizer (a symmetric positive definite matrix) for the closure system, and derive constraints that guarantee the system is globally symmetrizable hyperbolic. Moreover, we show that the closure system also inherits the dissipativeness of the RTE by checking the structural stability condition. A variety of benchmark tests including Gaussian source problem and the two-material problem show good accuracy, correct diffusion limit and generalization ability of our ML closure model. The new approach also offers long time stability by ensuring mathematical consistency between the closure system and the macroscopic model. Further, the new method demonstrates the plausibility of capturing kinetic effects in a moment system with a handful of moments and an appropriate closure model.
There are several issues that are worthy of further investigations. First, one could in principle generalize the result in Theorem 2.1 to more than 4 degrees of freedom, by following the same lines in the proof given in Appendix B. In this general case, the hyperbolicity constraints will be a set of implicit inequalities. How to incorporate the constraint of implicit inequalities into the architecture of the neural network would be an interesting topic to explore. Moreover, the moment closure model is expected to have better accuracy with more degrees of freedom. Another topic is that the characteristic speed of the current model may exceed the physical bound, as we observed in the numerical tests. This unphysical characteristic speed results in a small time step size in the computation, and is something we would like to address in our future work. Another interesting topic is to extend the current approach to the two dimensional case. These issues constitute the body of our ongoing work.

Acknowledgment
We thank Michael M. Crockatt from Sandia National Laboratories for providing a numerical solver for the radiative transfer equations. We also would like to acknowledge the High Performance Computing Center (HPCC) at Michigan State University for providing computational resources that have contributed to the research results reported within this paper.

Appendix A Hyperbolicity
In this part, we review the definition of the hyperbolicity and some equivalent conditions. See also the details in Chapter 3 in [38].
Consider the first-order system of equations in 1D: The following classical conclusion holds for these two definitions: Proposition A.3. The symmetrizable hyperbolicity is a necessary and sufficient condition for the first-order system (A.1) to be hyperbolic.
Proof. We start by the proof of sufficiency. By the definition of symmetrizable hyperbolic, there exists a symmetric positive definite (SPD) matrix A 0 such that A 0 A are symmetric. Since A 0 is a SPD matrix, there exists an invertible symmetric matrix B satisfying A 0 = B 2 . Then, we compute Thus, BAB −1 is symmetric since A 0 A and B are symmetric. Then we have A is real diagonlizable.
On the other hand, suppose that A is real diagonlizable, there exist an invertible real matrix P such that A = P −1 DP with D a real diagnoal matrix. Take A 0 = P T P . Then A 0 is a SPD matrix and A 0 A = P T DP is symmetric.

Appendix B Proof of Theorem 2.1
Proof. To save space, we only present the proof in the case of a N −3 = 0. One also can prove the case of a N −3 = 0 by following the same line.
We first write A in (2.11) into a block matrix: . Then we compute A 0 A: Here, 0 m×n denote the zero matrix of size m × n. It is easy to see that DA 1 is symmetric.
Next, we compute other blocks in A 0 A: Then, we compute (B.6) For A 0 A to be symmetric, we need BA 4 is symmetric and DA 2 = (BA 3 ) T . Thus, we have the (B.7) We solve for b 11 , b 12 and b 22 from the above linear system: This is just the condition (2.12) with a N −3 = 0 (or equivalently (2.15)). Finally, (2.13) can be easily derived, and the proof is omitted for brevity.

Appendix C Structural stability condition
In this part, we review the structural stability condition proposed in [41], which is fundamental for the quasilinear first-order hyperbolic system with source terms: The stability condition in [41] reads as: (i). There exists an invertible n × n matrix P (U ) and an invertible r × r matrix S(U ), defined on the equilibrium manifold E, such that Here Q U = ∂Q ∂U and I r is the identity matrix of order r.
This set of conditions has been tacitly respected by many well-developed physical models [42].
The first condition is classical for initial value problems of the system of ordinary differential equations (ODEs), while the second one means the symmetrizable hyperbolicity of the PDE system.
The third condition characterizes a kind of coupling between the convection term and the source term. This set of conditions implies the existence and stability of the zero relaxation limit of the corresponding initial value problems [41].