Sparse Representation Based Efficient Radiation Symmetry Analysis Method for Cylindrical Model of Inertial Confinement Fusion

TL;DR

This paper introduces a sparse representation method for efficient radiation symmetry analysis in inertial confinement fusion models, significantly reducing computation time while maintaining accuracy.

Contribution

The paper proposes a novel sparse representation approach using specific polynomial bases and a new algorithm to efficiently solve nonlinear radiation flux equations in ICF models.

Findings

Achieves tenfold reduction in computation time for radiation flux calculation.

Maintains high accuracy with significantly fewer equations.

Efficiency increases as mesh element size decreases.

Abstract

Radiation symmetry evaluation is critical to the laser driven Inertial Confinement Fusion (ICF), which is usually done by solving a view-factor equation model. The model is nonlinear, and the number of equations can be very large when the size of discrete mesh element is very small to achieve a prescribed accuracy, which may lead to an intensive equation solving process. In this paper, an efficient radiation symmetry analysis approach based on sparse representation is presented, in which, 1) the Spherical harmonics, annular Zernike polynomials and Legendre-Fourier polynomials are employed to sparsely represent the radiation flux on the capsule and cylindrical cavity, and the nonlinear energy equilibrium equations are transformed into the equations with sparse coefficients, which means there are many redundant equations, 2) only a few equations are selected to recover such sparse coefficients with Latin hypercube sampling, 3) a Conjugate Gradient Subspace Thresholding Pursuit (CGSTP) algorithm is then given to rapidly obtain such sparse coefficients equation with as few iterations as possible. Finally, the proposed method is validated with two experiment targets for Shenguang II and Shenguang III laser facility in China. The results show that only one tenth of computation time is required to solve one tenth of equations to achieve the radiation flux with comparable accuracy. Further more, the solution is much more efficient as the size of discrete mesh element decreases, in which, only 1.2\% computation time is required to obtain the accurate result.

Figures

Tables

Table 1. TABLE I: The first 9 expansion terms of the LF polynomials sequence

Terms	Order	Polynomials	Expression
$W_0$	0	$L_0{F}_0$	1
$W_1$	1	$L_0{F}_{-1}$	$\sin k_w{\phi }_w$
$W_2$	1	$L_0{F}_1$	$\cos k_w{\phi }_w$
$W_3$	1	$L_1{F}_0$	$z$
$W_4$	2	$L_0{F}_{-2}$	$\sin 2k_w{\phi }_w$
$W_5$	2	$L_0{F}_2$	$\cos 2k_w{\phi }_w$
$W_6$	2	$L_1{F}_{-1}$	$z\sin k_w{\phi }_w$
$W_7$	2	$L_1{F}_1$	$z\cos k_w{\phi }_w$
$W_8$	2	$L_2{F}_0$	$\left(3z^2-1\right)/2$

Table 2. TABLE II: The size and number of discrete elements of four simulation models. ”H”, ”C”, ” Δ ​ r Δ 𝑟 \Delta r ” and ” Δ ​ h Δ ℎ \Delta h ” represent the minimum dimensions of the longitude direction of capsule, the latitudinal direction of capsule, the radial direction of cavity, and the axial direction of cavity, respectively.

Model	Capsule	The end face of cavity	The cylinder surface	Total number
S2-1	2592 (H:$5^{\circ }$, C:$5^{\circ }$)	2016$\times$2 (H:${2.5}^{\circ }$, $\Delta r$:15um)	3152 (H:$5^{\circ }$, $\Delta h$:30um)	9776
S2-2	10368 (H:${2.5}^{\circ }$, C:${2.5}^{\circ }$)	8064$\times$2 (H:${1.25}^{\circ }$, $\Delta r$:7.5um)	12456 (H:${2.5}^{\circ }$, $\Delta h$:15um)	38952
S3-1	2592 (H:$5^{\circ }$, C:$5^{\circ }$)	4320 $\times$2(H:${2.5}^{\circ }$, $\Delta r$:15um)	9504 (H:$5^{\circ }$, $\Delta h$:30um)	20736
S3-2	10368 (H:${2.5}^{\circ }$, C:${2.5}^{\circ }$)	17280$\times$2 (H:${1.25}^{\circ }$, $\Delta r$:7.5um)	38016 (H:${2.5}^{\circ }$, $\Delta h$:15um)	82944

Table 3. TABLE III: Number of samples of four simulation models

Model	Capsule	The end face of cavity	The cylinder surface	Total sampling number	Sampling rate
S2-1	150	300	400	850	8.7%
S2-2	150	300	450	900	2.3%
S3-1	150	300	350	800	3.9%
S3-2	150	400	400	950	1.1%

Table 4. TABLE IV: Running time of conventional method (unit: second),”–” denotes out of memory

		Newton-Raphson method		PCG method
Model	View-factor calculation time	Iteration time	Total time	Iteration time	Total time
S2-1	4341.5	27.3	4368.8	21.4	4362.9
S2-2	69045.8	2203.6	71249.4	313.2	69359.1
S3-1	16253.1	1288.4	17541.5	96.8	16349.9
S3-2	255843.2	—	—	1572.6	257415.8

Table 5. TABLE V: Running time of Compressed Sensing method (unit: second).

			IHT		CGIHT		CGSTP
Model	View-factor calculation time	Sparse basis calculation time	Iteration time	Total time	Iteration time	Total time	Iteration time	Total time
S2-1	407.1	11.8	332.4	751.3	167.9	586.8	13.7	432.6
S2-2	1518.4	40.2	1142.6	2701.2	683.1	2241.7	38.6	1597.2
S3-1	658.6	22.7	671.8	1353.1	322.4	1003.7	20.7	705.0
S3-2	3014.2	103.5	2248.3	5366.0	1131.2	4248.9	84.5	3202.2

Table 6. TABLE VI: Number of samples of four simulation models

Model	Newton-Raphson	PCG	IHT	CGIHT	CGSTP
S2-1	3	10	1273	841	38
S2-2	2	8	1124	729	35
S3-1	3	9	837	657	31
S3-2	2	8	811	594	27

Table 7. TABLE VII: Reconstruction RMSE of different methods

Model	NR	PCG	IHT	CGIHT	CGSTP
S2-1	$3.14\times {10}^{-10}$	$5.52\times {10}^{-6}$	$3.71\times {10}^{-3}$	$1.53\times {10}^{-3}$	$7.84\times {10}^{-4}$
S2-2	$2.33\times {10}^{-10}$	$3.91\times {10}^{-6}$	$1.27\times {10}^{-3}$	$9.84\times {10}^{-4}$	$6.20\times {10}^{-4}$
S3-1	$1.65\times {10}^{-10}$	$4.13\times {10}^{-6}$	$9.81\times {10}^{-4}$	$8.92\times {10}^{-4}$	$5.75\times {10}^{-4}$
S3-2	$1.08\times {10}^{-10}$	$1.64\times {10}^{-6}$	$9.55\times {10}^{-4}$	$8.51\times {10}^{-4}$	$5.24\times {10}^{-4}$

Equations

$${\lambda }_i(E_i+\sum _{j=1}^N{F}_{ji}{\zeta }_j{B}_j)=B_i,$$

$$F_{ji}={\frac{1}{\pi }\left[{\mathbf{n}}_i\cdot \left({\mathbf{p}}_j-{\mathbf{p}}_i\right)\right]\cdot \left[{\mathbf{n}}_j\cdot \left({\mathbf{p}}_i-{\mathbf{p}}_j\right)\right]/\Vert {\mathbf{p}}_i-{\mathbf{p}}_j\Vert }^4,$$

$${\lambda }_i=\frac{1}{1+{\upsilon }^{-{\beta }^{-1}}{B}_i^{{\beta }^{-1}-1}{t}^{-\alpha {\beta }^{-1}}},$$

$$(\mathbf{I}-\mathbf{V})\mathbf{B}+C\mathbf{B}{{}^{\circ }}^{(1/\beta )}=\mathbf{E}$$

$$\mathbf{f}=\sum _{i=0}^N{c}_i{\psi }_i=\mathbf{\Psi }\mathbf{c}$$

$$\mathbf{y}=\mathbf{\Phi }\mathbf{f}=\mathbf{\Phi }\mathbf{\Psi }\mathbf{c}=\mathbf{Ac}$$

$$(1-{\delta }_K)||\mathbf{c}|{|}_2^2\le ||\mathbf{\Psi }\mathbf{c}|{|}_2^2\le (1+{\delta }_K)||\mathbf{c}|{|}_2^2.$$

$$\arg \underset{\mathbf{c}}{\min }{\Vert \mathbf{y}-\mathbf{Ac}\Vert }_{2}s.t.{\Vert \mathbf{c}\Vert }_0=K$$

$$B_s\left({\theta }_s^{(i)},{\phi }_s^{(i)}\right)=\sum _{n_s=0}^{+\infty }{c}_{n_s}{Y}_{\left(n_s\right)}\left({\theta }_s^{(i)},{\phi }_s^{(i)}\right),$$

$$Y_{{m_{S}}}^{{k_{S}}}({\theta_{s}},{\varphi_{s}})=N_{{m_{S}}}^{{k_{S}}}P_{{m_{S}}}^{{k_{S}}}(\cos{\theta_{s}})\left\{{\begin{array}[]{*{20}{c}}{\cos{k_{s}}{\varphi_{s}}{\rm{}}\quad if\quad{\rm{}}{k_{S}}\geq 0}\\ {\sin{k_{s}}{\varphi_{s}}{\rm{}}\quad if\quad{\rm{}}{k_{S}}<0}\\ \end{array}}\right.{\rm{,}}\quad(-{m_{S}}\leq{k_{S}}\leq{m_{S}})$$

$$B_S\left({\theta }_s,{\phi }_s\right)=\sum _{j=1}^{N_s}\sum _{n_s=0}^{+\infty }{c}_{n_s}{Y}_{n_s}\left({\theta }_s^{(j)},{\phi }_s^{(j)}\right)\approx \sum _{j=1}^{N_s}\sum _{n_s=0}^{L_s}{c}_{n_s}{Y}_{n_s}\left({\theta }_s^{(j)},{\phi }_s^{(j)}\right)=\mathbf{Y}{\mathbf{c}}^s,$$

$${{\bm{Y}}_{{N_{S}}\times{L_{S}}}}=\left[{\begin{array}[]{*{20}{c}}{{Y_{0}}(\theta_{s}^{(1)},\varphi_{s}^{(1)})}&{{Y_{1}}(\theta_{s}^{(1)},\varphi_{s}^{(1)})}&\cdots&{{Y_{{L_{S}}}}(\theta_{s}^{(1)},\varphi_{s}^{(1)})}\\ {{Y_{0}}(\theta_{s}^{(1)},\varphi_{s}^{(2)})}&{{Y_{1}}(\theta_{s}^{(1)},\varphi_{s}^{(2)})}&\cdots&{{Y_{{L_{S}}}}(\theta_{s}^{(1)},\varphi_{s}^{(2)})}\\ \vdots&\vdots&\ddots&\vdots\\ {{Y_{0}}(\theta_{s}^{({N_{\theta}})},\varphi_{s}^{({N_{{\varphi_{s}}}})})}&{{Y_{1}}(\theta_{s}^{({N_{\theta}})},\varphi_{s}^{({N_{{\varphi_{s}}}})})}&\cdots&{{Y_{{L_{S}}}}(\theta_{s}^{({N_{\theta}})},\varphi_{s}^{({N_{{\varphi_{s}}}})})}\\ \end{array}}\right],$$

$$U_{l_{u}}^{k_{u}}\left(r,\tilde{r};\varphi_{u}\right)=\left\{\begin{array}[]{ll}{R_{lu}^{k_{u}}(r,\tilde{r})\cos k_{u}\varphi_{u}}&{\text{ for }\quad k_{u}\geq 0}\\ {R_{l_{u}}^{k_{u}}(r,\tilde{r})\sin k_{u}\varphi_{u}}&{\text{ for }\quad k_{u}<0}\end{array}\right.$$

$$\underset{\tilde{r}}{\overset{1}{\int }}\underset{0}{\overset{2\pi }{\int }}{U}_{l_u}^{k_u}(r,\tilde{r};{\phi }_u)U_{l{}_u^{\prime }}^{k{}_u^{\prime }}(r,\tilde{r};{\phi }_u)rdrd{\phi }_u/\underset{\epsilon }{\overset{1}{\int }}\underset{0}{\overset{2\pi }{\int }}rdrd{\phi }_u={\delta }_{l_{u}l{}_u^{\prime }}{\delta }_{k_{u}k{}_u^{\prime }},$$

$$R_{2j}^0(r,\tilde{r})=R_{2j}^0{\left(\frac{r^2-{\tilde{r}}^2}{1-{\tilde{r}}^2}\right)}^{1/2}=P_j\left[\frac{2\left(r^2-{\tilde{r}}^2\right)}{1-{\tilde{r}}^2}-1\right],$$

$$R_{2j+k_n}^{k_u}(r,\tilde{r})={\left[\frac{1-{\tilde{r}}^2}{2\left(2j+k_u+1\right)H_j^{k_u}}\right]}^{1/2}{r}^{k_i}{Q}_j^{k_u}\left(r^2\right),$$

$$Q_j^{k_u}\left(r^2\right)=\frac{2\left(2j+2k_u-1\right)}{\left(j+k_u\right)\left(1-{\tilde{r}}^2\right)}\frac{H_j^{k_u-1}}{Q_j^{k_u-1}(0)}\sum _{i=0}^j\frac{Q_i^{k_u-1}(0)Q_i^{k_u-1}\left(r^2\right)}{H_i^{k_u-1}}$$

$$H_j^{k_u}=-\frac{2(2j+2k_u-1)}{(j+k_u)(1-{\tilde{r}}^2)}\frac{Q_{j+1}^{k_u-1}(0)}{Q_j^{k_u-1}(0)}{H}_j^{k_u-1}$$

$$Q_j^0(r^2)=R_{2j}^0(r,\tilde{r});H_j^0=\frac{1-{\tilde{r}}^2}{2(2j+1)}.$$

$$R_{2j+k_u}^{k_u}(r,\tilde{r})=R_{2j+k_u}^{|k_u|}(r,\tilde{r}).$$

$$R_{k_u}^{k_u}(r,\tilde{r})=\frac{r^{k_u}}{{(\sum _{i=0}^{k_u}{r}^{2i})}^{1/2}}.$$

$$B_u\left(r,{\phi }_u\right)=\sum _{i=1}^{N_u}\sum _{n_u=0}^{+\infty }{c}_{n_u}{U}_{n_u}\left(r^{(i)},\tilde{r};{\phi }_u^{(i)}\right)\approx \sum _{j=1}^{N_u}\sum _{n_u=0}^{L_u}{c}_{n_u}{U}_{n_u}\left(r^{(i)},\tilde{r};{\phi }_u^{(i)}\right)=\mathbf{U}{\mathbf{c}}^u,$$

$${{\bm{U}}_{{N_{u}}\times{L_{u}}}}=\left[{\begin{array}[]{*{20}{c}}{{U_{0}}({r^{(1)}},\tilde{r};\varphi_{u}^{(1)})}&{{U_{1}}({r^{(1)}},\tilde{r};\varphi_{u}^{(1)})}&\cdots&{{U_{{L_{u}}}}({r^{(1)}},\tilde{r};\varphi_{u}^{(1)})}\\ {{U_{0}}({r^{(1)}},\tilde{r};\varphi_{u}^{(2)})}&{{U_{1}}({r^{(1)}},\tilde{r};\varphi_{u}^{(2)})}&\cdots&{{U_{{L_{u}}}}({r^{(1)}},\tilde{r};\varphi_{u}^{(2)})}\\ \vdots&\vdots&\ddots&\vdots\\ {{U_{0}}({r^{({N_{r}})}},\tilde{r};\varphi_{u}^{({N_{{\varphi_{u}}}})})}&{{U_{1}}({r^{({N_{r}})}},\tilde{r};\varphi_{u}^{({N_{{\varphi_{u}}}})})}&\cdots&{{U_{{L_{u}}}}({r^{({N_{r}})}},\tilde{r};\varphi_{u}^{({N_{{\varphi_{u}}}})})}\\ \end{array}}\right]$$

$$F_{k_{w}}\left(\varphi_{w}\right)=\left\{\begin{array}[]{ll}{\cos k_{w}\varphi_{w}}&{k_{w}\geq 0}\\ {\sin k_{w}\varphi_{w}}&{k_{w}<0}\end{array}\right.$$

$$W_{n_w}\left(z,{\phi }_w\right)=P_{l_w}(z)F_{k_w}\left({\phi }_w\right),$$

$$B_w(z,{\phi }_w)=\sum _{i=1}^{N_w}\sum _{n_w=0}^{+\infty }{c}_{n_w}{W}_{n_w}(z^{(i)},{\phi }_w^{(i)})\approx \sum _{j=1}^{N_w}\sum _{n_w=0}^{L_w}{c}_{n_w}{W}_{n_w}(z^{(i)},{\phi }_w^{(i)})=\mathbf{W}{\mathbf{c}}^w,$$

$${{\bm{W}}_{{N_{w}}\times{L_{w}}}}=\left[{\begin{array}[]{*{20}{c}}{{W_{0}}({z^{(1)}},\varphi_{w}^{(1)})}&{{W_{1}}({z^{(1)}},\varphi_{w}^{(1)})}&\cdots&{{W_{{L_{w}}}}({z^{(1)}},\varphi_{w}^{(1)})}\\ {{W_{0}}({z^{(1)}},\varphi_{w}^{(2)})}&{{W_{1}}({z^{(1)}},\varphi_{w}^{(2)})}&\cdots&{{W_{{L_{w}}}}({z^{(1)}},\varphi_{w}^{(2)})}\\ \vdots&\vdots&\ddots&\vdots\\ {{W_{0}}({z^{({N_{z}})}},\varphi_{w}^{({N_{{\varphi_{w}}}})})}&{{W_{1}}({z^{({N_{z}})}},\varphi_{w}^{({N_{{\varphi_{w}}}})})}&\cdots&{{W_{{L_{w}}}}({z^{({N_{z}})}},\varphi_{w}^{({N_{{\varphi_{w}}}})})}\\ \end{array}}\right],$$

$$B=B_s+B_u+B_w={\mathbf{Y}}_{N_s\times L_s}{\mathbf{c}}^s+{\mathbf{U}}_{N_u\times L_u}{\mathbf{c}}^u+{\mathbf{W}}_{N_w\times L_w}{\mathbf{c}}^w,$$

$${\bm{B}}=\left[{\begin{array}[]{*{20}{c}}{{{\bm{Y}}_{{N_{s}}\times{L_{s}}}}}&{}\hfil&{}\hfil\\ {}\hfil&{{{\bm{U}}_{{N_{u}}\times{L_{u}}}}}&{}\hfil\\ {}\hfil&{}\hfil&{{{\bm{W}}_{{N_{w}}\times{L_{w}}}}}\\ \end{array}}\right]\times\left[{\begin{array}[]{*{20}{c}}{{{\bm{c}}^{s}}}\\ {{{\bm{c}}^{u}}}\\ {{{\bm{c}}^{w}}}\\ \end{array}}\right]={\bm{\Psi c}}.$$

$$f(\mathbf{c})=(\mathbf{I}-\mathbf{V})(\mathbf{\Psi }\mathbf{c})+C(\mathbf{\Psi }\mathbf{c}){\circ }^{1/\beta }-\mathbf{E}$$

Full text

I Introduction

Controllable nuclear fusion is potential to solve the energy crisis in the future, and laser driven Inertial Confinement Fusion (ICF) is supposed to be one of the most promising ways [1]. To achieve this goal, the driven asymmetry on the centrally located capsule needs to be evaluated and remains no more than 1% during the fusion process. Due to limited laser beams for NIF laser facility, i.e., 2 laser entrance holes and a few discrete diagnose holes, the radiation asymmetry may be much larger than such prescribed threshold [2]. Therefore, we need to efficiently compute the driven flux reached the capsule to evaluate its driven asymmetry. In the ICF, the radiation flux on the fuel capsule is related to the laser-plasma interactions (LPI) and transport of X-rays from the cavity wall to the capsule, which involves the solving of complex kinetic and hydrodynamic equations implemented in the codes such as LASNEX [3]. In practice, simple mathematical models such as view-factor codes are usually employed to compute the radiation flux distributed on the capsule, especially for the preliminary design and optimization of thermonuclear target structure and shape [3, 4]. As described in [1], the radiation flux is usually obtained by solving a non-linear view-factor based equation model, which can be solved by utilizing Newton-Rapshon [5], Jacobi iteration [6], Cholesky decomposition [7], or Preconditioned Conjugate Gradient method [8]. However, in order to improve the evaluation accuracy, the size of discrete element is usually very small, the number of elements or equations will increase significantly, which may lead to much time required to solve large scale nonlinear equations. The computation time may be unacceptable for researchers. Therefore, a new efficient computation approach is essentially required to significantly improve the efficiency of radiation symmetry evaluation and optimization.

Compressed sensing (CS), proposed by Donoho and Candes [9, 10], is a new method to reconstruct signals from significantly fewer measurements than that required by traditional methods, which has attracted considerable attention and achieves successful applications such as Medical imaging (MI) [11], Analog-digital Conversion [12], Computational Biology [13], and Computer Graphics [14]. In the field of radiation symmetry evaluation of ICF, compressed analysis approach has been applied to efficiently obtain the radiation flux distribution on the capsule [15]. However, the radiation flux computation model is simplified, and the accuracy of the radiation asymmetry evaluation is limited. The non-linear Time Dependent Energy Balance Model (TDEBM) presented in [16] is often used to compute the radiation flux. Therefore, compressive sensing approach for non-linear equation solving approach, such as Iterative Hard Thresholding (IHT) [17], can be applied. However, the step size is often taken as a constant, which may lead to more times of iteration convergence. Normalized Iterative Hard Thresholding (NIHT) [18] is further proposed, in which, fixed step length is replaced with a descending factor to accelerate the convergence rate. Nevertheless, the gradient based search direction in NIHT, may lead to the convergence is very slow when it approaches the minimum. Conjugate Gradient Iterative Hard Thresholding (CGIHT) is proposed in [19], to replace the search direction with the conjugate gradient direction, which can benefit accelerate the convergence. However, the search direction of CGIHT may not be conjugate, which may lead to more times of iterations being required. As one of the iterative greedy algorithms, Subspace Pursuit (SP) algorithm has attracted much attention for its backtracking idea in solving compressed sensing reconstruction problems, which requires only a small number of iterations [20]. However, SP algorithm can not be directly used to solve nonlinear problems. In this paper, we absorb the advantages of CGIHT algorithm and SP algorithm, and propose a Conjugate Gradient Subspace Thresholding Pursuit (CGSTP) algorithm to largely reduce the TDEBM and efficiently obtain the radiation flux. The core ideas include:

1. (1)

The spherical harmonics, annular Zernike polynomials and Legendre-Fourier polynomials are employed to sparsely represent the radiation flux on the capsule and spherical cavity, which enable only a few equations are required to recover the radiation flux with high accuracy, and significantly reduce the radiation computation model. 2. (2)

Over the three sets of polynomials domain, such reduced nonlinear equation model with sparse coefficients is formulated to enable compressed sensing algorithms, such as CGIHT, be used to efficiently obtain the radiation flux. 3. (3)

Conjugate Gradient Subspace Thresholding Pursuit algorithm is then presented to enable the adjacent search direction be conjugate, which facilitate the convergence of iteration, rapidly obtain the solution, and then efficiently evaluate the radiation symmetry for ICF experiments design.

The outline of this paper as follows. Section 2 introduces the nonlinear view-factor based equation model in a cylindrical cavity model. Radiation sparse representation model constructed with three orthogonal bases is described in section 3. Then section 4 gives a description of the proposed reconstruction algorithm, i.e. Conjugate Gradient Subspace Thresholding Pursuit (CGSTP) algorithm. At last, section 6 concludes following experiments validation in section 5.

II Sparse representation of radiation flux

II-A Radiation flux computation model based on view-factor

As shown in Figure 1 , eight laser beams are injected into the cylindrical cavity through two entrance holes. By interacting with the high-Z material on the inner wall of the cavity, the laser incident on the cavity is converted into X-ray, the X-ray irradiates evenly and drives the centrally located capsule implosion. With generated discrete mesh elements, the radiation flux on the wall of the cavity and surface of the capsule can be computed by using the view-factor function described in [15], i.e. the TDEBM. For a discrete mesh element, according to the energy conservation law, the energy received from other elements equals the energy emits to other elements. It can be formulated as:

[TABLE]

where $N$ is the total number of mesh elements, ${E_{i}}$ is the radiation flux comes directly from primary source converted from laser on the $i$-th mesh element which have a centroid ${{\bm{p}}_{i}}$ and normal ${{\bm{n}}_{i}}$, can be computed as ${E_{i}}=\sum\nolimits_{i=0}^{n}{S_{j}^{0}}{F_{ji}}{\zeta_{j}}$, here, $S_{j}^{0}$ is the radiation flux converted directly from laser beams of the $j$-th element. ${\zeta_{j}}$ is the area of the $j$-th mesh element with a centroid ${{\bm{p}}_{j}}$ and normal ${{\bm{n}}_{j}}$. ${F_{ji}}$ is the view-factor between the $j$-th mesh element and the $i$-th mesh element, which can be computed by

[TABLE]

and ${\lambda_{i}}$ is the albedo of the $i$-th mesh element, which is relating to the wall material, radiation time $t$ and the emit energy ${B_{i}}$, which can be given as

[TABLE]

where $\upsilon$ , $\alpha$ and $\beta$ are the material parameters in which $\upsilon=4.87$, $\alpha=8/13$, $\beta=16/13$. ${B_{i}}$ and ${B_{j}}$ are the radiation flux of the $i$-th mesh element and the $j$-th mesh element respectively. Let ${\bm{E}}={\left[{{E_{1}},{E_{2}},\cdot\cdot\cdot,{E_{n}}}\right]^{\rm T}}$ , ${\bm{S}}={\left[{S_{1}^{0},S_{2}^{0},\cdot\cdot\cdot,S_{N}^{0}}\right]^{\rm T}}$, ${\bm{F}}=({F_{ji}})$, ${\bm{\zeta}}={\left[{{\zeta_{1}},{\zeta_{2}},\cdot\cdot\cdot,{\zeta_{N}}}\right]^{\rm T}}$, ${\bm{B}}={\left[{{B_{1}},{B_{2}},\cdot\cdot\cdot,{B_{N}}}\right]^{\rm T}}$, $C={\upsilon^{-{\beta^{-1}}}}{t^{-\alpha{\beta^{-1}}}}$, where ${\rm T}$ indicates the matrix transpose. Let ${\bm{V}}=({F_{ji}}{\zeta_{j}})$ be an intermediate variable, then the matrix form of ${E_{i}}=\sum\nolimits_{i=0}^{n}{S_{j}^{0}}{F_{ji}}{\zeta_{j}}$ can be written as ${\bm{E}}={\bm{VS}}$ and Equation (3) can be written in a compact form as

[TABLE]

where $\bm{I}$ denotes an unit matrix. In Equation (4), $\bm{V}$ and $\bm{E}$ are known, $\bm{B}$ is the unique unknown variable. ∘ denote the hadamard power which means for any vector $\bm{{\cal V}}$ and constant $a$, $\bm{{\cal V}}{{{}^{\circ}}^{a}}={\left[{{\cal V}_{1}^{a},{\cal V}_{2}^{a},\cdot\cdot\cdot,{\cal V}_{n}^{a}}\right]^{\rm T}}$ holds. Since $\beta$ do not equal one, Equation (4) is nonlinear. Such nonlinear equations can be solved by utilizing some conventional iterative methods such as Newton-Raphson method. However, these methods is time consuming even unusable when the number of the discrete elements reaches ${10^{7}}$ or more.

II-B Compressed Sensing Method

Compressed Sensing (CS) also called Compressed Sampling, is an emerging technique, which has shown that sparse or compressible signal can be reconstructed from far fewer measurements than what is considered by Nyquist theorem [9]. Compressed sensing theory usually consists of sparse representation, measurement matrix, and sparse coefficient reconstruction.

II-B1 Sparse representation

Consider a $N\times 1$ signal ${\bm{f}}={\left[{f_{1},f_{2},\cdots,f_{N}}\right]^{\rm T}}$, which can be expressed as a linear combination of a complete basis ${\bm{\Psi}}$:

[TABLE]

where ${\bm{c}}={\left[{{c_{1}},{c_{2}},\cdot\cdot\cdot,{c_{N}}}\right]^{\rm T}}$, ${{\bm{\Psi}}}=\left[{{\bm{\psi}_{1}},{\bm{\psi}_{2}},\cdot\cdot\cdot,{\bm{\psi}_{N}}}\right]$. If there are $K$ non-zero elements in $\bm{c}$, then $\bm{c}$ is called $K$-sparse and we say $\bm{f}$ is sparse in the ${\bm{\Psi}}$ domain.

II-B2 Measurement matrix

$M$ projections of the signal $\bm{f}$ constitute a $M\times N(M\ll N)$ matrix ${\bm{\Phi}}$ called measurement matrix, then the sampling process can be described as:

[TABLE]

where ${\bm{A}}={\bm{\Phi}}{{\bm{\Psi}}}$ is known as measurement matrix. It is clear that the process of recovering $\bm{c}$ from $\bm{y}$ is an ill-posed problem because $M\ll N$ and seems impossible to solve such underdetermined linear equation. However, if the sparse coefficients $\bm{c}$ is $K$-sparse, and $\bm{\Phi}$ satisfy the Restricted Isometry Property (RIP) criterion [21], then $\bm{c}$ can be accurately recovering with a high probability from $M$ measurements. The RIP criterion can be expressed as

[TABLE]

The restricted isometry constant ${\delta_{K}}(0<{\delta_{K}}<1)$ is defined as the smallest constant for which this property holds for all $K$-sparse vectors $\bm{c}$.

II-B3 Sparse coefficients reconstruction algorithm

As can be seen from above, the sparse coefficient vector $\bm{c}$ can be recovered with high probability from $\bm{y}$ by solving the following ${\ell_{0}}$-norm optimization:

[TABLE]

where ${\ell_{0}}$-norm, $||{\bf{x}}|{|_{0}}$, denote the number of non-zero entries in the argument and $||{\bf{x}}|{|_{2}}$ indicates its Euclidean norm. Some iterative greedy algorithms have been developed to solve such combinatorial optimization problem, such as Iterative hard thresholding (IHT) [17], Subspace Pursuit (SP) [20], and a series of algorithms derived from them.

III Sparse representation of the radiation flux

In the cylindrical cavity model of ICF, as shown in Figure 2, the radiation region located in the cavity consists of three parts, i.e., spherical surface of the capsule, the bottom and the top surface and the side surface of the cylindrical cavity. Here, $B_{s}$, $B_{u}$ and $B_{w}$ are used to represent the radiation flux of theses three parts respectively and the total radiation flux can be expressed as $B=B_{S}+B_{u}+B_{w}$. Since the X-ray converted from incident laser are uniformly and continuously distributed on the radiation regions, the radiation flux could be represented by linear combinations of some complete orthogonal polynomials. In this section, three orthogonal bases which are constructed by three sets of orthogonal polynomials are used to sparsely represent the radiation flux.

III-A Sparse representation of the radiation flux on capsule

Let $\left(\theta_{s}^{(i)},\varphi_{s}^{(i)}\right)$ be a spherical coordinate of the $i$-th element, then the radiation $B_{s}\left(\theta_{s}^{(i)},\varphi_{s}^{(i)}\right)$ can be represented by a linear combinations of the spherical harmonics as follows [15]:

[TABLE]

where the subscript $n_{s}$ is a polynomial ordering index starting with ${n_{s}}=0$, ${c_{{n_{s}}}}$ is its coefficient and $Y({\theta}_{s},{\varphi}_{s})$ is the spherical harmonics which can be represented as

[TABLE]

where $N_{m_{S}}^{k_{S}}=\sqrt{\frac{\left(\left(2m_{s}+1\right)\left(m_{s}-k_{s}\right)!\right)}{\left(4\pi\left(m_{s}+k_{S}\right)!\right)}}$ is the normalization constant and $P_{m_{s}}^{k_{s}}\left(\cos\theta_{s}\right)$ is the Associated Legendre polynomial. As shown in Figure 2, let $\Delta{\theta}_{s}$ and $\Delta{\varphi}_{s}$ be the minimal angles of subdividing the capsule over the domain $\Omega({\theta}_{s},{\varphi}_{s})=[0,\pi]\times[0,2\pi)$. Equation (10) can be transformed into an approximate form as

[TABLE]

where ${L_{S}}$ is the number of spherical harmonics expansion terms, $\bm{c}^{s}=\left[c_{1},c_{2},\cdots,c_{L_{S}}\right]^{\mathrm{T}}$ is coefficient vector and $\bm{Y}$ is spherical harmonics basis which can be expressed as

[TABLE]

where $N_{{\theta}_{s}}=\frac{\pi}{\Delta\theta_{s}}$ and $N_{\varphi_{s}}=\frac{2\pi}{\Delta\varphi_{s}}$ denote the number of elements divided along the longitude and latitude respectively, and ${N_{S}}$ is the total number of discrete elements on the capsule which can be calculated by $N_{S}=N_{\theta}\times N_{\varphi}$. The 3D pictures of the first 16 terms of the spherical harmonics are shown in Figure 3.

III-B Annular Zernike Polynomials

The bottom and the top surface of the cylindrical cavity is annular because it has two laser entrance holes. As described in [22], Zernike annular polynomials, defined on an annular region, is obtained from the Zernike circle polynomials [22] by using the Gram-Schmidt orthogonalization process. The Zernike annular polynomials are widely used to mathematically describe wavefront aberrations of optical systems due to it is orthogonal over an annular domain, and it can uniquely describe any continuous function with the same definition domain [21]. Let $\left(r_{u}^{(i)},\varphi_{u}^{(i)}\right)$ be a polar coordinate of the $i$-th element of the bottom surface, the annular Zernike polynomials have the following forms in a polar coordinate:

[TABLE]

where $r(0\leq r\leq 1)$ and $\tilde{r}(\tilde{r}\leq r\leq 1)$ are the outer and inner semi-diameter of the annular region respectively, $\varphi_{u}\left(0\leq\varphi_{u}\leq 2\pi\right)$ is the circumferential angle, $l_{u}\geq 0$, $k_{u}\geq 0$ are integers, $l_{u}-\left|k_{u}\right|\geq 0$ and even, $R_{l_{u}}^{k_{u}}(r,\tilde{r})$ is the radial components of the polynomials. The orthonormality of the annular Zernike polynomials can be expressed as

[TABLE]

where $\delta_{l_{u}l^{\prime}_{u}}\left(\delta_{k_{u}k_{u}^{\prime}}\right)$ stands for the Kronecker delta, equal to one if $l_{Z}=l_{Z}^{\prime}\left(k_{Z}=k_{Z}^{\prime}\right)$ and to zero otherwise. when $\tilde{r}=0$, the annular polynomials reduce to the Zernike circle polynomials. Let $l_{u}=2j+k_{u}$, the radial polynomials expression of the annular Zernike polynomials can be obtained by the recurrence formula: [22] when $k_{u}=0$

[TABLE]

where $P_{j}(\cdot)$ is the Legendre polynomials. When $k_{u}>0$, the radial polynomials are given by the recurrence relationship

[TABLE]

where $Q_{j}^{k_{u}}(\tau)$ is a set of orthogonal polynomials obtained by orthogonalizing the sequence $1,\tau,\cdots,\tau^{j}$ over the interval $\left({\tilde{r}}^{2},1\right)$ with a weight function $\tau^{k_{u}}$, $H_{j}^{k_{u}}$ are the normalization constants. the recurrence formula of the $Q_{j}^{k_{u}}(\tau)$ and $H_{j}^{k_{u}}$ as follows:

[TABLE]

[TABLE]

especially, when $k_{u}=0$,

[TABLE]

When $k_{u}<0$, the radial polynomials can be obtained by the relationship:

[TABLE]

When $k_{u}=l_{u}$

[TABLE]

In practice, the outer diameter of the bottom surface is not equal to 1, so we require to map the reality interval $\left[r_{\text{outer}},r_{\text{inner}}\right]$ to the interval $[\tilde{r},1]$. Similar to the radiation flux representation on the capsule, the radiation flux $B_{u}\left(r_{u}^{(i)},\varphi_{u}^{(i)}\right)$ can be expressed by a linear combinations of the annular Zernike polynomials as

[TABLE]

where $n_{u}$ is a polynomial ordering index starting with $n_{u}=0$, $\boldsymbol{U}$ is the annular Zernike polynomials with matrix form which presented as

[TABLE]

where the subscript $N_{r}$ and $N_{\varphi_{u}}$ are the number of elements divided along the radial direction and rotation direction respectively. $N_{u}$ is the total number of elements of one of the end face of the cylinder cavity which can be calculated by $N_{u}=N_{r}\times N_{\varphi_{u}}$. $L_{u}$ is the number of the Zernike annular polynomials expansion terms which can be obtained by the formula $L_{u}=\left(n_{u}+2\right)\left(n_{u}+1\right)/2$. The 3D pictures of the first 15 terms of the annular Zernike polynomials are shown in Figure 4.

III-C Legendre Fourier polynomials

Since the side of the cylindrical cavity is a cylindrical surface, it is no doubt that cylindrical coordinates $\left(\rho,\varphi_{w},z\right)$ seem to be the optimal choice for describing $B_{w}$, where $\rho$, $\varphi_{w}$ and $Z$ denote radius, azimuth angle and height of the cylindrical coordinate. Obviously, $B_{w}$ is only related to azimuth angle and height and not to radius.

According to the construction principle of multivariable orthogonal polynomials, two sets of one-dimensional polynomials which are defined in the azimuth direction and height direction respectively, can be employed to construct a two-dimensional orthogonal polynomials for describing $B_{w}$. Since the cylinder surface profile along the azimuthal direction is closed and periodic, a Fourier series is suitable for characterizing the azimuthal coordinate, which can be expressed as

[TABLE]

where $k_{w}=0,\pm 1,\pm 2,\cdots$.

In height direction, coordinate $z$ is restricted in the interval $[-h,h]$, where $h$ is half the height of the cylindrical cavity. A nature choice for $z$ coordinate is a Fourier series [24]. Similar to the Associated Legendre polynomials used for spherical harmonics, some one-dimensional-polynomials such as Legendre polynomials and Chebyshev polynomials are suitable for the $z$ coordinate. Since Legendre polynomials are orthogonal across the interval $[-1,1]$, which are conveniently mapped to the interval $[-h,h]$, here we consider use Legendre polynomials for characterizing the $z$ coordinate. Then a two-dimensional polynomials can be defined as the tensor products of Fourier series and Legendre polynomials in the variables $z$ and $\varphi_{w}$:

[TABLE]

where $P_{l_{w}}(z)$ is Legendre polynomials with $l_{w}$ degree, $W_{n_{w}}\left(z,\varphi_{w}\right)$ termed Legendre Fourier (LF) polynomials defined in the domain $\Omega=\left\{\left(z,\varphi_{w}\right)\in[0,2\pi]\times[-1,1]\right\}$, $n_{w}$ denotes the ordering index of the LF polynomials starting with $n_{w}=0$. An integer $n_{w}$ represents the degree of the LF polynomials which can be calculated as $n_{w}=l_{w}+\left|k_{w}\right|$, where the operator $|\cdot|$ means taking the absolute value. It is clearly that the LF polynomials are separable in the cylindrical coordinates $z$ and $\varphi_{w}$.

$B_{w}\left(z,\varphi_{w}\right)$ can be approximately presented by the LF polynomials as

[TABLE]

where $c_{n_{w}}$ is the expansion coefficient of the LF polynomials.$\boldsymbol{W}$ denotes the LF orthogonal basis which can be written as

[TABLE]

where the subscript $N_{z}$ and $N_{\varphi_{w}}$ are the number of elements divided along the radial direction and rotation direction respectively. $N_{w}$ is the total number of elements of the side surface of the cylindrical cavity which can be calculated by $N_{w}=N_{z}\times N_{\varphi_{w}}$. $L_{w}$ is the number of the LT polynomials expansion terms. The top 9 expansion terms of the LT polynomials are given in Table I. The 3D pictures of the first 16 terms of the LF polynomials are shown in Figure 5.

III-D Sparse Representation Model of Radiation Flux

From Subsection 3.1-3.3, the total radiation flux in the cavity can be expressed in the form of sum:

[TABLE]

or matrix form:

[TABLE]

Equation (29) can be written as a functional form

[TABLE]

which has a unique unknown variable, i.e., sparse coefficient $\bm{c}$.

In this section, three sets of orthogonal polynomials are used to represent the radiation flux in a cylindrical cavity, and the radiation energy balance equation is transformed into an equation about the sparse coefficients.

IV Radiation flux reconstruction under Compressive Sensing framework

IV-A Compressed Observation Model of Radiation Flux

Compressed Sensing is usually works on linear system. First-order Taylor expansion around a point approximation is usually applied to linearization in nonlinear CS that is similar to the case in Equation (30). Equation (30) can be expanded at point $c^{*}$ as

[TABLE]

where $f_{c^{\star}}$ is a linear operator such as Jacobian matrix which can be expressed as

[TABLE]

where $\operatorname{diag}(\bullet)$ denotes a diagonal matrix whose elements on the diagonal are composed of $\bullet$. Then Equation (32) can be rewritten as:

[TABLE]

As the experimental results shown in Table 2, in the case of smaller grid, the computation is very large for solving view-factor matrix via Newton-Raphson method. In accordance to the CS theory, a measurement matrix of $M\times N(M\ll N)$ can be used for obtaining linear and non-adaptive measurements from Equation (33), which can be expressed as

[TABLE]

The commonly used measurement matrices have stochastic Gaussian matrices, random Bernoulli matrices, partial orthogonal matrices, structured random matrices and cyclic matrices and sparse random matrices [10]. Since the basis $\boldsymbol{\Psi}$ is orthogonal, from which the randomly sampled elements must meet the orthogonal, and also to meet the Restricted Isometry Property (RIP), the sampling points above the elements can be directly selected. This means only needs to calculation the sampling elements composed of measurement matrix. In practical, the measurement matrix is a similar diagonal matrix consisting of the sampling of row index of $f_{c^{\star}}$.

Let $A=\boldsymbol{\Phi}f_{c^{\star}}$, $\boldsymbol{y}=\boldsymbol{\Phi}\left(f_{c^{\star}}c^{\star}-f\left(\boldsymbol{c}^{\star}\right)\right)$, Equation (34) can be written as

[TABLE]

Since $M<N$, in general, recovery of $\bm{c}$ is a combination problem which is known to be NP-hard. Fortunately, under stricter conditions on $\bm{A}$ (i.e., Restricted Isometry Property), a range of different algorithms can be used to recover $\bm{c}$ from roughly $O(K\log(N))$ measurements [9]. Based on the results of Section 3, we can roughly estimate the number of samples for capsule, the bottom and the top surface of the cavity and the side of the cavity, respectively. For example, for a model with a total elements number of 9776 (where the the number of elements for capsule, the bottom and top surface of the cavity, and the side of the cavity are 2592, 2016, 2016 and 3152, respectively), the number of samples for capsule is $30\times\log(2592)=102.4$, the bottom and top surface have the same samples of $35\times\log(2016)=115.6$, and the number of samples on the side of cavity is $100\times\log(3152)=349.9$, and the total number of samples can be estimated as $150+150\times 2+400=850$. The number of samples in other models is shown in Table III.

IV-B Algorithms for solving model

Sparse coefficient $\bm{c}$ can be recovered by optimizing the following sparse constraint problems

[TABLE]

some greedy iterative algorithms such as hard thresholding algorithms and orthogonal matching algorithms can be used for solving Equation (36). In this section, we proposed a new iterative greedy algorithm for solving the radiation flux compressed observation model based on two typical CS reconstruction algorithms.

IV-B1 Conjugate Gradient Iterative Hard Thresholding Algorithm

Iterative hard thresholding (IHT) algorithm is usually used for CS reconstruction due to its simple operation and low computational complexity, and it’s iterative formula as

[TABLE]

where $H_{K}(\bullet)$ is a nonlinear operator that sets all but the absolute value of largest $K$ elements of $\bm{c}$ to zero, $\mu$ is the step size and usually taken as a constant in IHT algorithm. The convergence of this algorithm was proven in [dd] under the condition that $\|\bm{A}\|_{2}<1$, in which case, the IHT algorithm converges to a local minimum. Since step size and the scale of sensing matrix have a great influence on the convergence of IHT [17], Normalized Iterative Hard Thresholding algorithm (NIHT) was proposed in to improve IHT algorithm [18]. The restriction of the snesing matrix $\bm{A}$ in NIHT is replaced by a descending iteration factor, which greatly improves the convergence and stability. However, the search direction in NIHT algorithm is gradient direction, it may result slow convergence near the minimum. Conjugate Gradient Iterative Hard Thresholding (CGIHT) algorithm was proposed in [19], in which, conjugate gradient direction was adopted as the search direction to speed up the convergence. The CGIHT algorithm can be stated as Algorithm 1, in which, $\emptyset$ denotes an empty set and $\bm{d}$ denotes search direction, $\gamma$ is an orthogonalization weight, the support of $\bm{c}$ denoted by $T_{n}$ and $\bm{A}_{T_{n}}$ indicates a submatrix of $\bm{A}$. By allowing the search direction change after each low complexity iteration, CGIHT is able to quickly identify the correct descending direction while using the computational advantages of conjugate gradient method when the support set remains stable.

Although the CGIHT algorithm can accelerate convergence, the Jacobian matrix needs updating at each iteration in non-linear problems, which may result in search directions that do not satisfy the conjugate condition, resulting in slower convergence. In the process of solving nonlinear radiation flux calculation model, it is time-consuming to update the Jacobian matrix, and we should to devise means to reduce the number of iterations.

IV-B2 Subspace Pursuit Algorithm

Orthogonal matching pursuit (OMP) algorithm is wildly used for CS reconstruction owing to its simple implementation and lower computational complexity [20]. In each iteration, OMP selects the column of the measurement matrix which is most strongly correlated with the current residuals and adds this column into the set of selected columns, then the residual is updated by projecting the measurement vector $\bm{y}$ onto the linear subspace spanned by the selected columns. The reconstruction complexity of the standard OMP is roughly $O(KMN)$ because it always needs $s$ iterations. Unlike the OMP algorithm, at each iteration, subspace pursuit algorithm selected $s$ candidates based on the correlation values between the columns of $\bm{A}$ and the measurement vector $\bm{y}$, and then a back-tracking operation of the candidate set in the previous iteration is implemented. Compared to the OMP algorithm, the computational complexity of the SP algorithm can be further reduced to $O(MN\log(K))$ when the nonzero entries of the sparse signal decays slowly. The SP algorithm can be described as follows.

IV-B3 Conjugate Gradient Subspace Thresholding Pursuit Algorithm

Inspired by CGIHT algorithm and SP algorithm, a new greedy iterative algorithm named conjugate gradient subspace thresholding pursuit (CGSTP) is proposed for reconstructing radiation flux. This algorithm is based on SP algorithm, drawing on CGIHT algorithm and using conjugate gradient as the search direction to accelerate convergence. The pseudo-code of CGSTP algorithm is shown in Algorithm 3.

The CGHTP algorithm is a simple combination of the CGIHT algorithm and the SP algorithm, and absorbs the advantages of both algorithms. SP algorithm introduces backtracking operation in the iteration process. In each iteration, the support set of the previous iteration is reconsidered, and the new support set is updated by the least square method. For compressible signals, SP algorithm can find the best $s$ atoms in $log(K)$ iteration. However, SP algorithm has strict requirements for RIP conditions.

The CGIHT algorithm optimizes the search direction, and the conjugate direction is used as one of its optimization directions, which has a superlinear convergence rate associated with $(\sqrt{\kappa}-1)/(\sqrt{\kappa}+1)$ [25], where $\kappa$ is the condition number of matrix $\bm{A}$. However, it does not make further accurate estimation of the sparse coefficient after updating the support set, and the CGIHT algorithm only shows its advantage in convergence speed when the correlation of the columns in $A$ is very small.

The CGIHT algorithm uses the negative gradient direction and the conjugate gradient direction as the search direction. If the sensing matrix remains unchanged in all steps, only a small number of iterations are needed to converge. In fact, since the incomplete orthogonality between the columns of the sensing matrix, the support set selected in each iteration may be different, and the sub-matrix formed by the columns of $\bm{A}$ with the support set as the index is different in each iteration, which will lead to the search directions obtained in the previous iteration and the current iteration are not conjugate, thus the superlinear convergence speed of the algorithm can not be ensured. What’s more, in the nonlinear problem, because the CGIHT algorithm does not accurately estimate the sparse coefficients in each iteration, the solution process requires more iterations.

The CGSTP algorithm absorbs the advantages of SP algorithm and CGIHT algorithm in support set update strategy and search direction optimization respectively. CGSTP provides two alternative search directions in the iterative process, including the gradient descent direction and the conjugate gradient direction, and adaptively selects the optimal direction according to the criterion of whether the support set is the same in the previous iteration and the current iteration. In theory, if the sensing matrix in all iterations remains unchanged, CGSTP can achieve superlinear convergence rate [25]. However, in solving the non-linear compressed observation model of radiation flux, the Jacobian matrix changes in each iteration, which results in the change of the sensing matrix, and the acceleration effect is not obvious. Fortunately, since CGSTP uses a backtracking operation similar to SP algorithm to accurately estimate the sparse coefficients after each updating of the support set, resulting in a small change in the Jacobian matrix obtained from the previous iteration and the current iteration, which increases the probability that the sensing matrices in different iterations are the same.

IV-C Overview of radiation flux computation process with CGSTP

According to the above analysis, solving the non-linear radiation energy balance equation in the Compressed Sensing framework can be summarized into three modules, namely, the sparse representation module of radiation flux, the compressive sampling module and the reconstruction module. Figure 6 summarizes the basic contents of each module.

As can be seen from Figure 6, in the sparse representation module, all the radiation regions of the cylindrical cavity model are first divided into discrete elements according to a given precision, and parameter data such as coordinate values, normal vectors, and spatial angles are obtained. Then, the sparse basis composed of orthogonal polynomials of different expansion orders is calculated according to the elements parameters. The representation error of radiation flux over the sparse basis is analyzed to determine the number of expansion terms of polynomials, and the sparsity of radiation flux in the sparse transformation domain is determined by simulation experiments. Finally, a sparse representation model of radiation energy for a small number of sparse coefficients is constructed.

In the compressive sampling module, firstly, all discrete mesh elements are indexed, and $M$ elements are collected by random uniform sampling method (Latin Hypercube Sampling (LHS) method used in this paper). Then, the sparse representation model of radiation flux is transformed into a compressed observation model by using the equation corresponding to $M$ sampling indexes, so as to realize the compression and dimensionality reduction of the model. At the same time, the occlusion condition of each element in the cylindrical cavity is analyzed, and the partial view-factor matrix is calculated according to the sampling index.

In the reconstruction module, the linear approximation system is first constructed by using the first-order Taylor expansion at the current iteration point to construct a linear approximation system. Then, the partial Jacobian matrix is calculated according to the sampling index, and the sensing matrix is calculated according to the Jacobian matrix, The minimization problem based on ${\ell_{0}}$ norm is established when the sparsity is known as $s$. Finally, the above ${\ell_{0}}$ norm minimization problem is solved by CGSTP algorithm.

V Simulation Experiments and Discussion

V-A Experimental model and simulation environment

In this section, two cylindrical cavity models for Shenguang II and Shenguang-III laser devices are used for simulation experiments to verify the performance of the proposed method. The cylindrical cavity model in Shenguang II is shown in Figure 7. The models of Shenguang-II and Shenguang-III are represented by ”S2” and ”S3” respectively. All experiments are tested in MATLAB R2015b [26], running on a Windows 10 machine with Intel I5-7500 CPU 3.4GHz and 8GB RAM.

The elements parameters of the S2 model and the S3 model are shown in Table II, where S2-1, S2-2, S3-1, and S3-2 represent models having different elements sizes and their total elements numbers are respectively 9776, 38952, 20736 and 82944.

The random uniform sampling method used in this section is the Latin Hypercube Sampling (LHS) method. The sampling number is estimated according to the formula $s\log(N)$. Taking the S21 model as an example, the sampling numbers of three radiation regions are calculated as $30\times\log(2592)\approx 102.4$ (take 150), $35\times\log(2016)=115.6$ (take 150) and $100\times\log(3152)=349.8$ (take 400). The sampling rates of different models are shown in Table III.

As shown in the Table II and Table III, we can see that (1) the number of discrete mesh elements significantly increases with the decrease of the size of mesh elements, and (2) the sampling rate for Compressed Sensing based approaches gradually decreases as the number of mesh elements increases.

V-B Sparse representation error analysis of radiation flux

Firstly, we verify the sparsity level of radiation flux ${\bm{B}}_{s}$, ${\bm{B}}_{u}$ and ${\bm{B}}_{w}$ over three orthogonal bases $\bm{Y}$, $\bm{U}$ and $\bm{W}$, respectively. The expansion coefficients of three sets of polynomials and the representation errors (root mean squared error), between the original flux and the recovered radiation flux with $n_{y}$, $n_{u}$ and $n_{w}$ expansion order are shown in Figure 8(a) and (b) respectively. The relationship between the sparsity level of the coefficients and the representation error of the radiant energy flow is shown in Figure 9(a) and the magnitude of the sparse coefficients are presented in 9(b).

As shown in Figure 8(a) and (b), with the increase of the number of expansion terms of orthogonal polynomials, the representation error decreases gradually. When the expansion terms of Spherical Harmonic polynomials, Annular Zernike polynomials and Legendre-Fourier polynomials reaches 400, 325 and 1225 respectively, the representation errors approach to zero and the curve remains stable. This means that only a small number of expansion terms are needed to accurately represent the complete radiation flux. From Figure 9 (a) we can see that the number of the absolute value of coefficients larger than $10^{-3}$ over three sets of orthogonal bases are $s1\approx 30$, $s2\approx 35$ and $s3\approx 100$ respectively.

Since the top 35 order energy accounts for more than 99.9% of the total energy and the energy of more than 35th order is little. That is the analysis of radiation flux’s harmonic expansion situation under various irradiation conditions. The drive asymmetry must be less than 2% and the calculation accuracy should be setting at least an order for magnitude, equal to 0.1%.

V-C Solution results of radiation energy balance equations

The resulting radiation flux on the capsule is computed with TDEBM and then compared with Newton-Raphson (NR), Preconditioned Conjugate Gradient (PCG), IHT, CGIHT and CGSTP on computation time, accuracy, iteration number for different kinds of mesh elements sizes. In all tests, the IHT, cgiht, and cgstp algorithms were run independently 20 times, and the average of 20 results was used as the test result. Then the computation time list for the four calculation models using different algorithms is shown in Table IV and V, and the iteration steps and accuracy of different algorithms are shown in Table VI and Table VII respectively.

As shown in Table IV and Table V, we can see that (1) as the number of non-linear equations exponentially increases, the computation times thus increase sharply for traditional approaches, (2) less computation time is required for compressed sensing based approaches since the computation model is significantly reduced, only almost 8.7%, or 1.1% of equations are utilized to obtain the solution, and (3) since the search direction is optimized and the updating strategy of support set is improved, less time required for CGSTP than that of IHT and CGIHT to obtain the solution. Figure 10(a) shows the time ratio between the CGSTP algorithm and the PCG algorithm for calculating the four models. Figure 11(b) shows the iteration time difference between the CGSTP algorithm and the IHT and CGIHT algorithms. It can be seen from Figure 10(a) that as the model scale increases, the time acceleration ratio of the CGSTP algorithm relative to the PCG method becomes larger and larger, and the maximum acceleration ratio can reach 80. This is because as the scale of the equation increases, the sampling rate becomes smaller and smaller, and the required calculation amount is less and less than the total calculation amount. Figure 10 (b) shows that the difference between CGSTP and IHT and CGIHT is larger and larger with the increase of model size. This is because the radiation flux tends to be continuous with the increase of the number of discrete mesh elements, which reduces the number of updates of Jacobian matrix. This means that CGSTP has more advantages when the scale is larger.

The convergence performance of the above algorithms is shown in Table VI.It can be seen that (1) a small number of iterations is required for the NR, PCG and CGSTP approaches, whereas hundreds of iterations is required for the IHT, CGIHT approaches, (2) Even though less time of iterations are required for the NR and PCG approaches, the computation time for each iteration is much longer than that of IHT, CGIHT and CGSTP algorithms, and (3) as compared with the IHT and CGIHT algorithms, less iterations is required for the CGSTP approach since optimal search direction and step length is applied. Figure 11 shows the relationship between IHT, CGIHT and CGSTP algorithm iteration times and reconstruction error. From Figure 11, we can see that the error curve of CGSTP algorithm decreases very fast, and it only needs about six iterations to approach zero and keep stable, which means that the CGSTP algorithms has a faster convergence speed than the other algorithms.

The accuracy of the different algorithms is listed in Table VII. It can be seen from Table VII that (1) the NR algorithm and the PCG algorithm have the highest accuracy, and the root mean square error (RMSE) is lower than $10^{-5}$, (2) the accuracy of CGSTP algorithm is higher than IHT algorithm and CGIHT algorithm in solving different models, and (3) the accuracy of all algorithms is improved with the increase of the number of mesh elements. The RMSE of CGSTP algorithm in calculating different models is less than $10^{-3}$, which satisfies the radiation symmetry analysis.

VI Conclusions

The radiation flux distribution symmetry evaluation is very important in ICF experiments, which involve a time consumption process in solving the nonlinear energy equilibrium computation model. In order to accelerate such equation solving process, an efficient radiation computation approach is presented. Firstly, we investigated the distribution characteristics of radiation flux on a cylindrical cavity model and three sets of orthogonal polynomials are employed to accurately and sparsely represent the radiation flux. Then, only a few mesh elements are sampled to formulate the sparse equation model, which enable the non-linear equation model largely reduced. Finally, a greedy iterative algorithm named CGSTP is presented to efficiently solve the non-linear radiation energy balance equations. Some experimental targets is utilized to validate the efficiency of the presented approach. The result show that the radiation flux can be efficiently obtained in almost ${1\mathord{\left/{\vphantom{1{80}}}\right.\kern-1.2pt}{80}}$ time of traditional approach. In future, we will continue to explore the sparse representation for recent free-form hohlraums, to enable the present approach be used for more targets design, analysis and optimization.