Introduction
Heat conduction departs from Fourier’s law when device dimensions become comparable to phonon mean free paths [1 , 2 , 3 ] . The phonon Boltzmann transport equation (BTE) describes this regime by following the phonon distribution in physical and velocity space. Combined with the Callaway relaxation model [4 , 5 ] , which distinguishes momentum-destroying resistive (R) scattering from momentum-conserving normal (N) scattering, a single equation covers ballistic, diffusive and hydrodynamic transport, including phonon Poiseuille flow [6 , 7 , 8 , 9 , 10 ] . This generality has a price: deterministic solutions of the BTE are expensive, and design studies or uncertainty quantification, which repeat them over ranges of scattering rates, multiply that expense.
The cost has two distinct sources. The first is the size of the state. Even for a gray model in two spatial dimensions, a discrete-ordinate method stores the distribution at every spatial degree of freedom and every ordinate, so storage grows with the product of the spatial and angular resolutions. Monte Carlo methods avoid this storage at the price of statistical noise, which variance-reduced formulations reduce substantially [11 ] . The second source is slow convergence. Source iteration alternates a transport sweep with an update of the lagged scattering source and propagates information by about one mean free path per sweep, so its error-reduction factor approaches one as collisions become frequent [12 ] . Diffusion synthetic acceleration [13 ] , multiscale kinetic schemes [14 , 15 , 16 , 17 , 18 ] and the general synthetic iterative scheme (GSIS) [19 , 20 , 21 ] remove this stiffness. Between kinetic sweeps, GSIS solves macroscopic moment equations whose higher-order terms are evaluated from the current kinetic state. For the phonon Callaway model it reaches the steady state in tens to a few hundred iterations over the whole range of Knudsen numbers, where source iteration may need thousands [21 ] . It does not, however, reduce the size of the state.
Low-rank methods target this first source of cost. Proper generalized decomposition (PGD) [22 , 23 , 24 , 25 ] builds a sum of separated products directly from the governing equation, without precomputed snapshots: a greedy procedure adds one product at a time, and alternating Galerkin updates determine its factors. Model parameters can be treated as additional coordinates, so that one representation covers a whole parameter domain [26 , 27 ] . PGD has been applied to neutron transport [28 ] and to rarefied gas flow [29 ] ; for the parametrized Shakhov model, adding synthetic equations to the PGD iteration makes the reduced solver recover the Navier–Stokes limit [30 ] . Reduced-basis methods [31 ] instead build low-dimensional spaces from full-order solutions, and dynamical low-rank approximation evolves factorized solutions in time [32 , 33 , 34 , 35 ] ; its asymptotic-preserving variant based on a micro–macro decomposition [36 ] shares our aim of obtaining the macroscopic limit independently of the rank.
Combining PGD with synthetic acceleration raises two questions. The first is which splitting to use. In transport problems the upwind direction ties space to angle, and a splitting that separates space from direction altogether loses this structure, because the spatial factor then sees only direction-averaged streaming coefficients. The second is how the separated iteration converges. A separated representation does not change how macroscopic information propagates. If the scattering source is taken from the moments of the current representation, the iteration remains a source iteration in disguise: even when this coupling is treated implicitly in the polar factor, macroscopic information crosses the domain only through transport sweeps of the spatial factor, and enrichment adds modes to an unconverged state. This stall occurs precisely where the solution is most compressible: the leading Fourier and Guyer–Krumhansl equilibria have rank at most two in the splitting used below. The synthetic step of GSIS, in contrast, needs only a few angular moments of the kinetic state, which a separated representation supplies cheaply. Whether the truncation and alternation errors of PGD interfere with that step requires analysis of the PGD iteration itself and of the matching between reconstructed and synthetic fields.
In this paper we combine PGD with GSIS for the steady linearized gray Callaway model on two-dimensional domains. Our contributions are the following.
A splitting with a Knudsen-uniform rank bound. The distribution is separated into a space–azimuth factor, which carries the upwind structure and is computed by standard discontinuous Galerkin (DG) sweeps, and a polar factor. Because the polar angle only stretches the in-plane optical depth, the separation rank is bounded, uniformly in the Knudsen numbers, by the numerical rank of the kernel e − τ / sin ϕ e^{-\tau/\sin\phi} e − τ / s i n ϕ P N P_N P N N + 1 N+1 N + 1
Synthetic coupling to the separated state. PGD–GSIS freezes the macroscopic source within each alternating pair and refreshes it by the synthetic equations, discretized by a hybridizable DG (HDG) method. Their high-order loads, including DG face terms, and the half-space wall closure are assembled from only four polar moments of each mode (Section 4.2 and Appendix D).
A convergence analysis of pure PGD and PGD–GSIS. Pure PGD, although implicit in the polar factor, moves macroscopic information only through lagged spatial sweeps and inherits the stall of source iteration. Each PGD–GSIS step splits into the full-order GSIS error propagation [21 ] and the image of the representation error, which the synthetic map reads only through third-order angular moments. In the ballistic limit the roles reverse (Section 4.3).
An analysis of the asymptotic limits. A formal expansion of the iteration identifies the required polar space; explicit conditions on the polar space, the projected residual and the source increment then preserve the leading Fourier flux and the viscous GK balance, and in the ballistic limit the vanishing collision rate cancels the O ( K n C ) \mathcal{O}(\mathrm{Kn}_{C}) O ( Kn C )
A parametric extension, in which the polar factor also depends on both inverse Knudsen numbers, so that one spatial basis covers four decades of each. A joint update of all polar factors, one small linear system per polar and parameter node, is required for accuracy at large Knudsen numbers. We identify what sets the rank and where the error concentrates (Section 5).
Section 2 introduces the model, its macroscopic structure and the test problems. Sections 3–5 develop and analyze the separated representation, the alternating iteration and the parametric extension. Section 6 reports results for an isothermal cavity, cross-plane conduction and phonon Poiseuille flow, and Section 7 concludes. The appendices contain the asymptotic analysis, the Fourier form of the iteration, and the DG and HDG discretizations.
Kinetic model and its macroscopic structure
Linearized gray Callaway equation
Positions x = ( x , y ) \mathbf{x}=(x,y) x = ( x , y ) D D D
Ω = ( sin ϕ cos θ , sin ϕ sin θ , cos ϕ ) = ( sin ϕ s ( θ ) , cos ϕ ) , d Ω = sin ϕ d ϕ d θ , (1) \boldsymbol{\Omega}=(\sin\phi\cos\theta,\sin\phi\sin\theta,\cos\phi)=(\sin\phi\,\mathbf{s}(\theta),\cos\phi),\qquad\mathrm{d}\Omega=\sin\phi\,\mathrm{d}\phi\,\mathrm{d}\theta, \tag*{(1)} Ω = ( sin ϕ cos θ , sin ϕ sin θ , cos ϕ ) = ( sin ϕ s ( θ ) , cos ϕ ) , d Ω = sin ϕ d ϕ d θ , (1) with azimuth θ ∈ [ 0 , 2 π ) \theta\in[0,2\pi) θ ∈ [ 0 , 2 π ) ϕ ∈ [ 0 , π ] \phi\in[0,\pi] ϕ ∈ [ 0 , π ] T 0 T_0 T 0 I ~ ( x ~ , Ω ) \widetilde{I}(\widetilde{\mathbf{x}},\boldsymbol{\Omega}) I ( x , Ω ) ∫ I ~ d Ω = C v ( T ~ − T 0 ) \int\widetilde{I}\,\mathrm{d}\Omega=C_v(\widetilde{T}-T_0) ∫ I d Ω = C v ( T − T 0 ) q ~ = v g ∫ Ω I ~ d Ω \widetilde{\mathbf{q}}=v_g\int\boldsymbol{\Omega}\widetilde{I}\,\mathrm{d}\Omega q = v g ∫ Ω I d Ω [4 , 5 , 21 ] reads
v g Ω ⋅ ∇ I ~ = I ~ R e q − I ~ τ R + I ~ N e q − I ~ τ N , I ~ R e q = C v ( T ~ − T 0 ) 4 π , I ~ N e q = I ~ R e q + 3 Ω ⋅ q ~ 4 π v g . (2) v_g\boldsymbol{\Omega}\cdot\nabla\widetilde{I}=\frac{\widetilde{I}^{\mathrm{eq}}_R-\widetilde{I}}{\tau_R}+\frac{\widetilde{I}^{\mathrm{eq}}_N-\widetilde{I}}{\tau_N},\qquad\widetilde{I}^{\mathrm{eq}}_R=\frac{C_v(\widetilde{T}-T_0)}{4\pi},\qquad\widetilde{I}^{\mathrm{eq}}_N=\widetilde{I}^{\mathrm{eq}}_R+\frac{3\boldsymbol{\Omega}\cdot\widetilde{\mathbf{q}}}{4\pi v_g}. \tag*{(2)} v g Ω ⋅ ∇ I = τ R I R eq − I + τ N I N eq − I , I R eq = 4 π C v ( T − T 0 ) , I N eq = I R eq + 4 π v g 3 Ω ⋅ q . (2) With a length L L L Δ T \Delta T Δ T x = x ~ / L \mathbf{x}=\widetilde{\mathbf{x}}/L x = x / L I = I ~ / ( C v Δ T ) I=\widetilde{I}/(C_v\Delta T) I = I / ( C v Δ T ) T = ( T ~ − T 0 ) / Δ T T=(\widetilde{T}-T_0)/\Delta T T = ( T − T 0 ) /Δ T q = q ~ / ( C v v g Δ T ) \mathbf{q}=\widetilde{\mathbf{q}}/(C_vv_g\Delta T) q = q / ( C v v g Δ T ) K n R , N = v g τ R , N / L \mathrm{Kn}_{R,N}=v_g\tau_{R,N}/L Kn R , N = v g τ R , N / L δ R = K n R − 1 \delta_R=\mathrm{Kn}_R^{-1} δ R = Kn R − 1 δ N = K n N − 1 \delta_N=\mathrm{Kn}_N^{-1} δ N = Kn N − 1 δ C = δ R + δ N = K n C − 1 \delta_C=\delta_R+\delta_N=\mathrm{Kn}_C^{-1} δ C = δ R + δ N = Kn C − 1
Ω ⋅ ∇ I + δ C I = δ C 4 π T + 3 δ N 4 π Ω ⋅ q , T = ∫ I d Ω , q = ∫ Ω I d Ω . (3) \boldsymbol{\Omega}\cdot\nabla I+\delta_C I=\frac{\delta_C}{4\pi}T+\frac{3\delta_N}{4\pi}\boldsymbol{\Omega}\cdot\mathbf{q},\qquad T=\int I\,\mathrm{d}\Omega,\qquad\mathbf{q}=\int\boldsymbol{\Omega}I\,\mathrm{d}\Omega. \tag*{(3)} Ω ⋅ ∇ I + δ C I = 4 π δ C T + 4 π 3 δ N Ω ⋅ q , T = ∫ I d Ω , q = ∫ Ω I d Ω. (3) The right-hand side is δ C \delta_C δ C ( T + 3 χ Ω ⋅ q ) / ( 4 π ) (T+3\chi\,\boldsymbol{\Omega}\cdot\mathbf{q})/(4\pi) ( T + 3 χ Ω ⋅ q ) / ( 4 π ) χ = δ N / δ C ∈ [ 0 , 1 ) \chi=\delta_N/\delta_C\in[0,1) χ = δ N / δ C ∈ [ 0 , 1 ) ϕ ↦ π − ϕ \phi\mapsto\pi-\phi ϕ ↦ π − ϕ q = ( q x , q y ) \mathbf{q}=(q_x,q_y) q = ( q x , q y )
Moment equations and the synthetic system
With the trace-free moments
N i j = ∫ ( Ω i Ω j − δ i j 3 ) I d Ω , M i j k = ∫ K i j k ( Ω ) I d Ω , (4) N_{ij}=\int\left(\Omega_i\Omega_j-\frac{\delta_{ij}}{3}\right)I\,\mathrm{d}\Omega,\qquad M_{ijk}=\int K_{ijk}(\boldsymbol{\Omega})I\,\mathrm{d}\Omega, \tag*{(4)} N ij = ∫ ( Ω i Ω j − 3 δ ij ) I d Ω , M ij k = ∫ K ij k ( Ω ) I d Ω , (4) where K i j k = Ω i Ω j Ω k − ( Ω i δ j k + Ω j δ i k + Ω k δ i j ) / 5 K_{ijk}=\Omega_i\Omega_j\Omega_k-(\Omega_i\delta_{jk}+\Omega_j\delta_{ik}+\Omega_k\delta_{ij})/5 K ij k = Ω i Ω j Ω k − ( Ω i δ j k + Ω j δ ik + Ω k δ ij ) /5 δ i j \delta_{ij} δ ij 1 1 1 Ω i \Omega_i Ω i Ω i Ω j − δ i j / 3 \Omega_i\Omega_j-\delta_{ij}/3 Ω i Ω j − δ ij /3 (3) are, exactly (repeated indices summed, ∂ z = 0 \partial_z=0 ∂ z = 0
∂ k q k = 0 , 1 3 ∂ i T + ∂ k N i k = − δ R q i , 2 5 ∂ ⟨ j q i ⟩ + ∂ k M i j k = − δ C N i j , (5) \partial_kq_k=0,\qquad\frac{1}{3}\partial_iT+\partial_kN_{ik}=-\delta_Rq_i,\qquad\frac{2}{5}\partial_{\langle j}q_{i\rangle}+\partial_kM_{ijk}=-\delta_CN_{ij}, \tag*{(5)} ∂ k q k = 0 , 3 1 ∂ i T + ∂ k N ik = − δ R q i , 5 2 ∂ ⟨ j q i ⟩ + ∂ k M ij k = − δ C N ij , (5) with ∂ ⟨ j q i ⟩ = ( ∂ j q i + ∂ i q j ) / 2 − δ i j ∂ k q k / 3 \partial_{\langle j}q_{i\rangle}=(\partial_jq_i+\partial_iq_j)/2-\delta_{ij}\partial_kq_k/3 ∂ ⟨ j q i ⟩ = ( ∂ j q i + ∂ i q j ) /2 − δ ij ∂ k q k /3 N i j N_{ij} N ij H i j = ∂ k M i j k H_{ij}=\partial_kM_{ijk} H ij = ∂ k M ij k
1 3 ∇ T + δ R q − K n C 5 [ ∇ 2 q + 1 3 ∇ ( ∇ ⋅ q ) ] = K n C ∇ ⋅ H , ∇ ⋅ q = 0. (6) \frac{1}{3}\nabla T+\delta_R\mathbf{q}-\frac{\mathrm{Kn}_C}{5}\left[\nabla^2\mathbf{q}+\frac{1}{3}\nabla(\nabla\cdot\mathbf{q})\right]=\mathrm{Kn}_C\nabla\cdot\mathbf{H},\qquad\nabla\cdot\mathbf{q}=0. \tag*{(6)} 3 1 ∇ T + δ R q − 5 Kn C [ ∇ 2 q + 3 1 ∇ ( ∇ ⋅ q ) ] = Kn C ∇ ⋅ H , ∇ ⋅ q = 0. (6) No closure has been made: if H H H (6) is an exact rearrangement of Eq. (5) .
This is the macroscopic (synthetic) system of GSIS [20 , 21 ] . It has the form of the Brinkman equations, with q \boldsymbol{q} q T / 3 T/3 T /3 δ R \delta_R δ R K n C / 5 \mathrm{Kn}_C/5 Kn C /5 H H H K n C \mathrm{Kn}_C Kn C H = 0 H=0 H = 0 (6) by K n R \mathrm{Kn}_R Kn R [6 , 7 ]
q − ℓ G K 2 [ ∇ 2 q + 1 3 ∇ ( ∇ ⋅ q ) ] = − K n R 3 ∇ T , ℓ G K 2 = K n R K n C 5 , (7) \boldsymbol{q}-\ell_{\mathrm{GK}}^2\left[\nabla^2\boldsymbol{q}+\frac{1}{3}\nabla(\nabla\cdot\boldsymbol{q})\right]=-\frac{\mathrm{Kn}_R}{3}\nabla T,\qquad\ell_{\mathrm{GK}}^2=\frac{\mathrm{Kn}_R\mathrm{Kn}_C}{5}, \tag*{(7)} q − ℓ GK 2 [ ∇ 2 q + 3 1 ∇ ( ∇ ⋅ q ) ] = − 3 Kn R ∇ T , ℓ GK 2 = 5 Kn R Kn C , (7) for smooth bulk fields with K n C ≪ L g \mathrm{Kn}_C\ll L_g Kn C ≪ L g L g L_g L g § 8 ); at K n C ≳ L g \mathrm{Kn}_C\gtrsim L_g Kn C ≳ L g ℓ G K / L g \ell_{\mathrm{GK}}/L_g ℓ GK / L g K n R \mathrm{Kn}_R Kn R K n N → 0 \mathrm{Kn}_N\to0 Kn N → 0 ℓ G K → 0 \ell_{\mathrm{GK}}\to0 ℓ GK → 0
Boundary conditions
Every boundary prescribes the incoming half-range Ω ⋅ n < 0 \boldsymbol{\Omega}\cdot\boldsymbol{n}<0 Ω ⋅ n < 0 n \boldsymbol{n} n thermalizing wall at temperature T b T_b T b I = T b / ( 4 π ) I=T_b/(4\pi) I = T b / ( 4 π ) diffuse adiabatic wall emits an isotropic intensity whose level cancels the outgoing normal energy flux,
I b = ∫ Ω ⋅ n > 0 ( Ω ⋅ n ) I d Ω − ∫ Ω ⋅ n < 0 ( Ω ⋅ n ) d Ω ; (8) I_b=\frac{\int_{\boldsymbol{\Omega}\cdot\boldsymbol{n}>0}(\boldsymbol{\Omega}\cdot\boldsymbol{n})I\,\mathrm{d}\Omega}{-\int_{\boldsymbol{\Omega}\cdot\boldsymbol{n}<0}(\boldsymbol{\Omega}\cdot\boldsymbol{n})\,\mathrm{d}\Omega}; \tag*{(8)} I b = − ∫ Ω ⋅ n < 0 ( Ω ⋅ n ) d Ω ∫ Ω ⋅ n > 0 ( Ω ⋅ n ) I d Ω ; (8) numerator and denominator are evaluated with the same discrete ordinates, so the discrete wall flux vanishes exactly. On periodic faces paired by a translation Δ \boldsymbol{\Delta} Δ Δ T b / ( 4 π ) \Delta T_b/(4\pi) Δ T b / ( 4 π ) Δ T b ≠ 0 \Delta T_b\ne0 Δ T b = 0 Δ \boldsymbol{\Delta} Δ § 3.2 . Without a thermalizing wall the solution is defined up to I ↦ I + c / ( 4 π ) I\mapsto I+c/(4\pi) I ↦ I + c / ( 4 π )
Test problems
The three problems of Figure 1 are the gray-Callaway benchmarks used to validate full-order GSIS [21 ] ; each isolates a different macroscopic mechanism.
Figure 1. Test problems (not to scale). Solid edges are thermalizing (red hot, blue cold) or diffusely adiabatic (olive); dashed edges are periodic. The channels have height H = 0.2 H = 0.2 H = 0.2
Isothermal cavity. The unit square has its top wall at T b = 1 T_b=1 T b = 1 T b = 0 T_b=0 T b = 0
Cross-plane conduction. The channel [ 0 , 1 ] × [ 0 , H ] [0,1]\times[0,H] [ 0 , 1 ] × [ 0 , H ] H = 0.2 H=0.2 H = 0.2 x = 0 x=0 x = 0 T b = 0 T_b=0 T b = 0 x = 1 x=1 x = 1 T b = 1 T_b=1 T b = 1 y = 0 , H y=0,H y = 0 , H x x x ∇ ⋅ q = 0 \nabla\cdot\boldsymbol{q}=0 ∇ ⋅ q = 0 q x q_x q x (7) vanishes identically. The bulk gradient is then ∂ x T = − 3 δ R q x \partial_x T = -3\delta_R q_x ∂ x T = − 3 δ R q x δ R \delta_R δ R
Phonon Poiseuille flow. The same channel has diffusely adiabatic walls at y = 0 y = 0 y = 0 H H H x = 0 , 1 x = 0, 1 x = 0 , 1 Δ T b = 1 \Delta T_b = 1 Δ T b = 1 x x x (7) gives q x ( y ) ∝ 1 − cosh [ ( y − H / 2 ) / ℓ G K ] / cosh [ H / ( 2 ℓ G K ) ] q_x(y) \propto1 - \cosh[(y - H/2)/\ell_{\rm GK}]/\cosh[H/(2\ell_{\rm GK})] q x ( y ) ∝ 1 − cosh [( y − H /2 ) / ℓ GK ] / cosh [ H / ( 2 ℓ GK )] ℓ G K ≪ H \ell_{\rm GK} \ll H ℓ GK ≪ H ℓ G K ≫ H \ell_{\rm GK} \gg H ℓ GK ≫ H ≈ 5 H 2 / ( 24 K n C ) \approx5H^2/(24\mathrm{Kn}_{\mathrm{C}}) ≈ 5 H 2 / ( 24 Kn C ) K n R \mathrm{Kn}_{\mathrm{R}} Kn R
Separated representation and its rank
Space–azimuth by polar splitting
For fixed ( K n R , K n N ) (\mathrm{Kn}_{\mathrm{R}}, \mathrm{Kn}_{\mathrm{N}}) ( Kn R , Kn N )
I m ( x , θ , ϕ ) = ∑ r = 1 m X r ( x , θ ) Φ r ( ϕ ) . (9) I_m(x,\theta,\phi) = \sum_{r=1}^{m} X_r(x,\theta)\,\Phi_r(\phi). \tag*{(9)} I m ( x , θ , ϕ ) = r = 1 ∑ m X r ( x , θ ) Φ r ( ϕ ) . (9) The in-plane streaming operator is sin ϕ s ( θ ) ⋅ ∇ \sin\phi\,\boldsymbol{s}(\theta)\cdot\nabla sin ϕ s ( θ ) ⋅ ∇ sin ϕ ≥ 0 \sin\phi\ge0 sin ϕ ≥ 0 θ \theta θ cos ϕ ∂ z \cos\phi\,\partial_z cos ϕ ∂ z ϕ \phi ϕ
The macroscopic fields follow from azimuthal contractions of X r X_r X r Φ r \Phi_r Φ r
Y r = ∫ 0 2 π X r d θ , Z r = ∫ 0 2 π s X r d θ , M k ( Φ ) = ∫ 0 π Φ sin k ϕ d ϕ , (10) Y_r = \int_0^{2\pi} X_r\,\mathrm{d}\theta,\qquad Z_r = \int_0^{2\pi} \boldsymbol{s}X_r\,\mathrm{d}\theta,\qquad M_k(\Phi) = \int_0^\pi\Phi\sin^k\phi\,\mathrm{d}\phi, \tag*{(10)} Y r = ∫ 0 2 π X r d θ , Z r = ∫ 0 2 π s X r d θ , M k ( Φ ) = ∫ 0 π Φ sin k ϕ d ϕ , (10) with T r = M 1 ( Φ r ) T_r = M_1(\Phi_r) T r = M 1 ( Φ r ) Q r = M 2 ( Φ r ) Q_r = M_2(\Phi_r) Q r = M 2 ( Φ r )
T m k i n = ∑ r = 1 m Y r T r , q m k i n = ∑ r = 1 m Z r Q r . (11) T_m^{\mathrm{kin}} = \sum_{r=1}^{m} Y_rT_r,\qquad q_m^{\mathrm{kin}} = \sum_{r=1}^{m} Z_rQ_r. \tag*{(11)} T m kin = r = 1 ∑ m Y r T r , q m kin = r = 1 ∑ m Z r Q r . (11) The factor equations below use the product measure d x d θ d ϕ dx\,d\theta\,d\phi d x d θ d ϕ sin ϕ \sin\phi sin ϕ (10) come from d Ω d\Omega d Ω N i j N_{ij} N ij M 1 M_1 M 1 M 3 M_3 M 3 M i j k M_{ijk} M ij k M 2 M_2 M 2 M 4 M_4 M 4 § 4.2 therefore sees each polar factor through four numbers.
Low-rank structure and a Knudsen-uniform bound
The rank of the splitting is set by how ϕ \phi ϕ (3) along the characteristic through x x x Ω \Omega Ω in-plane distance t t t t / sin ϕ t/\sin\phi t / sin ϕ
I ( x , θ , ϕ ) = I b e − δ C d / sin ϕ + ∫ 0 d [ A ( x − t s ) + sin ϕ B ( x − t s , θ ) ] e − δ C t / sin ϕ d t sin ϕ , (12) I(x,\theta,\phi)=I_b e^{-\delta_C d/\sin\phi}+\int_0^d\left[A(x-t s)+\sin\phi\,B(x-t s,\theta)\right]e^{-\delta_C t/\sin\phi}\frac{dt}{\sin\phi}, \tag*{(12)} I ( x , θ , ϕ ) = I b e − δ C d / s i n ϕ + ∫ 0 d [ A ( x − t s ) + sin ϕ B ( x − t s , θ ) ] e − δ C t / s i n ϕ sin ϕ d t , (12) with A = δ C T / ( 4 π ) A=\delta_C T/(4\pi) A = δ C T / ( 4 π ) B = 3 δ N s ⋅ q / ( 4 π ) B=3\delta_N s\cdot q/(4\pi) B = 3 δ N s ⋅ q / ( 4 π ) d ( x , θ ) d(x,\theta) d ( x , θ ) x b x_b x b ( x , θ ) (x,\theta) ( x , θ ) § 2.3 do not depend on ϕ \phi ϕ A A A t t t 1 / sin ϕ 1/\sin\phi 1/ sin ϕ
I ( x , θ , ϕ ) = ∫ [ 0 , ∞ ) κ τ ( ϕ ) d m x , θ ( τ ) , κ τ ( ϕ ) = e − τ / sin ϕ , (13) I(x,\theta,\phi)=\int_{[0,\infty)}\kappa_\tau(\phi)\,dm_{x,\theta}(\tau),\qquad\kappa_\tau(\phi)=e^{-\tau/\sin\phi}, \tag*{(13)} I ( x , θ , ϕ ) = ∫ [ 0 , ∞ ) κ τ ( ϕ ) d m x , θ ( τ ) , κ τ ( ϕ ) = e − τ / s i n ϕ , (13) where τ = δ C t \tau=\delta_C t τ = δ C t m x , θ m_{x,\theta} m x , θ T ( x ) / ( 4 π ) T(x)/(4\pi) T ( x ) / ( 4 π ) τ = 0 \tau=0 τ = 0 I b − T ( x b ) / ( 4 π ) I_b-T(x_b)/(4\pi) I b − T ( x b ) / ( 4 π ) τ = δ C d \tau=\delta_Cd τ = δ C d ( ∂ τ T + 3 χ s ⋅ q ) / ( 4 π ) (\partial_\tau T+3\chi\,s\cdot q)/(4\pi) ( ∂ τ T + 3 χ s ⋅ q ) / ( 4 π ) { κ τ } τ ≥ 0 \{\kappa_\tau\}_{\tau\ge0} { κ τ } τ ≥ 0 p η p_\eta p η p p p ∥ κ τ − ∑ k = 1 p u k ( τ ) v k ∥ ≤ η ∥ κ 0 ∥ \left\|\kappa_\tau-\sum_{k=1}^{p}u_k(\tau)v_k\right\|\le\eta\|\kappa_0\| ∥ κ τ − ∑ k = 1 p u k ( τ ) v k ∥ ≤ η ∥ κ 0 ∥ τ ≥ 0 \tau\ge0 τ ≥ 0 L 2 ( 0 , π ) L^2(0,\pi) L 2 ( 0 , π )
Proposition 1 (Knudsen-uniform separation rank). Let T T T q q q T T T I I I (3) with this gain term and the boundary conditions of § 2.3 . For every η ∈ ( 0 , 1 ) \eta\in(0,1) η ∈ ( 0 , 1 ) X 1 , … , X p η X_1,\ldots,X_{p_\eta} X 1 , … , X p η ( x , θ ) (x,\theta) ( x , θ ) ( x , θ ) (x,\theta) ( x , θ )
∥ I ( x , θ , ⋅ ) − ∑ k = 1 p η X k ( x , θ ) v k ∥ ≤ 2 η ( 1 + ln η − 1 ) ∥ κ 0 ∥ [ ∣ I b ∣ + ∥ T ∥ ∞ + V ( T ) + 3 χ ∥ q ∥ ∞ 4 π ] , (14) \left\|I(x,\theta,\mathord\cdot)-\sum_{k=1}^{p_\eta}X_k(x,\theta)v_k\right\|\le2\eta(1+\ln\eta^{-1})\|\kappa_0\|\left[|I_b|+\frac{\|T\|_\infty+V(T)+3\chi\|q\|_\infty}{4\pi}\right], \tag*{(14)} I ( x , θ , ⋅ ) − k = 1 ∑ p η X k ( x , θ ) v k ≤ 2 η ( 1 + ln η − 1 ) ∥ κ 0 ∥ [ ∣ I b ∣ + 4 π ∥ T ∥ ∞ + V ( T ) + 3 χ ∥ q ∥ ∞ ] , (14) where V ( T ) V(T) V ( T ) T T T v k v_k v k x x x θ \theta θ δ R \delta_R δ R δ N \delta_N δ N p η p_\eta p η
Proof. Since 0 ≤ κ τ ≤ e − τ 0\le\kappa_\tau\le e^{-\tau} 0 ≤ κ τ ≤ e − τ u k = 0 u_k=0 u k = 0 τ ≥ τ η = ln η − 1 \tau\ge\tau_\eta=\ln\eta^{-1} τ ≥ τ η = ln η − 1 κ τ \kappa_\tau κ τ (13) defines X k = ∫ [ 0 , τ η ) u k d m x , θ X_k=\int_{[0,\tau_\eta)}u_k\,dm_{x,\theta} X k = ∫ [ 0 , τ η ) u k d m x , θ ∥ κ 0 ∥ [ η ∣ m x , θ ∣ ( [ 0 , τ η ) ) + ∫ [ τ η , ∞ ) e − τ d ∣ m x , θ ∣ ] \|\kappa_0\|[\eta|m_{x,\theta}|([0,\tau_\eta))+\int_{[\tau_\eta,\infty)}e^{-\tau}\,d|m_{x,\theta}|] ∥ κ 0 ∥ [ η ∣ m x , θ ∣ ([ 0 , τ η )) + ∫ [ τ η , ∞ ) e − τ d ∣ m x , θ ∣ ] m x , θ m_{x,\theta} m x , θ (14) . Without the cut-off at τ η \tau_\eta τ η ∫ 0 d ∣ B ∣ d t \int_0^d|B|\,dt ∫ 0 d ∣ B ∣ d t m x , θ m_{x,\theta} m x , θ δ N \delta_N δ N
The proposition covers the exact transport solution for any given source, including the frozen-source problems of PGD–GSIS. For the DG solution ϕ \phi ϕ δ C / sin ϕ \delta_C/\sin\phi δ C / sin ϕ
( sin ϕ G θ + δ C ) − 1 = 1 sin ϕ ∫ 0 ∞ e − δ C t / sin ϕ e − t G θ d t , (\sin\phi\,G_\theta+\delta_C)^{-1}=\frac{1}{\sin\phi}\int_0^\infty e^{-\delta_C t/\sin\phi}e^{-tG_\theta}\,\mathrm{d}t, ( sin ϕ G θ + δ C ) − 1 = sin ϕ 1 ∫ 0 ∞ e − δ C t / s i n ϕ e − t G θ d t , where the DG streaming matrix G θ G_\theta G θ
The kernel rank is small in practice. A weighted SVD of κ τ \kappa_\tau κ τ τ = 0 \tau=0 τ = 0 [ 10 − 8 , 10 3 ] [10^{-8},10^3] [ 1 0 − 8 , 1 0 3 ] p η ≈ 5 p_\eta\approx5 p η ≈ 5 η = 10 − 2 \eta=10^{-2} η = 1 0 − 2 10 − 4 10^{-4} 1 0 − 4 10 − 6 10^{-6} 1 0 − 6 p η p_\eta p η
Rank in the transport regimes
Proposition 1 holds in every regime. How much of the bound p η p_\eta p η m x , θ m_{\mathbf{x},\theta} m x , θ
Diffusive and hydrodynamic regimes. In the collision-dominated bulk, repeated integration by parts turns Eq. (12) into a series in sin ϕ / δ C \sin\phi/\delta_C sin ϕ / δ C I = E − K n C Ω ⋅ ∇ E + O ( K n C 2 ) I=E-\mathrm{Kn}_C\boldsymbol{\Omega}\cdot\nabla E+O(\mathrm{Kn}_C^2) I = E − Kn C Ω ⋅ ∇ E + O ( Kn C 2 ) E = ( T + 3 χ Ω ⋅ q ) / ( 4 π ) E=(T+3\chi\boldsymbol{\Omega}\cdot\mathbf{q})/(4\pi) E = ( T + 3 χ Ω ⋅ q ) / ( 4 π )
Proposition 2. Let f ( Ω ) f(\boldsymbol{\Omega}) f ( Ω ) N N N Ω \boldsymbol{\Omega} Ω Ω z = cos ϕ \Omega_z=\cos\phi Ω z = cos ϕ f = ∑ k = 1 N + 1 a k ( θ ) g k ( ϕ ) f=\sum_{k=1}^{N+1}a_k(\theta)g_k(\phi) f = ∑ k = 1 N + 1 a k ( θ ) g k ( ϕ ) { g k } \{g_k\} { g k }
G N = span { cos 2 j ϕ : 2 j ≤ N } ⊕ sin ϕ ⋅ span { cos 2 j ϕ : 2 j ≤ N − 1 } , dim G N = N + 1. G_N=\operatorname{span}\{\cos^{2j}\phi:2j\le N\}\oplus\sin\phi\cdot\operatorname{span}\{\cos^{2j}\phi:2j\le N-1\},\qquad\dim G_N=N+1. G N = span { cos 2 j ϕ : 2 j ≤ N } ⊕ sin ϕ ⋅ span { cos 2 j ϕ : 2 j ≤ N − 1 } , dim G N = N + 1. Consequently a distribution whose angular dependence is, at every x x x N + 1 N+1 N + 1 ( N + 1 ) ( N + 2 ) / 2 (N+1)(N+2)/2 ( N + 1 ) ( N + 2 ) /2
Proof. A monomial Ω x a Ω y b Ω z c \Omega_x^a\Omega_y^b\Omega_z^c Ω x a Ω y b Ω z c c c c cos a θ sin b θ ⋅ sin n ϕ cos c ϕ \cos^a\theta\sin^b\theta\cdot\sin^n\phi\cos^c\phi cos a θ sin b θ ⋅ sin n ϕ cos c ϕ n = a + b n=a+b n = a + b n n n sin n ϕ = ( 1 − cos 2 ϕ ) n / 2 \sin^n\phi=(1-\cos^2\phi)^{n/2} sin n ϕ = ( 1 − cos 2 ϕ ) n /2 cos ϕ \cos\phi cos ϕ n + c ≤ N n+c\le N n + c ≤ N n n n sin ϕ \sin\phi sin ϕ n − 1 + c ≤ N − 1 n-1+c\le N-1 n − 1 + c ≤ N − 1 G N G_N G N ⌊ N / 2 ⌋ + ⌊ ( N − 1 ) / 2 ⌋ + 2 = N + 1 \lfloor N/2\rfloor+\lfloor(N-1)/2\rfloor+2=N+1 ⌊ N /2 ⌋ + ⌊( N − 1 ) /2 ⌋ + 2 = N + 1 G N G_N G N
For N = 0 N=0 N = 0 T / ( 4 π ) = X ⋅ 1 T/(4\pi)=X\cdot1 T / ( 4 π ) = X ⋅ 1 X = T / ( 4 π ) X=T/(4\pi) X = T / ( 4 π ) N = 1 N=1 N = 1 P 1 P_1 P 1
T 4 π + 3 Ω ⋅ q 4 π = T ( x ) 4 π ⏟ X 1 ⋅ 1 ⏟ Φ 1 + s ⋅ q ( x ) ⏟ X 2 ⋅ 3 sin ϕ 4 π ⏟ Φ 2 . (15) \frac{T}{4\pi}+\frac{3\boldsymbol{\Omega}\cdot\mathbf{q}}{4\pi}
=
\underbrace{\frac{T(x)}{4\pi}}_{X_1}\cdot\underbrace{1}_{\Phi_1}
+
\underbrace{\mathbf{s}\cdot\mathbf{q}(x)}_{X_2}\cdot\underbrace{\frac{3\sin\phi}{4\pi}}_{\Phi_2}. \tag*{(15)} 4 π T + 4 π 3 Ω ⋅ q = X 1 4 π T ( x ) ⋅ Φ 1 1 + X 2 s ⋅ q ( x ) ⋅ Φ 2 4 π 3 sin ϕ . (15) The azimuthal integrals ∫ s ( s ⋅ q ) d θ = π q \int\mathbf{s}(\mathbf{s}\cdot\mathbf{q})\,\mathrm{d}\theta=\pi\mathbf{q} ∫ s ( s ⋅ q ) d θ = π q M 2 ( Φ 2 ) = 1 / π M_2(\Phi_2)=1/\pi M 2 ( Φ 2 ) = 1/ π (15) returns T T T q \mathbf{q} q (11) . The two collision-dominated regimes differ in the size of the flux. In the diffusive regime q = O ( K n R ) \mathbf{q}=O(\mathrm{Kn}_{R}) q = O ( Kn R ) P 1 P_1 P 1 q = O ( 1 ) \mathbf{q}=O(1) q = O ( 1 ) − 3 χ K n C sin 2 ϕ s ⋅ ∇ ( s ⋅ q ) / ( 4 π ) -3\chi\mathrm{Kn}_{C}\sin^2\phi\,\mathbf{s}\cdot\nabla(\mathbf{s}\cdot\mathbf{q})/(4\pi) − 3 χ Kn C sin 2 ϕ s ⋅ ∇ ( s ⋅ q ) / ( 4 π ) P 2 P_2 P 2
Ballistic regime. When δ C d → 0 \delta_Cd\to0 δ C d → 0 τ = 0 \tau=0 τ = 0 I I I I b I_b I b ϕ \phi ϕ Φ 1 = 1 \Phi_1=1 Φ 1 = 1 κ 0 \kappa_0 κ 0 sin ϕ ≲ τ \sin\phi\lesssim\tau sin ϕ ≲ τ ∥ κ τ − κ 0 ∥ / ∥ κ 0 ∥ ≈ ( 4 ln 2 / π ) 1 / 2 τ 1 / 2 \lVert\kappa_\tau-\kappa_0\rVert/\lVert\kappa_0\rVert\approx(4\ln2/\pi)^{1/2}\tau^{1/2} ∥ κ τ − κ 0 ∥ / ∥ κ 0 ∥ ≈ ( 4 ln 2/ π ) 1/2 τ 1/2 τ \tau τ § 6.3 ).
Transition regime and wall layers. At intermediate Knudsen numbers, and in the wall layers and corners of the collision-dominated regimes, the measure spreads over optical depths of order one. Neither expansion applies, and the rank is limited only by the kernel rank p η p_\eta p η § 6.3 confirms this ordering. A similar low-rank structure in the hydrodynamic and free-molecular limits was noted for the Shakhov model in [30 ] ; Proposition 1 extends it to the transition regime with an explicit bound.
PGD–GSIS construction and analysis
The factors are computed greedily: modes 1 , … , m − 1 1,\ldots,m-1 1 , … , m − 1 ( X m , Φ m ) (X_m,\Phi_m) ( X m , Φ m ) Pure PGD is the Galerkin PGD of the complete linear equation (3) : the gain term is evaluated from the kinetic moments of the current separated state, implicitly in the polar update and lagged by one sweep in the spatial update. PGD–GSIS evaluates the gain term from a macroscopic source U s = ( T s , q s ) U_s=(T_s,\mathbf{q}_s) U s = ( T s , q s ) (6) . Section 4.1 derives the factor equations of both solvers, Section 4.2 gives the PGD–GSIS algorithm, Section 4.3 analyzes the convergence of both iterations, and Section 4.4 states a conditional AP property of PGD–GSIS.
Factor equations of pure PGD and PGD–GSIS
Let ( a , b ) x θ = ∫ D ∫ 0 2 π a b d θ d x (a,b)_{x\theta}=\int_{\mathcal{D}}\int_0^{2\pi}ab\,\mathrm{d}\theta\,\mathrm{d}x ( a , b ) x θ = ∫ D ∫ 0 2 π ab d θ d x b h ( V , U ) b_h(V,U) b h ( V , U ) ( T , q ) (T,\mathbf{q}) ( T , q ) (3) . Testing Eq. (3) with X m Φ ∗ X_m\Phi^* X m Φ ∗ ϕ \phi ϕ
( δ C α m m + β m m sin ϕ ) Φ m = F T , m + sin ϕ ( F q , m − μ m ) − ∑ r = 1 m − 1 ( δ C α m r + β m r sin ϕ ) Φ r , (16) (\delta_C\alpha_{mm}+\beta_{mm}\sin\phi)\Phi_m=F_{T,m}+\sin\phi\,(F_{q,m}-\mu_m)-\sum_{r=1}^{m-1}(\delta_C\alpha_{mr}+\beta_{mr}\sin\phi)\Phi_r, \tag*{(16)} ( δ C α mm + β mm sin ϕ ) Φ m = F T , m + sin ϕ ( F q , m − μ m ) − r = 1 ∑ m − 1 ( δ C α m r + β m r sin ϕ ) Φ r , (16) with α m r = ( X m , X r ) x θ \alpha_{mr}=(X_m,X_r)_{x\theta} α m r = ( X m , X r ) x θ β m r = b h ( X m , X r ) \beta_{mr}=b_h(X_m,X_r) β m r = b h ( X m , X r ) F T , m = δ C ( X m , T ) x θ / ( 4 π ) F_{T,m}=\delta_C(X_m,T)_{x\theta}/(4\pi) F T , m = δ C ( X m , T ) x θ / ( 4 π ) F q , m = 3 δ N ( X m , s ⋅ q ) x θ / ( 4 π ) F_{q,m}=3\delta_N(X_m,\boldsymbol{s}\cdot\boldsymbol{q})_{x\theta}/(4\pi) F q , m = 3 δ N ( X m , s ⋅ q ) x θ / ( 4 π ) μ m \mu_m μ m X m X_m X m X ∗ Φ m X^*\Phi_m X ∗ Φ m
a m m X m + b m m s ⋅ ∇ h X m = c m T + d m s ⋅ q − ∑ r = 1 m − 1 ( a m r X r + b m r s ⋅ ∇ h X r ) , (17) a_{mm}X_m+b_{mm}\boldsymbol{s}\cdot\nabla_hX_m=c_mT+d_m\boldsymbol{s}\cdot\boldsymbol{q}-\sum_{r=1}^{m-1}(a_{mr}X_r+b_{mr}\boldsymbol{s}\cdot\nabla_hX_r), \tag*{(17)} a mm X m + b mm s ⋅ ∇ h X m = c m T + d m s ⋅ q − r = 1 ∑ m − 1 ( a m r X r + b m r s ⋅ ∇ h X r ) , (17) a m r = δ C ∫ 0 π Φ m Φ r d ϕ , b m r = ∫ 0 π sin ϕ Φ m Φ r d ϕ , c m = δ C 4 π ∫ 0 π Φ m d ϕ , d m = 3 δ N 4 π T m , (18) a_{mr}=\delta_C\int_0^\pi\Phi_m\Phi_r\,\mathrm{d}\phi,\qquad b_{mr}=\int_0^\pi\sin\phi\,\Phi_m\Phi_r\,\mathrm{d}\phi,\qquad c_m=\frac{\delta_C}{4\pi}\int_0^\pi\Phi_m\,\mathrm{d}\phi,\qquad d_m=\frac{3\delta_N}{4\pi}T_m, \tag*{(18)} a m r = δ C ∫ 0 π Φ m Φ r d ϕ , b m r = ∫ 0 π sin ϕ Φ m Φ r d ϕ , c m = 4 π δ C ∫ 0 π Φ m d ϕ , d m = 4 π 3 δ N T m , (18) where ∇ h \nabla_h ∇ h ∑ r ≤ m b m r X r = T m T b / ( 4 π ) \sum_{r\le m}b_{mr}X_r=T_mT_b/(4\pi) ∑ r ≤ m b m r X r = T m T b / ( 4 π ) θ j \theta_j θ j n p × n p n_p\times n_p n p × n p n p = ( p + 1 ) ( p + 2 ) / 2 n_p=(p+1)(p+2)/2 n p = ( p + 1 ) ( p + 2 ) /2 X m X_m X m N s w N_{\mathrm{sw}} N sw X m ← X m / ∥ X m ∥ X_m\leftarrow X_m/\lVert X_m\rVert X m ← X m / ∥ X m ∥ Φ m ← ∥ X m ∥ Φ m \Phi_m\leftarrow\lVert X_m\rVert\Phi_m Φ m ← ∥ X m ∥ Φ m ( T , q ) (T,\boldsymbol{q}) ( T , q ) (16) and (17) .
Pure PGD. The gain term is evaluated from the kinetic moments (11) of the current state, ( T , q ) = K [ I m ] = ( T m k i n , q m k i n ) (T,\boldsymbol{q})=\mathcal{K}[I_m]=(T_m^{\mathrm{kin}},\boldsymbol{q}_m^{\mathrm{kin}}) ( T , q ) = K [ I m ] = ( T m kin , q m kin ) α m r T = ( X m , Y r ) x θ / ( 4 π ) \alpha_{mr}^{T}=(X_m,Y_r)_{x\theta}/(4\pi) α m r T = ( X m , Y r ) x θ / ( 4 π ) α m r q = 3 ( X m , s ⋅ Z r ) x θ / ( 4 π ) \alpha_{mr}^{q}=3(X_m,\boldsymbol{s}\cdot Z_r)_{x\theta}/(4\pi) α m r q = 3 ( X m , s ⋅ Z r ) x θ / ( 4 π ) F T , m = δ C ∑ r ≤ m α m r T T r F_{T,m}=\delta_C\sum_{r\le m}\alpha_{mr}^{T}T_r F T , m = δ C ∑ r ≤ m α m r T T r F q , m = δ N ∑ r ≤ m α m r q Q r F_{q,m}=\delta_N\sum_{r\le m}\alpha_{mr}^{q}Q_r F q , m = δ N ∑ r ≤ m α m r q Q r r = m r=m r = m T m = M 1 ( Φ m ) T_m=M_1(\Phi_m) T m = M 1 ( Φ m ) Q m = M 2 ( Φ m ) Q_m=M_2(\Phi_m) Q m = M 2 ( Φ m ) (16) becomes
( δ C α m m + β m m sin ϕ ) Φ m − δ C α m m T M 1 ( Φ m ) − δ N α m m q sin ϕ M 2 ( Φ m ) = ∑ r = 1 m − 1 [ δ C α m r T T r + δ N α m r q sin ϕ Q r − ( δ C α m r + β m r sin ϕ ) Φ r ] − μ m sin ϕ . (19) \begin{aligned}
(\delta_C\alpha_{mm}+\beta_{mm}\sin\phi)\Phi_m-\delta_C\alpha_{mm}^{T}M_1(\Phi_m)-\delta_N\alpha_{mm}^{q}\sin\phi\,M_2(\Phi_m)\\
=\sum_{r=1}^{m-1}\left[\delta_C\alpha_{mr}^{T}T_r+\delta_N\alpha_{mr}^{q}\sin\phi\,Q_r-(\delta_C\alpha_{mr}+\beta_{mr}\sin\phi)\Phi_r\right]-\mu_m\sin\phi. \tag*{(19)}
\end{aligned} ( δ C α mm + β mm sin ϕ ) Φ m − δ C α mm T M 1 ( Φ m ) − δ N α mm q sin ϕ M 2 ( Φ m ) = r = 1 ∑ m − 1 [ δ C α m r T T r + δ N α m r q sin ϕ Q r − ( δ C α m r + β m r sin ϕ ) Φ r ] − μ m sin ϕ . (19) The polar nodes are coupled through the moments M 1 M_1 M 1 M 2 M_2 M 2 (19) is an N ϕ × N ϕ N_\phi\times N_\phi N ϕ × N ϕ (17) with ( T , q ) = K [ I m ] (T,\boldsymbol{q})=\mathcal{K}[I_m] ( T , q ) = K [ I m ] Y m T m Y_mT_m Y m T m Z m Q m Z_mQ_m Z m Q m
PGD–GSIS. The gain term is evaluated from the frozen source, ( T , q ) = U s (T,\boldsymbol{q})=U_s ( T , q ) = U s § 4.2 ). The N ϕ N_\phi N ϕ (16) then decouple and the polar update is a division. It is well defined: upwinding makes b h ( X , X ) ≥ 0 b_h(X,X)\ge0 b h ( X , X ) ≥ 0 δ C α m m > 0 \delta_C\alpha_{mm}>0 δ C α mm > 0 (17) is fixed during the update; repeated sweeps only update the lagged periodic and adiabatic traces.
PGD–GSIS coupling and algorithm
PGD–GSIS couples each polar/spatial pair to a synthetic solve that refreshes U s U_s U s I m I_m I m
U s ← G [ I m ] = ( T s y n , q s y n ) , (21) U_s \leftarrow\mathcal{G}[I_m] = (T^{\mathrm{syn}},\boldsymbol{q}^{\mathrm{syn}}), \tag*{(21)} U s ← G [ I m ] = ( T syn , q syn ) , (21) where the synthetic map G \mathcal{G} G (6) with H \boldsymbol{H} H I m I_m I m K \mathcal{K} K (11) , is used instead inside the factor updates of pure PGD; both solvers report K [ I m ] \mathcal{K}[I_m] K [ I m ]
High-order moments from the PGD factors. The synthetic system needs H i j = ∂ k M i j k H_{ij} = \partial_k M_{ijk} H ij = ∂ k M ij k v h v_h v h K K K
( v h , H i j ) K = − ( ∂ k v h , ∫ K K i j k I m d Ω ) K + ⟨ v h , ∫ ∂ K n k K i j k I m u p d Ω ⟩ ∂ K , (22) (v_h,H_{ij})_K = -\left(\partial_k v_h,\int_K K_{ijk}I_m\,\mathrm{d}\Omega\right)_K + \left\langle v_h,\int_{\partial K} n_k K_{ijk}I_m^{\mathrm{up}}\,\mathrm{d}\Omega\right\rangle_{\partial K}, \tag*{(22)} ( v h , H ij ) K = − ( ∂ k v h , ∫ K K ij k I m d Ω ) K + ⟨ v h , ∫ ∂ K n k K ij k I m up d Ω ⟩ ∂ K , (22) where I m u p I_m^{\mathrm{up}} I m up ∫ K i j k X r Φ r d Ω = ∫ 0 2 π X r [ M 4 ( Φ r ) s i s j s k − 1 5 M 2 ( Φ r ) ( s i δ j k + s j δ i k + s k δ i j ) ] d θ \int K_{ijk}X_r\Phi_r\,\mathrm{d}\Omega= \int_0^{2\pi} X_r\left[M_4(\Phi_r)s_i s_j s_k - \frac{1}{5}M_2(\Phi_r)(s_i\delta_{jk} + s_j\delta_{ik} + s_k\delta_{ij})\right]\,\mathrm{d}\theta ∫ K ij k X r Φ r d Ω = ∫ 0 2 π X r [ M 4 ( Φ r ) s i s j s k − 5 1 M 2 ( Φ r ) ( s i δ j k + s j δ ik + s k δ ij ) ] d θ θ \theta θ O ( m N θ N x ) O(mN_\theta N_x) O ( m N θ N x ) N ϕ N_\phi N ϕ N x = N e l n p N_x = N_{\mathrm{el}}n_p N x = N el n p
Synthetic solve and wall closure. Equation [6 ] is written as a first-order system for ( T , q ) (T,\boldsymbol{q}) ( T , q ) [37 , 38 , 39 ] on the kinetic mesh with the same degree (Appendix § 11 ). Element unknowns are eliminated locally and the global problem involves only the traces ( T ^ , q ^ ) (\widehat{T},\widehat{\boldsymbol{q}}) ( T , q ) F ‾ ( U ) = F k i n [ I m ] + J F ( U − U k i n ) \overline{\boldsymbol{F}}(U) = \boldsymbol{F}^{\mathrm{kin}}[I_m] + J_F(U-U^{\mathrm{kin}}) F ( U ) = F kin [ I m ] + J F ( U − U kin ) F k i n \boldsymbol{F}^{\mathrm{kin}} F kin I m I_m I m M 2 M_2 M 2 M 3 M_3 M 3 J F J_F J F [19 , 21 ] . The synthetic fields are the first three element components of the solution.
Algorithm 1 Fixed-parameter PGD–GSIS
Input: mesh, angular quadrature, boundary data, ( δ R , δ N ) (\delta_R,\delta_N) ( δ R , δ N )
Assemble DG data; factorize the local and global HDG systems and the wall Jacobians.
Set U s = ( T s , q s ) = 0 U_s=(T_s,q_s)=0 U s = ( T s , q s ) = 0
for m = 1 , … , M max m=1,\ldots,M_{\max} m = 1 , … , M m a x do
Initialize X m X_m X m
repeat
Polar update (16) with U s U_s U s M k ( Φ m ) M_k(\Phi_m) M k ( Φ m ) k = 1 , … , 4 k=1,\ldots,4 k = 1 , … , 4 (18) .
Spatial update (17) with U s U_s U s
Normalize the pair.
Synthetic refresh from the complete retained state: U s ← G [ I m ] U_s\leftarrow\mathcal{G}[I_m] U s ← G [ I m ]
until the amplitude test holds or L max L_{\max} L m a x
Retain ( X m , Φ m ) (X_m,\Phi_m) ( X m , Φ m ) r e n ( m ) < ϵ e n r_{\mathrm{en}}^{(m)}<\epsilon_{\mathrm{en}} r en ( m ) < ϵ en m > 1 m>1 m > 1
end for
Output: factors { X r , Φ r } \{X_r,\Phi_r\} { X r , Φ r } K [ I m ] \mathcal{K}[I_m] K [ I m ] G [ I m ] \mathcal{G}[I_m] G [ I m ]
Algorithm and stopping criteria. Algorithm 1 has three nested loops, and none of its tests uses a reference solution. Spatial sweeps stop on the relative change of X m X_m X m ∥ ⋅ ∥ h \lVert\cdot\rVert_h ∥ ⋅ ∥ h ϵ s w \epsilon_{\mathrm{sw}} ϵ sw N s w N_{\mathrm{sw}} N sw ∥ Y m ∥ h ∣ T m ∣ \lVert Y_m\rVert_h\lvert T_m\rvert ∥ Y m ∥ h ∣ T m ∣ ∥ Z m , i ∥ h ∣ Q m ∣ \lVert Z_{m,i}\rVert_h\lvert Q_m\rvert ∥ Z m , i ∥ h ∣ Q m ∣ ϵ A L S \epsilon_{\mathrm{ALS}} ϵ ALS L max L_{\max} L m a x
r e n ( m ) = ∥ W m ∥ h ∥ ∑ r ≤ m W r ∥ h < ϵ e n , W r = ( Y r T r , Z r Q r ) , (23) r_{\mathrm{en}}^{(m)}=\frac{\lVert W_m\rVert_h}{\left\lVert\sum_{r\leq m}W_r\right\rVert_h}<\epsilon_{\mathrm{en}},\qquad W_r=(Y_rT_r,Z_rQ_r), \tag*{(23)} r en ( m ) = ∑ r ≤ m W r h ∥ W m ∥ h < ϵ en , W r = ( Y r T r , Z r Q r ) , (23) for some m > 1 m>1 m > 1 M max M_{\max} M m a x (19) , and the spatial sweeps lag the current-mode moments; there is no refresh step.
Storage and computational cost. For a single parameter the distribution needs N x N θ N ϕ N_xN_\theta N_\phi N x N θ N ϕ m ( N x N θ + N ϕ ) m(N_xN_\theta+N_\phi) m ( N x N θ + N ϕ ) N x = 2000 N_x=2000 N x = 2000 N θ = 80 N_\theta=80 N θ = 80 N ϕ = 40 N_\phi=40 N ϕ = 40 m = 30 m=30 m = 30 6.4 × 10 6 6.4\times10^6 6.4 × 1 0 6 4.8 × 10 6 4.8\times10^6 4.8 × 1 0 6 N ϕ N_\phi N ϕ N δ N_\delta N δ m ( N x N θ + N ϕ N δ ) m(N_xN_\theta+N_\phi N_\delta) m ( N x N θ + N ϕ N δ ) N x N θ N ϕ N δ N_xN_\theta N_\phi N_\delta N x N θ N ϕ N δ N δ = 441 N_\delta=441 N δ = 441 5.3 × 10 6 5.3 \times10^{6} 5.3 × 1 0 6 2.8 × 10 9 2.8 \times10^{9} 2.8 × 1 0 9 N ϕ N_{\phi} N ϕ N δ N_{\delta} N δ § 5.1 ).
Per alternating pair, the spatial update costs O ( N θ N e l n p 3 ) O(N_{\theta}N_{\mathrm{el}}n_{p}^{3}) O ( N θ N el n p 3 ) O ( N θ N e l n p 2 ) O(N_{\theta}N_{\mathrm{el}}n_{p}^{2}) O ( N θ N el n p 2 ) O ( m N θ N x ) O(mN_{\theta}N_{x}) O ( m N θ N x ) N ϕ N_{\phi} N ϕ O ( m N ϕ ) O(mN_{\phi}) O ( m N ϕ ) N ϕ × N ϕ N_{\phi} \times N_{\phi} N ϕ × N ϕ O ( m N θ N x ) O(mN_{\theta}N_{x}) O ( m N θ N x ) O ( N θ N ϕ N e l n p 2 ) O(N_{\theta}N_{\phi}N_{\mathrm{el}}n_{p}^{2}) O ( N θ N ϕ N el n p 2 ) N ϕ N_{\phi} N ϕ N ϕ = 40 N_{\phi}=40 N ϕ = 40 § 6.4 ).
Convergence and kinetic–synthetic consistency
Consider a Fourier error mode e i k x e^{\mathrm{i}kx} e i k x § 9 ). Define S U = ( δ C T + 3 δ N Ω ⋅ q ) / ( 4 π ) SU = (\delta_{C}T + 3\delta_{N}\boldsymbol{\Omega} \cdot\boldsymbol{q})/(4\pi) S U = ( δ C T + 3 δ N Ω ⋅ q ) / ( 4 π ) T = Ω ⋅ ∇ + δ C \mathcal{T} = \boldsymbol{\Omega} \cdot\nabla+ \delta_{C} T = Ω ⋅ ∇ + δ C A M = M ∘ T − 1 ∘ S A_{M} = \mathcal{M} \circ\mathcal{T}^{-1} \circ\mathcal{S} A M = M ∘ T − 1 ∘ S M ∈ { K , G } \mathcal{M} \in\{\mathcal{K},\mathcal{G}\} M ∈ { K , G } [21 ] ; they serve as references for the two PGD iterations.
Pure PGD. The implicit polar update (19) balances the current mode against its own moments only as seen through its fixed spatial shape X m X_{m} X m β m m = b h ( X m , X m ) \beta_{mm} = b_{h}(X_{m}, X_{m}) β mm = b h ( X m , X m ) § 10 ) and vanishes for a smooth mode in a periodic bulk. The polar update therefore adjusts the angular shape and the amplitude of the mode, but it does not move macroscopic information across the domain. Only the spatial update does, and there the moments of the current mode are lagged. For a fixed polar factor Φ \Phi Φ (17) is an in-plane transport problem with extinction a / b a/b a / b c T Y + d Q s ⋅ Z cT Y + dQ \boldsymbol{s} \cdot Z c T Y + d Q s ⋅ Z e i k x e^{\mathrm{i}kx} e i k x ( Y ^ , Z ^ ) (\widehat{Y},\widehat{Z}) ( Y , Z ) ρ s w \rho_{\mathrm{sw}} ρ sw § 9 ). Table 1 lists 1 − ρ s w 1-\rho_{\mathrm{sw}} 1 − ρ sw Φ = 1 \Phi=1 Φ = 1 Φ = sin φ \Phi=\sin\varphi Φ = sin φ 1 − ρ ( A K ) 1-\rho(A_{\mathcal{K}}) 1 − ρ ( A K ) K n C K n R k 2 / 3 \mathrm{Kn}_{C}\mathrm{Kn}_{R}k^{2}/3 Kn C Kn R k 2 /3 K n C ( δ R + K n C k 2 / 5 ) \mathrm{Kn}_{C}(\delta_{R}+\mathrm{Kn}_{C}k^{2}/5) Kn C ( δ R + Kn C k 2 /5 ) [21 , 12 ] . The temperature factor is the slower one at the diffusive pairs; the flux factor is the slower one at the hydrodynamic pairs, where normal scattering nearly conserves momentum.
Reducing a slow error component by one decade therefore takes ln 10 / ( 1 − ρ s w ) \ln10/(1-\rho_{\mathrm{sw}}) ln 10/ ( 1 − ρ sw ) 3 × 10 4 3 \times10^{4} 3 × 1 0 4 10 4 10^{4} 1 0 4 8 × 10 3 8 \times10^{3} 8 × 1 0 3 ( 100 , 0.01 ) (100, 0.01) ( 100 , 0.01 ) N s w L max = 50 N_{\mathrm{sw}}L_{\max}=50 N sw L m a x = 50 1 − ρ s w 1-\rho_{\mathrm{sw}} 1 − ρ sw ϵ s w = 10 − 3 \epsilon_{\mathrm{sw}}=10^{-3} ϵ sw = 1 0 − 3 1 − ρ s w < 10 − 3 1-\rho_{\mathrm{sw}}<10^{-3} 1 − ρ sw < 1 0 − 3 ℓ ≤ 1 \ell\le1 ℓ ≤ 1 K \mathcal{K} K Proposition 3 (i) below). The slow components are close to P 1 P_1 P 1 1 − ρ s w ≈ 0.99 1-\rho_{\mathrm{sw}}\approx0.99 1 − ρ sw ≈ 0.99
PGD–GSIS. Let U n U^n U n n n n I e x n = T − 1 S U n I_{\mathrm{ex}}^n=\mathcal{T}^{-1}\mathcal{S}U^n I ex n = T − 1 S U n e n = I m n − I e x n e^n=I_{\mathrm{m}}^n-I_{\mathrm{ex}}^n e n = I m n − I ex n § 9 ) gives
U n + 1 = A G U n + G [ e n ] . (24) U^{n+1}=\mathcal{A}_{\mathcal{G}}U^n+\mathcal{G}[e^n]. \tag*{(24)} U n + 1 = A G U n + G [ e n ] . (24) The same identity with K \mathcal{K} K G \mathcal{G} G § 5.1 ). The effect of PGD errors on PGD–GSIS is therefore determined by their image under the synthetic map.
Proposition 3. (i) K [ e ] = ( ∫ e d Ω , ∫ Ω e d Ω ) \mathcal{K}[e]=\left(\int e\,\mathrm{d}\Omega,\int\Omega e\,\mathrm{d}\Omega\right) K [ e ] = ( ∫ e d Ω , ∫ Ω e d Ω ) e e e ℓ ≤ 1 \ell\le1 ℓ ≤ 1
(ii) In a periodic domain, G [ e ] \mathcal{G}[e] G [ e ] e e e M i j k [ e ] = ∫ K i j k e d Ω M_{ijk}[e]=\int K_{ijk}e\,\mathrm{d}\Omega M ij k [ e ] = ∫ K ij k e d Ω K i j k K_{ijk} K ij k Ω ⊗ Ω ⊗ Ω \Omega\otimes\Omega\otimes\Omega Ω ⊗ Ω ⊗ Ω ℓ = 3 \ell=3 ℓ = 3 e e e e i k x e^{ikx} e ik x
G [ e ] = ( 3 i k K n C M x x x [ e ] , 0 , − K n C k 2 δ R + K n C k 2 / 5 M y x x [ e ] ) . (25) \mathcal{G}[e]=\left(3ik\mathrm{Kn}_{\mathrm{C}}M_{xxx}[e],\,0,\,-\frac{\mathrm{Kn}_{\mathrm{C}}k^2}{\delta_{\mathrm{R}}+\mathrm{Kn}_{\mathrm{C}}k^2/5}M_{yxx}[e]\right). \tag*{(25)} G [ e ] = ( 3 ik Kn C M xxx [ e ] , 0 , − δ R + Kn C k 2 /5 Kn C k 2 M y xx [ e ] ) . (25) Proof. Part (i) is the definition (11). For (ii), G \mathcal{G} G H = ∇ ⋅ M H=\nabla\cdot M H = ∇ ⋅ M ∇ → i k e x \nabla\to ik e_x ∇ → ik e x (6) gives q ^ x = 0 \widehat{q}_x=0 q x = 0 ∇ ⋅ q = 0 \nabla\cdot q=0 ∇ ⋅ q = 0 (25) from the two momentum equations. □
Since K n C k 2 / ( δ R + K n C k 2 / 5 ) = 5 ( ℓ G K k ) 2 / [ 1 + ( ℓ G K k ) 2 ] \mathrm{Kn}_{\mathrm{C}}k^2/(\delta_{\mathrm{R}}+\mathrm{Kn}_{\mathrm{C}}k^2/5)=5(\ell_{\mathrm{GK}}k)^2/[1+(\ell_{\mathrm{GK}}k)^2] Kn C k 2 / ( δ R + Kn C k 2 /5 ) = 5 ( ℓ GK k ) 2 / [ 1 + ( ℓ GK k ) 2 ] ℓ = 3 \ell=3 ℓ = 3 3 k K n C 3k\mathrm{Kn}_{\mathrm{C}} 3 k Kn C k K n C ≪ 1 k\mathrm{Kn}_{\mathrm{C}}\ll1 k Kn C ≪ 1 ℓ G K k ≪ 1 \ell_{\mathrm{GK}}k\ll1 ℓ GK k ≪ 1 ℓ = 3 \ell=3 ℓ = 3 § 4.2 adds half-range moment inputs with no explicit K n C \mathrm{Kn}_{\mathrm{C}} Kn C γ G n = ∥ G [ e n ] ∥ / ∥ U n ∥ \gamma_{\mathcal{G}}^n=\|\mathcal{G}[e^n]\|/\|U^n\| γ G n = ∥ G [ e n ] ∥/∥ U n ∥ (24) gives ∥ U n + 1 ∥ ≤ ( ∥ A G ∥ + γ G n ) ∥ U n ∥ \|U^{n+1}\|\le(\|\mathcal{A}_{\mathcal{G}}\|+\gamma_{\mathcal{G}}^n)\|U^n\| ∥ U n + 1 ∥ ≤ ( ∥ A G ∥ + γ G n ) ∥ U n ∥ Table 1 describe the observed iteration but do not establish a uniform contraction bound.
Table 1. Fourier analysis for an error mode e i π x e^{i\pi x} e iπ x N θ × N ϕ = 80 × 40 N_{\theta}\times N_{\phi}=80\times40 N θ × N ϕ = 80 × 40 1 − ρ ( A K ) 1-\rho(A_{K}) 1 − ρ ( A K ) 1 − ρ s w T 1-\rho_{\mathrm{sw}}^{T} 1 − ρ sw T 1 − ρ s w q 1-\rho_{\mathrm{sw}}^{q} 1 − ρ sw q Φ = 1 \Phi=1 Φ = 1 Φ = sin ϕ \Phi=\sin\phi Φ = sin ϕ ρ ( A G ) \rho(A_{G}) ρ ( A G ) ∥ e ∥ \lVert e\rVert ∥ e ∥ ∥ I m n − I e x n ∥ / ∥ I e x n ∥ \lVert I_{m}^{n}-I_{\mathrm{ex}}^{n}\rVert/\lVert I_{\mathrm{ex}}^{n}\rVert ∥ I m n − I ex n ∥ / ∥ I ex n ∥ M max = 20 M_{\max}=20 M m a x = 20 L max = 10 L_{\max}=10 L m a x = 10 γ G \gamma_{G} γ G ∥ G [ e n ] ∥ / ∥ U n ∥ \lVert G[e^{n}]\rVert/\lVert U^{n}\rVert ∥ G [ e n ]∥ / ∥ U n ∥ n 12 n_{12} n 12 10 12 10^{12} 1 0 12
Effect of the synthetic refresh. At the collision-dominated pairs, Table 1 gives ρ ( A G ) ≤ 1.4 × 10 − 2 \rho(\mathcal{A}_{\mathcal{G}})\le1.4\times10^{-2} ρ ( A G ) ≤ 1.4 × 1 0 − 2 10 12 10^{12} 1 0 12 ∥ G [ e n ] ∥ / ∥ K [ e n ] ∥ \|\mathcal{G}[e^n]\|/\|\mathcal{K}[e^n]\| ∥ G [ e n ] ∥/∥ K [ e n ] ∥ 5.0 × 10 − 3 5.0\times10^{-3} 5.0 × 1 0 − 3 1.1 × 10 − 2 1.1\times10^{-2} 1.1 × 1 0 − 2 8.7 × 10 − 2 8.7\times10^{-2} 8.7 × 1 0 − 2 Proposition 3 predict. At the two hydrodynamic pairs ℓ G K k ≈ 1.4 \ell_{\mathrm{GK}}k\approx1.4 ℓ GK k ≈ 1.4 ρ ( A G ) \rho(\mathcal{A}_{\mathcal{G}}) ρ ( A G )
PGD error, γ G ≤ 2.8 × 10 − 3 \gamma_{G}\leq2.8\times10^{-3} γ G ≤ 2.8 × 1 0 − 3 ρ ( A G ) \rho(A_{G}) ρ ( A G ) 4 / 5 4/5 4/5 § 9 ) and PGD–GSIS needs 133 pairs, whereas the margins of source iteration and of the pure-PGD sweep are close to one. This small GSIS margin at large K n C \mathrm{Kn}_{C} Kn C § 5.1 .
Consistency of kinetic and synthetic fields. At a self-consistent solution of the continuous kinetic problem with consistent wall data, G [ I m ] = K [ I m ] G[I_{m}]=K[I_{m}] G [ I m ] = K [ I m ] (5) by r 0 = ∇ ⋅ q k i n r_{0}=\nabla\cdot q^{\mathrm{kin}} r 0 = ∇ ⋅ q kin r 1 = 1 3 ∇ T k i n + ∇ ⋅ N k i n + δ R q k i n r_{1}=\frac{1}{3}\nabla T^{\mathrm{kin}}+\nabla\cdot N^{\mathrm{kin}}+\delta_{R}q^{\mathrm{kin}} r 1 = 3 1 ∇ T kin + ∇ ⋅ N kin + δ R q kin r 2 = 2 5 ∂ ⟨ j q i ⟩ k i n + ∇ ⋅ M k i n + δ C N k i n r_{2}=\frac{2}{5}\partial_{\langle j}q_{i\rangle}^{\mathrm{kin}}+\nabla\cdot M^{\mathrm{kin}}+\delta_{C}N^{\mathrm{kin}} r 2 = 5 2 ∂ ⟨ j q i ⟩ kin + ∇ ⋅ M kin + δ C N kin N k i n N^{\mathrm{kin}} N kin G [ I m ] − K [ I m ] G[I_{m}]-K[I_{m}] G [ I m ] − K [ I m ] − r 0 -r_{0} − r 0 − r 1 + K n C ∇ ⋅ r 2 -r_{1}+\mathrm{Kn}_{C}\nabla\cdot r_{2} − r 1 + Kn C ∇ ⋅ r 2 r 0 = 0 r_{0}=0 r 0 = 0 Δ T = − 3 r 1 x / ( i k ) + 3 K n C r 2 x x \Delta T=-3r_{1x}/(ik)+3\mathrm{Kn}_{C}r_{2xx} Δ T = − 3 r 1 x / ( ik ) + 3 Kn C r 2 xx K n C \mathrm{Kn}_{C} Kn C § 4.4 quantifies bulk matching at small K n C \mathrm{Kn}_{C} Kn C § 5.2 discusses node-dependent residuals in a shared parametric basis.
Conditional asymptotic-preserving property
Appendix A first expands the frozen-source pair formally in K n C \mathrm{Kn}_{C} Kn C [30 ] did for the Shakhov model. To leading order the pair is an alternating least-squares iteration for the singular value decomposition of the local equilibrium, with a contraction factor independent of K n C \mathrm{Kn}_{C} Kn C O ( K n C ) O(\mathrm{Kn}_{C}) O ( Kn C ) η n \eta_{n} η n ζ n \zeta_{n} ζ n ω n \omega_{n} ω n Proposition 4 derives
∥ K [ I m n ] − G [ I m n ] ∥ 2 ≤ C ( K n C 2 + η n + ζ n + ω n ) . (26) \left\lVert\mathcal{K}[I_m^n]-\mathcal{G}[I_m^n]\right\rVert_2 \le C\left(\mathrm{Kn}_C^2+\eta_n+\zeta_n+\omega_n\right). \tag*{(26)} ∥ K [ I m n ] − G [ I m n ] ∥ 2 ≤ C ( Kn C 2 + η n + ζ n + ω n ) . (26) Thus η n + ζ n + ω n = O ( K n C 1 + σ A P ) \eta_n+\zeta_n+\omega_n=O(\mathrm{Kn}_C^{1+\sigma_{\mathrm{AP}}}) η n + ζ n + ω n = O ( Kn C 1 + σ AP ) 0 < σ A P ≤ 1 0<\sigma_{\mathrm{AP}}\le1 0 < σ AP ≤ 1 K n R / 3 \mathrm{Kn}_R/3 Kn R /3 δ C \delta_C δ C K n C \mathrm{Kn}_C Kn C O ( 1 ) O(1) O ( 1 )
The collision-dominated results are conditional bulk statements. The greedy spatial update enforces only its current-mode row; frozen rows and polar-space completeness must also be controlled. The amplitude tests of Algorithm 1 do not imply these bounds, and the reported fixed tolerances and caps are not claimed to certify AP for arbitrarily small Knudsen numbers. Parametric certification requires the bounds node by node. Wall layers, compatibility of the discrete DG/HDG moments, angular quadrature and a uniform iteration budget require further analysis. Preservation of the limiting constitutive law is distinct from the acceleration measured in § 4.3 .
Parametric PGD–GSIS over the Knudsen numbers
Representation and factor equations
Let δ = ( δ R , δ N ) \boldsymbol{\delta}=(\delta_R,\delta_N) δ = ( δ R , δ N ) P \mathcal{P} P d ν \mathrm{d}\nu d ν
I m ( x , θ , ϕ ; δ ) = ∑ r = 1 m X r ( x , θ ) Φ r ( ϕ , δ ) , T m k i n = ∑ r = 1 m Y r ( x ) T r ( δ ) , q m k i n = ∑ r = 1 m Z r ( x ) Q r ( δ ) . (27) I_m(x,\theta,\phi;\boldsymbol{\delta})=\sum_{r=1}^{m}X_r(x,\theta)\Phi_r(\phi,\boldsymbol{\delta}),\qquad T_m^{\mathrm{kin}}=\sum_{r=1}^{m}Y_r(x)T_r(\boldsymbol{\delta}),\qquad q_m^{\mathrm{kin}}=\sum_{r=1}^{m}Z_r(x)Q_r(\boldsymbol{\delta}). \tag*{(27)} I m ( x , θ , ϕ ; δ ) = r = 1 ∑ m X r ( x , θ ) Φ r ( ϕ , δ ) , T m kin = r = 1 ∑ m Y r ( x ) T r ( δ ) , q m kin = r = 1 ∑ m Z r ( x ) Q r ( δ ) . (27) The spatial factors are shared by all parameter values; ( ϕ , δ R , δ N ) (\phi,\delta_R,\delta_N) ( ϕ , δ R , δ N ) δ C \delta_C δ C 1 / sin ϕ 1/\sin\phi 1/ sin ϕ (13) as a product. Each inverse Knudsen number takes N N N [ δ min , δ max ] [\delta_{\min},\delta_{\max}] [ δ m i n , δ m a x ] d ν \mathrm{d}\nu d ν d ( log δ ) / log ( δ max / δ min ) \mathrm{d}(\log\delta)/\log(\delta_{\max}/\delta_{\min}) d ( log δ ) / log ( δ m a x / δ m i n ) log δ \log\delta log δ K n R , K n N ∈ [ 10 − 2 , 10 2 ] \mathrm{Kn}_R,\mathrm{Kn}_N\in[10^{-2},10^2] Kn R , Kn N ∈ [ 1 0 − 2 , 1 0 2 ] N = 21 N=21 N = 21 N δ = 441 N_\delta=441 N δ = 441 P \mathcal{P} P
With the source frozen at every node, the polar update is Eq. (16) evaluated node by node: the nodes are independent and each update is again a division. The spatial update is Eq. (17) with parameter-averaged coefficients,
a ˉ m r = ∫ P ∫ 0 π δ C Φ m Φ r d ϕ d ν , b ˉ m r = ∫ P ∫ 0 π sin ϕ Φ m Φ r d ϕ d ν , \bar{a}_{mr}=\int_{\mathcal{P}}\int_{0}^{\pi}\delta_{\mathrm{C}}\Phi_m\Phi_r\,\mathrm{d}\phi\,\mathrm{d}\nu,\qquad\bar{b}_{mr}=\int_{\mathcal{P}}\int_{0}^{\pi}\sin\phi\,\Phi_m\Phi_r\,\mathrm{d}\phi\,\mathrm{d}\nu, a ˉ m r = ∫ P ∫ 0 π δ C Φ m Φ r d ϕ d ν , b ˉ m r = ∫ P ∫ 0 π sin ϕ Φ m Φ r d ϕ d ν , a ˉ m m X m + b ˉ m m s ⋅ ∇ h X m = S T , m + s ⋅ S q , m − ∑ r = 1 m − 1 ( a ˉ m r X r + b ˉ m r s ⋅ ∇ h X r ) . (28) \bar{a}_{mm}X_m+\bar{b}_{mm}\boldsymbol{s}\cdot\nabla_hX_m=S_{T,m}+\boldsymbol{s}\cdot\boldsymbol{S}_{q,m}-\sum_{r=1}^{m-1}\left(\bar{a}_{mr}X_r+\bar{b}_{mr}\boldsymbol{s}\cdot\nabla_hX_r\right). \tag*{(28)} a ˉ mm X m + b ˉ mm s ⋅ ∇ h X m = S T , m + s ⋅ S q , m − r = 1 ∑ m − 1 ( a ˉ m r X r + b ˉ m r s ⋅ ∇ h X r ) . (28) with S T , m = 1 4 π ∫ P ∫ 0 π δ C Φ m T s d ϕ d ν S_{T,m}=\frac{1}{4\pi}\int_{\mathcal{P}}\int_{0}^{\pi}\delta_{\mathrm{C}}\Phi_mT_s\,\mathrm{d}\phi\,\mathrm{d}\nu S T , m = 4 π 1 ∫ P ∫ 0 π δ C Φ m T s d ϕ d ν S q , m = 3 4 π ∫ P ∫ 0 π δ N sin ϕ Φ m q s d ϕ d ν S_{q,m}=\frac{3}{4\pi}\int_{\mathcal{P}}\int_{0}^{\pi}\delta_{\mathrm{N}}\sin\phi\,\Phi_m\boldsymbol{q}_s\,\mathrm{d}\phi\,\mathrm{d}\nu S q , m = 4 π 3 ∫ P ∫ 0 π δ N sin ϕ Φ m q s d ϕ d ν single effective extinction a ˉ m m / b ˉ m m \bar{a}_{mm}/\bar{b}_{mm} a ˉ mm / b ˉ mm Φ m 2 \Phi_m^2 Φ m 2 δ C \delta_{\mathrm{C}} δ C sin ϕ \sin\phi sin ϕ δ C \delta_{\mathrm{C}} δ C
The separated distribution is small (§ 4.2 ), but the synthetic operator is stored per node: the local HDG factors alone take ( 7 n p ) 2 N e l (7n_p)^2N_{\mathrm{el}} ( 7 n p ) 2 N el p = 3 p=3 p = 3 N e l = 200 N_{\mathrm{el}}=200 N el = 200 N δ = 441 N_{\delta}=441 N δ = 441 § 6.5 ). They are independent of the rank and could be computed on the fly or shared between neighboring nodes; we have not pursued this.
Joint polar update. With the spatial factors fixed, all retained polar factors are re-solved against the refreshed source:
∑ r = 1 m ( δ C α k r + β k r sin ϕ ) Φ r = F T , k + sin ϕ ( F q , k − μ k ) , k = 1 , … , m , (29) \sum_{r=1}^{m}(\delta_{\mathrm{C}}\alpha_{kr}+\beta_{kr}\sin\phi)\Phi_r=F_{T,k}+\sin\phi\,(F_{q,k}-\mu_k),\qquad k=1,\ldots,m, \tag*{(29)} r = 1 ∑ m ( δ C α k r + β k r sin ϕ ) Φ r = F T , k + sin ϕ ( F q , k − μ k ) , k = 1 , … , m , (29) an m × m m\times m m × m ( ϕ , δ ) (\phi,\boldsymbol{\delta}) ( ϕ , δ ) span { X 1 , … , X m } ⊗ L 2 ( 0 , π ) \operatorname{span}\{X_1,\ldots,X_m\}\otimes L^2(0,\pi) span { X 1 , … , X m } ⊗ L 2 ( 0 , π )
The joint update is required for accuracy over the whole parameter box. Without it, Φ r ( ⋅ , δ ) \Phi_r(\cdot,\boldsymbol{\delta}) Φ r ( ⋅ , δ ) r r r U s ( δ ) U_s(\boldsymbol{\delta}) U s ( δ ) Proposition 3 ), but at large K n C \mathrm{Kn}_{\mathrm{C}} Kn C O ( K n C ) O(\mathrm{Kn}_{\mathrm{C}}) O ( Kn C ) § 4.3 ) and the contraction margin of full-order GSIS shrinks to 1 − ρ ( A G ) → 1 / 5 1-\rho(A_G)\to1/5 1 − ρ ( A G ) → 1/5 Table 1 ). The perturbation γ C n \gamma_{\mathrm{C}}^n γ C n (24) can then exceed that margin: the construction loses accuracy at large Knudsen numbers and can diverge at the rarefied corner (§ 6.5 ). With the joint update, the nodal state is the Galerkin solution on the current spatial basis, so the representation error is governed by that basis rather than by stale polar factors. It is not guaranteed to decrease monotonically, because a Galerkin projection of a transport operator is not a best approximation.
Other remedies are possible. The parameter domain can be divided into subregions, for example by ranges of K n C \mathrm{Kn}_{C} Kn C α k r \alpha_{kr} α k r β k r \beta_{kr} β k r δ \delta δ m m m O ( m N x N δ ) \mathcal{O}(mN_xN_\delta) O ( m N x N δ ) m × m m \times m m × m O ( m 3 N ϕ N δ ) \mathcal{O}(m^3N_\phi N_\delta) O ( m 3 N ϕ N δ ) N δ N_\delta N δ § 6.5 ).
Parametric pure PGD. For comparison, a parametric pure PGD can use the joint update in two forms. With a frozen kinetic source, U s ← K [ I m ] U_s \leftarrow\mathcal{K}[I_m] U s ← K [ I m ] (29) applies unchanged, and Eq. (24) with K \mathcal{K} K G \mathcal{G} G Table 1 ). With the implicit coupling of § 4.1 , the joint update reads, at every parameter node,
∑ r = 1 m [ ( δ C α k r + β k r sin ϕ ) Φ r − δ C α k r T M 1 ( Φ r ) − δ N α k r q sin ϕ M 2 ( Φ r ) ] = − μ k sin ϕ , k = 1 , … , m . (30) \sum_{r=1}^{m}\left[\left(\delta_C\alpha_{kr}+\beta_{kr}\sin\phi\right)\Phi_r-\delta_C\alpha_{kr}^{T}M_1(\Phi_r)-\delta_N\alpha_{kr}^{q}\sin\phi M_2(\Phi_r)\right]=-\mu_k\sin\phi,\qquad k=1,\ldots,m. \tag*{(30)} r = 1 ∑ m [ ( δ C α k r + β k r sin ϕ ) Φ r − δ C α k r T M 1 ( Φ r ) − δ N α k r q sin ϕ M 2 ( Φ r ) ] = − μ k sin ϕ , k = 1 , … , m . (30) It couples all polar nodes through M 1 M_1 M 1 M 2 M_2 M 2 m N ϕ × m N ϕ mN_\phi\times mN_\phi m N ϕ × m N ϕ ϕ \phi ϕ 2 m 2m 2 m § 4.3 , which move macroscopic information by about one mean free path per sweep. In our tests, parametric pure PGD with the frozen kinetic source did not reach the correct solution with 30 modes; with the implicit joint update it did, but only with considerably more modes and alternating pairs than PGD–GSIS, whose macroscopic source comes from the synthetic equations and whose joint update needs only m × m m \times m m × m
Parametric rank and error sensitivity
Equation [27 ] places every T m k i n ( ⋅ , δ ) T_m^{\mathrm{kin}}(\cdot,\delta) T m kin ( ⋅ , δ ) span { Y 1 , … , Y m } \operatorname{span}\{Y_1,\ldots,Y_m\} span { Y 1 , … , Y m } M T = { T ( ⋅ ; δ ) } δ ∈ P \mathcal{M}_T=\{T(\cdot;\delta)\}_{\delta\in\mathcal{P}} M T = { T ( ⋅ ; δ ) } δ ∈ P ε T \varepsilon_T ε T d m ( M T ) ≤ ε T d_m(\mathcal{M}_T)\leq\varepsilon_T d m ( M T ) ≤ ε T d m d_m d m m m m
For Fourier states with ∥ q ∥ ≲ K n R \lVert q\rVert\lesssim\mathrm{Kn}_R ∥ q ∥ ≲ Kn R e = I m − I e=I_m-I e = I m − I ∥ e ∥ / K n R \lVert e\rVert/\mathrm{Kn}_R ∥ e ∥ / Kn R
Algorithm 2 Parametric PGD–GSIS
Input: parameter nodes and weights, mesh and quadrature, tolerances and caps
Factorize the local and global HDG systems and wall Jacobians at every parameter node.
for m = 1 , … , M max m = 1, \ldots, M_{\max} m = 1 , … , M m a x do
Initialize X m X_m X m
repeat
Polar update (16) at every node, source frozen.
Spatial update (28), source frozen; normalize the pair.
Synthetic refresh U s ( δ ) ← G δ [ I m ] U_s(\delta) \leftarrow\mathcal{G}_{\delta}[I_m] U s ( δ ) ← G δ [ I m ]
Joint polar update (29) of all retained polar factors at every node.
until the amplitude test holds or L max L_{\max} L m a x
Retain ( X m , Φ m ) (X_m,\Phi_m) ( X m , Φ m ) r e n ( m ) < ϵ e n r_{\mathrm{en}}^{(m)} < \epsilon_{\mathrm{en}} r en ( m ) < ϵ en
end for
Output: shared spatial factors X r X_r X r Φ r ( ⋅ , δ j ) \Phi_r(\cdot,\delta_j) Φ r ( ⋅ , δ j )
errors at individual nodes. This motivates nodewise error checks and, if needed, a parameter measure that emphasizes the small-K n R \mathrm{Kn}_R Kn R Proposition 4 must likewise hold node by node; a small parameter-weighted residual alone does not establish them.
At large K n C \mathrm{Kn}_C Kn C O ( K n C ) \mathcal{O}(\mathrm{Kn}_C) O ( Kn C ) § 4.3 ), making the ballistic corner another sensitive region. The joint update of § 5.1 removes the stale-factor part of these residuals; with it, all our constructions remained stable over the whole parameter box (§ 6.5 ).
Nodal macroscopic reconstruction costs O ( m N x ) \mathcal{O}(mN_x) O ( m N x ) T r T_r T r Q r Q_r Q r § 6.5 assesses it independently by withholding nodes.
Numerical examples
Numerical setup and error measures
The fixed-parameter tests use six pairs ( K n R , K n N ) (\mathrm{Kn}_R,\mathrm{Kn}_N) ( Kn R , Kn N ) ( 0.01 , 0.01 ) (0.01, 0.01) ( 0.01 , 0.01 ) ( 0.01 , 100 ) (0.01, 100) ( 0.01 , 100 ) ( 0.1 , 10 ) (0.1, 10) ( 0.1 , 10 ) ( 10 , 0.1 ) (10, 0.1) ( 10 , 0.1 ) ( 100 , 0.01 ) (100, 0.01) ( 100 , 0.01 ) ( 100 , 100 ) (100, 100) ( 100 , 100 ) K n C ≤ 0.1 \mathrm{Kn}_C \le0.1 Kn C ≤ 0.1 K n C = 50 \mathrm{Kn}_C = 50 Kn C = 50 (7) gives ℓ G K ≈ 0.44 \ell_{\mathrm{GK}} \approx0.44 ℓ GK ≈ 0.44 0.45 0.45 0.45 H = 0.2 H = 0.2 H = 0.2
The spatial factor uses polynomials of degree p = 3 p = 3 p = 3 N θ = 80 N_\theta= 80 N θ = 80 [ 0 , π ] [0,\pi] [ 0 , π ] [ π , 2 π ] [\pi,2\pi] [ π , 2 π ] N ϕ = 40 N_{\phi}=40 N ϕ = 40 [ 0 , π ] [0,\pi] [ 0 , π ] N s w = 5 N_{\mathrm{sw}}=5 N sw = 5 L max = 10 L_{\max}=10 L m a x = 10 M max = 20 M_{\max}=20 M m a x = 20 ϵ s w = ϵ A L S = 10 − 3 \epsilon_{\mathrm{sw}}=\epsilon_{\mathrm{ALS}}=10^{-3} ϵ sw = ϵ ALS = 1 0 − 3 ϵ e n = 10 − 3 \epsilon_{\mathrm{en}}=10^{-3} ϵ en = 1 0 − 3
E T = ∥ T m k i n − T F R ∥ h max ( ∥ T F R ∥ h , s T ) , E q = ∥ q m k i n − q F R ∥ h max ( ∥ q F R ∥ h , s q ) , (31) E_T=\frac{\lVert T_m^{\mathrm{kin}}-T_{\mathrm{FR}}\rVert_h}{\max(\lVert T_{\mathrm{FR}}\rVert_h,s_T)},\qquad E_q=\frac{\lVert q_m^{\mathrm{kin}}-q_{\mathrm{FR}}\rVert_h}{\max(\lVert q_{\mathrm{FR}}\rVert_h,s_q)}, \tag*{(31)} E T = max (∥ T FR ∥ h , s T ) ∥ T m kin − T FR ∥ h , E q = max (∥ q FR ∥ h , s q ) ∥ q m kin − q FR ∥ h , (31) with small positive floors s T , s q s_T,s_q s T , s q ∣ T − T F R ∣ / max ∣ T F R ∣ \lvert T-T_{\mathrm{FR}}\rvert/\max\lvert T_{\mathrm{FR}}\rvert ∣ T − T FR ∣ / max ∣ T FR ∣ E q E_q E q
PGD and PGD–GSIS across transport regimes
Figure 2 compares pure PGD and PGD–GSIS on the cavity. Pure PGD has maximum pointwise errors of 1.7 × 10 − 1 1.7\times10^{-1} 1.7 × 1 0 − 1 ( 0.01 , 0.01 ) (0.01,0.01) ( 0.01 , 0.01 ) 2.8 × 10 − 2 2.8\times10^{-2} 2.8 × 1 0 − 2 ( 0.01 , 100 ) (0.01,100) ( 0.01 , 100 ) 4.0 4.0 4.0 ( 100 , 0.01 ) (100,0.01) ( 100 , 0.01 ) [ 0 , 1 ] [0,1] [ 0 , 1 ] 6.0 × 10 − 6 6.0\times10^{-6} 6.0 × 1 0 − 6 1 − ρ s w ≤ 3 × 10 − 4 1-\rho_{\mathrm{sw}}\leq3\times10^{-4} 1 − ρ sw ≤ 3 × 1 0 − 4 Table 1 ), so that even the full budget of 50 sweeps per mode removes less than 2% of the slow error components; at the ballistic pair the margin is close to one. The hydrodynamic failure is the most severe. There the slow component is the heat flux, which normal scattering nearly conserves, and enrichment increases the error from 0.63 at rank one to 3.2 at rank twenty (Figure 5 ): new modes are fitted to an unconverged state instead of correcting it. PGD–GSIS reduces the maximum errors to 1.1 × 10 − 5 1.1\times10^{-5} 1.1 × 1 0 − 5 2.2 × 10 − 5 2.2\times10^{-5} 2.2 × 1 0 − 5 2.4 × 10 − 4 2.4\times10^{-4} 2.4 × 1 0 − 4 5.2 × 10 − 6 5.2\times10^{-6} 5.2 × 1 0 − 6 0.24 ≤ T ≤ 0.26 0.24\leq T\leq0.26 0.24 ≤ T ≤ 0.26 [ 0.08 , 0.49 ] [0.08,0.49] [ 0.08 , 0.49 ]
Figure 5. Relative L 2 L^2 L 2 T T T T T T q x q_x q x
The channel problems show the same dichotomy (Figure 3 and Figure 4 ), and PGD–GSIS lies on the full-order profiles in every regime. The profiles follow the mechanisms of § 2.4 . In cross-plane conduction the diffusive profile is the Fourier line T = x T=x T = x ( 100 , 0.01 ) (100,0.01) ( 100 , 0.01 ) ( 100 , 100 ) (100,100) ( 100 , 100 ) T = 0.5 T=0.5 T = 0.5 ∂ x T = − 3 δ R q x \partial_xT=-3\delta_Rq_x ∂ x T = − 3 δ R q x K n R / 3 = 3.3 × 10 − 3 \mathrm{Kn}_R/3=3.3\times10^{-3} Kn R /3 = 3.3 × 1 0 − 3 ℓ G K \ell_{\mathrm{GK}} ℓ GK K n C \mathrm{Kn}_C Kn C ( 100 , 0.01 ) (100,0.01) ( 100 , 0.01 ) ℓ G K = 0.45 > H \ell_{\mathrm{GK}}=0.45>H ℓ GK = 0.45 > H 5 H 2 / ( 24 K n C ) = 0.83 5H^2/(24\mathrm{Kn}_C)=0.83 5 H 2 / ( 24 Kn C ) = 0.83 ( 100 , 100 ) (100,100) ( 100 , 100 ) ( 0.01 , 0.01 ) (0.01,0.01) ( 0.01 , 0.01 )
Figure 3. Cross-plane conduction: T T T y = H / 2 y = H/2 y = H /2 T = x T = x T = x
Figure 4. Phonon Poiseuille flow: q x q_x q x Figure 3 .
Rank convergence and modal contributions
Figure 5 shows the error against the number of retained modes. With PGD–GSIS the diffusive pairs reach about 10 − 4 10^{-4} 1 0 − 4 7 × 10 − 6 7\times10^{-6} 7 × 1 0 − 6
Figure 2. Isothermal cavity at ( K n R , K n N ) = ( 0.01 , 0.01 ) (\mathrm{Kn}_{R},\mathrm{Kn}_{N})=(0.01,0.01) ( Kn R , Kn N ) = ( 0.01 , 0.01 ) ( 0.01 , 100 ) (0.01,100) ( 0.01 , 100 ) ( 100 , 0.01 ) (100,0.01) ( 100 , 0.01 ) ( 100 , 100 ) (100,100) ( 100 , 100 ) ∣ T − T F R ∣ / max ∣ T F R ∣ \lvert T-T_{\mathrm{FR}}\rvert/\max\lvert T_{\mathrm{FR}}\rvert ∣ T − T FR ∣ / max ∣ T FR ∣
end between 10 − 5 10^{-5} 1 0 − 5 10 − 4 10^{-4} 1 0 − 4 3.1 × 10 − 6 3.1 \times10^{-6} 3.1 × 1 0 − 6 7.9 × 10 − 4 7.9 \times10^{-4} 7.9 × 1 0 − 4 8.6 × 10 − 5 8.6 \times10^{-5} 8.6 × 1 0 − 5
Figure 6 shows the first six PGD–GSIS macroscopic contributions for the channel problems. In cross-plane conduction the first temperature contribution carries the nearly linear diffusive profile and the plateaus of the hydrodynamic and ballistic pairs. The sum of modes two to six, relative to the rank-six temperature field in L 2 L^2 L 2 2.5 × 10 − 4 2.5 \times10^{-4} 2.5 × 1 0 − 4 5.4 × 10 − 2 5.4 \times10^{-2} 5.4 × 1 0 − 2 3.0 × 10 − 2 3.0 \times10^{-2} 3.0 × 1 0 − 2 (15) concerns the full distribution. The ranks at which enrichment stops (Table 2 ) follow the ordering of § 3.3 : three in the diffusive limit, four to seven in the ballistic limit, and most, up to fourteen, at (10, 0.1), where optical depths of order one to ten span the domain and wall layers are thick.
Figure 6. First six PGD–GSIS macroscopic contributions in the channel problems. Top (a): cross-plane temperature contributions T r ( x ) = Y r ( x ) T r T_r(x) = Y_r(x)T_r T r ( x ) = Y r ( x ) T r q x , r ( y ) q_{x,r}(y) q x , r ( y ) ( K n R , K n N ) = ( 0.01 , 100 ) (\mathrm{Kn}_R,\mathrm{Kn}_N) = (0.01, 100) ( Kn R , Kn N ) = ( 0.01 , 100 ) ( 100 , 0.01 ) (100, 0.01) ( 100 , 0.01 ) ( 100 , 100 ) (100, 100) ( 100 , 100 ) r = 1 , … , 6 r = 1,\ldots,6 r = 1 , … , 6
Table 2. Rank at which enrichment stops, with the average number of alternating pairs per mode in parentheses; †: rank cap M max = 20 M_{\max}=20 M m a x = 20 ϵ e n \epsilon_{\mathrm{en}} ϵ en § 4.1 ).
The modal amplitudes in Figure 7 record the actual evolution of successive macroscopic contributions during enrichment and complement the reference-error curves of Figure 5 . They are the quantities monitored by the increment test (23); these diagnostic runs continue to twenty modes. For PGD–GSIS every curve falls by more than two decades over twenty modes; for pure PGD the hydrodynamic curve grows or stagnates.
Figure 7. Modal amplitudes during enrichment for pure PGD (top) and PGD–GSIS (bottom), with enrichment continued to twenty modes. Columns as in Figure 5 . The curves show the evolution of successive macroscopic contributions; only the decay along a curve is meaningful, and the amplitudes are not error bounds.
amplitude is nevertheless not an error estimate: the ballistic Poiseuille pair decays more slowly than several others yet reaches a small error, because the separated terms are not orthogonal.
Computational cost and memory scaling
Table 2 and Figure 8 compare four solvers on all pairs with 12 OpenMP threads: full-order source iteration (CIS), full-order GSIS, pure PGD, and PGD–GSIS with enrichment stopped at ϵ e n = 10 − 3 \epsilon_{\mathrm{en}} = 10^{-3} ϵ en = 1 0 − 3
Figure 8. Wall-clock time (top, logarithmic) and peak resident memory (bottom) of the four solvers at the six Knudsen pairs; N θ × N ϕ = 80 × 40 N_{\theta} \times N_{\phi}=80 \times40 N θ × N ϕ = 80 × 40
Full-order CIS needs 10 2 10^{2} 1 0 2 10 4 10^{4} 1 0 4 1 − ρ ( A K ) 1-\rho(A_{K}) 1 − ρ ( A K ) Table 1 is smallest. Full-order GSIS converges in 8–216 iterations and 1.1–17.7 s. PGD–GSIS stops at ranks 3–14 after 15–114 synthetic updates in total—the same order as full-order GSIS, as the analysis of § 4.3 predicts—and takes 1.2–40.4 s. Its heat-flux error stays below 6 × 10 − 3 6 \times10^{-3} 6 × 1 0 − 3 7 × 10 − 4 7 \times10^{-4} 7 × 1 0 − 4 Figure 9 ). At the diffusive pairs E q E_q E q E T E_T E T E q E_q E q ∥ e ∥ / K n R \lVert e\rVert/\mathrm{Kn}_{R} ∥ e ∥ / Kn R § 5.2 . Pure PGD takes 0.5–29.8 s, but at the stalled pairs this is the cost of exhausting the budget of 200 alternating pairs without converging (E q E_q E q 1 − ρ s w 1-\rho_{\mathrm{sw}} 1 − ρ sw
Figure 9. Errors E T E_T E T E q E_q E q E T E_T E T
At N ϕ = 40 N_{\phi}=40 N ϕ = 40 § 4.2 ),
and neither solver is limited by the number of updates. Its peak memory (58–131 MiB) is below that of full-order GSIS (93–162 MiB) but above that of full-order CIS (roughly 45–80 MiB), which stores no synthetic operator: at this angular resolution the peak is dominated by arrays that do not depend on N ϕ N_\phi N ϕ Figure 10 ). Full-order memory increases by about 1.2 MiB per polar node, the size of one polar slice of the distribution (N x N θ N_xN_\theta N x N θ ϕ \phi ϕ N ϕ = 160 N_\phi= 160 N ϕ = 160 2.7 × 10 2 2.7 \times10^2 2.7 × 1 0 2 2.0 × 10 3 2.0 \times10^3 2.0 × 1 0 3 N ϕ N_\phi N ϕ N ϕ N_\phi N ϕ § 4.2 removes that dependence.
Figure 10. Cavity: peak resident memory (top) and wall-clock time (bottom) versus the polar node count at N θ = 80 N_\theta= 80 N θ = 80
Accuracy and cost of parametric PGD–GSIS
Each test problem carries one parametric construction over K n R , K n N ∈ [ 10 − 2 , 10 2 ] \mathrm{Kn}_R,\mathrm{Kn}_N \in[10^{-2}, 10^2] Kn R , Kn N ∈ [ 1 0 − 2 , 1 0 2 ] 10 − 3 10^{-3} 1 0 − 3 Figure 11 pairs the first three temperature spatial factors of the cavity construction with their parameter-dependent coefficients. To expose shape and magnitude separately, set Y r ∗ = Y r / max ∣ Y r ∣ Y_r^* = Y_r/\max|Y_r| Y r ∗ = Y r / max ∣ Y r ∣ T r ∗ = max ∣ Y r ∣ T r T_r^* = \max|Y_r|T_r T r ∗ = max ∣ Y r ∣ T r Y r ∗ T r ∗ = Y r T r Y_r^*T_r^* = Y_rT_r Y r ∗ T r ∗ = Y r T r ∣ T r ∗ ( δ ) ∣ |T_r^*(\delta)| ∣ T r ∗ ( δ ) ∣ r r r
Figure 11. First three temperature factor–coefficient pairs of the parametric cavity construction, with r = 1 , 2 , 3 r = 1, 2, 3 r = 1 , 2 , 3 Y r ∗ Y_r^{*} Y r ∗ [ − 1 , 1 ] [-1, 1] [ − 1 , 1 ] T r ∗ T_r^{*} T r ∗ ∣ T r ∗ ∣ |T_r^{*}| ∣ T r ∗ ∣
The first spatial factor is one-signed and peaks at the hot wall; its coefficient varies from about 0.70 to 0.81 across the parameter box. On the Fourier side the first two temperature contributions dominate, at about 0.70 and 0.18. The second coefficient changes sign toward the kinetic side and reaches about −0.5 at K n R = 10 2 \mathrm{Kn}_R = 10^2 Kn R = 1 0 2
The reconstruction is scored against 225 independent full-order GSIS solutions on a 15 × 15 15 \times15 15 × 15 Figure 12 , Table 3 ); both channel constructions are more accurate. The error distribution follows § 5.2 . The largest E q E_q E q K n R = 10 − 2 \mathrm{Kn}_R = 10^{-2} Kn R = 1 0 − 2 O ( K n R ) O(\mathrm{Kn}_R) O ( Kn R ) E T E_T E T 10 − 2 10^{-2} 1 0 − 2
Figure 12. Parametric cavity reconstruction at rank 30: relative error of T T T ∣ q ∣ |q| ∣ q ∣ ℓ 2 \ell^2 ℓ 2 101 × 101 101 \times101 101 × 101 1.4 × 10 − 3 1.4 \times10^{-3} 1.4 × 1 0 − 3 1.1 × 10 − 2 1.1 \times10^{-2} 1.1 × 1 0 − 2 1.9 × 10 − 3 1.9 \times10^{-3} 1.9 × 1 0 − 3 3.9 × 10 − 2 3.9 \times10^{-2} 3.9 × 1 0 − 2
Table 3. Parametric constructions (21 × 21 nodes): smallest/largest relative error (%) over the 225 full-order references at three truncation ranks of the same stored factors. E q E_q E q ∣ q ∣ |q| ∣ q ∣ q x q_x q x E T E_T E T
These errors are one to two orders of magnitude above the fixed-parameter errors at the same pairs, and Table 3 shows why they should not be read as the limit of the method: from rank 10 to 30 the largest cavity errors fall by factors of 4.0 (E T E_T E T E q E_q E q E q E_q E q
of the thirtieth mode is still 0.68% (cavity), 0.15% (cross-plane) and 0.04% (Poiseuille) of the first. Every construction stopped at the rank cap. The larger ranks are consistent with the shared-basis requirement of § 5.2 : one space must approximate the temperature and flux fields across the whole parameter box.
Interpolation between nodes was assessed by withholding every other node in each direction: the coefficients T r T_r T r Q r Q_r Q r ( log 10 K n R , log 10 K n N ) (\log_{10} \mathrm{Kn}_R,\log_{10} \mathrm{Kn}_N) ( log 10 Kn R , log 10 Kn N ) 2.7 × 10 − 3 2.7 \times10^{-3} 2.7 × 1 0 − 3 4.9 × 10 − 3 4.9 \times10^{-3} 4.9 × 1 0 − 3 1.0 × 10 − 2 1.0 \times10^{-2} 1.0 × 1 0 − 2 6.9 × 10 − 2 6.9 \times10^{-2} 6.9 × 1 0 − 2
The three constructions take 1.26 min (cavity), 0.52 min (cross-plane) and 0.60 min (Poiseuille) on 32 threads, with 9–10 alternating pairs per mode. Peak memory is 12.0 GiB for the cavity and 4.9 GiB for the channels, dominated by the per-node HDG operators of § 5.1 .
Role of the joint polar update. Without the joint update, the synthetic refresh is unstable at the rarefied corner. The 21 × 21 construction then diverged at ( K n R , K n N ) = ( 63 , 100 ) (\mathrm{Kn}_R,\mathrm{Kn}_N) = (63, 100) ( Kn R , Kn N ) = ( 63 , 100 ) ( 100 , 63 ) (100, 63) ( 100 , 63 ) ( 100 , 100 ) (100, 100) ( 100 , 100 ) [ 0 , 1 ] [0, 1] [ 0 , 1 ] E T = 0.93 E_T = 0.93 E T = 0.93 ( 100 , 100 ) (100, 100) ( 100 , 100 ) E T = 4.0 × 10 − 3 E_T = 4.0 \times10^{-3} E T = 4.0 × 1 0 − 3 E q = 5.3 × 10 − 3 E_q = 5.3 \times10^{-3} E q = 5.3 × 1 0 − 3 E T = 9.0 × 10 − 3 E_T = 9.0 \times10^{-3} E T = 9.0 × 1 0 − 3 E q = 2.4 × 10 − 2 E_q = 2.4 \times10^{-2} E q = 2.4 × 1 0 − 2 E T E_T E T ( 100 , 100 ) (100, 100) ( 100 , 100 ) m × m m \times m m × m K n R = 10 \mathrm{Kn}_R = 10 Kn R = 10 K n N ∈ [ 10 − 2 , 10 2 ] \mathrm{Kn}_N \in[10^{-2},10^{2}] Kn N ∈ [ 1 0 − 2 , 1 0 2 ]
Conclusions
We have combined proper generalized decomposition with the general synthetic iterative scheme for the steady linearized gray Callaway model on two-dimensional domains. Three findings organize the results.
First, the splitting of the distribution into a space–azimuth factor and a polar factor is intrinsically low-rank. Upwinding depends on the azimuth alone, and along each characteristic the polar angle enters through the stretching of the in-plane optical depth. The kernel e − τ / sin φ e^{-\tau/\sin\varphi} e − τ / s i n φ N N N N + 1 N+1 N + 1
Second, the treatment of the macroscopic coupling decides convergence. Pure PGD keeps the scattering source implicit in the polar update, but its spatial update moves macroscopic information only through sweeps that lag the moments of the current mode. A Fourier analysis shows that these sweeps have the vanishing contraction margin of source iteration, and the reported pure-PGD experiments stall or grow at the collision-dominated pairs. PGD–GSIS freezes the source within each alternating pair and refreshes it by the synthetic equations; its error propagation splits into the full-order GSIS operator and the image of the representation error, which in the periodic bulk is read only through degree-three angular moments, with a gain that is small in the Fourier regime. The conditional AP analysis shows that smooth states with controlled projected residuals, polar approximation errors and source increments reconstruct the first streaming correction, suppress those moments and match the synthetic fields to o ( K n C ) o(\mathrm{Kn}_{C}) o ( Kn C ) 6 × 10 − 3 6 \times10^{-3} 6 × 1 0 − 3
Third, a single spatial basis can represent four decades of both Knudsen numbers, provided all polar factors are updated jointly. Without this update, polar factors fitted to an earlier source leave representation errors that the synthetic map amplifies at large K n C \mathrm{Kn}_{C} Kn C m × m m \times m m × m K n R \mathrm{Kn}_{R} Kn R O ( K n R ) \mathcal{O}(\mathrm{Kn}_{R}) O ( Kn R ) O ( K n C ) \mathcal{O}(\mathrm{Kn}_{C}) O ( Kn C )
These observations suggest the next steps. A representation of the deviation from the Callaway equilibrium, as in micro–macro low-rank schemes [36 ] , could reduce relative-flux sensitivity. Reporting the synthetic flux together with its matching defect is another option in the Fourier regime. The per-node synthetic operators, which dominate the parametric memory, could be assembled on the fly. A further direction is to introduce separate PGD factors for the two Knudsen numbers and extend this separation to the macroscopic HDG equations. Shared spatial systems could then replace per-node factorizations and solves, potentially yielding substantial reductions in memory and computational cost. Uniform attainment of the AP residual conditions with fixed rank and iteration budgets, and a fully discrete AP analysis including kinetic layers, remain open. Nongray scattering, where frequency is a natural additional coordinate, and three-dimensional geometries, where the upwind direction depends on the polar angle, are the main extensions.
CRediT authorship contribution statement
Jia Liu: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Visualization, Writing – original draft.
Wei Su: Conceptualization, Methodology, Supervision, Funding acquisition, Writing – review & editing.
Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Data availability
The code and data that support the findings of this study are available from the corresponding author upon reasonable request.
Declaration of generative AI and AI-assisted technologies in the manuscript preparation process
During the preparation of this work the authors used Claude (Anthropic) to edit the language and structure of the manuscript and to check derivations. After using this tool, the authors reviewed and edited the content as needed and take full responsibility for the content of the published article.
Asymptotic limits of PGD–GSIS
We analyze the frozen-source construction of Algorithm 1 in three limits. In the diffusive and hydrodynamic limits the collision term is stiff, and the question is asymptotic preservation. In the ballistic limit the question is whether the synthetic coupling degrades free streaming. As in the continuous-coordinate PGD analysis of [30 ] , spatial and angular discretization errors are excluded. Sections Appendix § 8 .§ 1 –Appendix § 8 .§ 3 develop the estimates common to the two collision-dominated limits, Sections Appendix § 8 .§ 4 –Appendix § 8 .§ 6 treat the three limits, and Section Appendix § 8 .§ 7 states the scope. The collision-dominated results concern smooth bulk fields on a fixed subdomain D 0 ⋐ D D_0 \Subset D D 0 ⋐ D
Formal expansion of the frozen-source pair
As in [30 ] , the mechanism can first be seen from a formal expansion of the iteration. Set λ = K n C \lambda= \mathrm{Kn}_C λ = Kn C
D f = Ω ⋅ ∇ f , E n = T s n + 3 χ Ω ⋅ q s n 4 π . (32) D f = \boldsymbol{\Omega} \cdot\nabla f,\qquad E_n = \frac{T_s^n + 3\chi\boldsymbol{\Omega} \cdot\boldsymbol{q}_s^n}{4\pi}. \tag*{(32)} D f = Ω ⋅ ∇ f , E n = 4 π T s n + 3 χ Ω ⋅ q s n . (32) Away from walls, Eqs. (16) and (17) of PGD–GSIS, divided by δ C \delta_C δ C I + λ D I = E n I + \lambda D I = E_n I + λ D I = E n R = I m + λ D I m − E n R = I_m + \lambda D I_m - E_n R = I m + λ D I m − E n
( X m , R ) x θ = 0 for every ϕ , ∫ 0 π Φ m R d ϕ = 0 for every ( x , θ ) . (33) (X_m,R)_{x\theta}=0 \quad\text{for every }\phi,\qquad\int_0^\pi\Phi_m R\,\mathrm{d}\phi=0 \quad\text{for every }(x,\theta). \tag*{(33)} ( X m , R ) x θ = 0 for every ϕ , ∫ 0 π Φ m R d ϕ = 0 for every ( x , θ ) . (33) At λ = 0 \lambda= 0 λ = 0 min ∥ E n − I m − 1 − X m Φ m ∥ \min\lVert E_n-I_{m-1}-X_m\Phi_m\rVert min ∥ E n − I m − 1 − X m Φ m ∥ Φ m \Phi_m Φ m X m X_m X m E n − I m − 1 E_n-I_{m-1} E n − I m − 1 E n − I m − 1 E_n-I_{m-1} E n − I m − 1 λ \lambda λ O ( λ 2 ) O(\lambda^2) O ( λ 2 ) § 4.3 ). By Proposition 2, E n E_n E n G 1 = span { 1 , sin ϕ } G_1=\operatorname{span}\{1,\sin\phi\} G 1 = span { 1 , sin ϕ }
Next expand I m = I ( 0 ) + λ I ( 1 ) + O ( λ 2 ) I_m=I^{(0)}+\lambda I^{(1)}+O(\lambda^2) I m = I ( 0 ) + λ I ( 1 ) + O ( λ 2 ) R R R I ( 0 ) = E n I^{(0)}=E_n I ( 0 ) = E n I ( 1 ) = − D E n I^{(1)}=-D E_n I ( 1 ) = − D E n
D E n = sin ϕ s ⋅ ∇ T s n + 3 χ sin 2 ϕ s ⋅ ∇ ( s ⋅ q s n ) 4 π . (34) D E_n=\frac{\sin\phi\,\boldsymbol{s}\cdot\nabla T_s^n+3\chi\sin^2\phi\,\boldsymbol{s}\cdot\nabla(\boldsymbol{s}\cdot\boldsymbol{q}_s^n)}{4\pi}. \tag*{(34)} D E n = 4 π sin ϕ s ⋅ ∇ T s n + 3 χ sin 2 ϕ s ⋅ ∇ ( s ⋅ q s n ) . (34) The state E n − λ D E n E_n-\lambda D E_n E n − λ D E n Ω \boldsymbol{\Omega} Ω G 2 = span { 1 , sin ϕ , sin 2 ϕ } G_2=\operatorname{span}\{1,\sin\phi,\sin^2\phi\} G 2 = span { 1 , sin ϕ , sin 2 ϕ } K i j k K_{ijk} K ij k O ( λ 2 ) O(\lambda^2) O ( λ 2 ) § 8 .§ 4 and Appendix § 8 .§ 5 .
This argument assumes three things that the algorithm does not guarantee. First, Eq. (A.2) enforces only the rows of mode m m m R R R
The first assumption matters even at a frozen source. We ran the iteration of Appendix § 9 with U s U_s U s k = π k = \pi k = π K n R = K n N = ε \mathrm{Kn}_{R} = \mathrm{Kn}_{N} = \varepsilon Kn R = Kn N = ε λ = ε / 2 \lambda= \varepsilon/2 λ = ε /2 10 − 11 10^{-11} 1 0 − 11 ε = 10 − 2 \varepsilon= 10^{-2} ε = 1 0 − 2 10 − 3 10^{-3} 1 0 − 3 3.4 × 10 − 2 3.4 \times10^{-2} 3.4 × 1 0 − 2 3.5 × 10 − 2 3.5 \times10^{-2} 3.5 × 1 0 − 2 ( σ 2 / σ 1 ) 2 = 3.5 × 10 − 2 (\sigma_2/\sigma_1)^2 = 3.5 \times10^{-2} ( σ 2 / σ 1 ) 2 = 3.5 × 1 0 − 2 L 2 L^2 L 2 7.9 × 10 − 4 7.9 \times10^{-4} 7.9 × 1 0 − 4 8.0 × 10 − 5 8.0 \times10^{-5} 8.0 × 1 0 − 5 λ \lambda λ G 2 G_2 G 2 ζ n ≤ 7 × 10 − 7 \zeta_n \leq7 \times10^{-7} ζ n ≤ 7 × 1 0 − 7 L 2 L^2 L 2 λ \lambda λ 1.4 × 10 − 4 1.4 \times10^{-4} 1.4 × 1 0 − 4 1.4 × 10 − 5 1.4 \times10^{-5} 1.4 × 1 0 − 5 η n \eta_n η n O ( λ ) O(\lambda) O ( λ ) 4.0 × 10 − 6 4.0 \times10^{-6} 4.0 × 1 0 − 6 4.1 × 10 − 8 4.1 \times10^{-8} 4.1 × 1 0 − 8 λ 2 \lambda^2 λ 2 ( K n R = 1 , K n N = ε ) (\mathrm{Kn}_{R} = 1, \mathrm{Kn}_{N} = \varepsilon) ( Kn R = 1 , Kn N = ε ) η n \eta_n η n ζ n \zeta_n ζ n § 8 .§ 2 –Appendix § 8 .§ 5 replace the three assumptions by measurable indicators and prove the corresponding estimates.
Retained-mode indicators
With λ \lambda λ D D D E n E_n E n I m n I_m^n I m n U s n U_s^n U s n I − n I_-^n I − n
U s n = G [ I − n ] , U s n + 1 = G [ I m n ] . (35) U_s^n = \mathcal{G}[I_-^n], \qquad U_s^{n+1} = \mathcal{G}[I_m^n]. \tag*{(35)} U s n = G [ I − n ] , U s n + 1 = G [ I m n ] . (35) The analysis applies after a synthetic refresh, with the same temperature gauge. Let V n = span { Φ 1 , … , Φ m } ⊂ L 2 ( 0 , π ; d ϕ ) V_n = \operatorname{span}\{\Phi_1,\ldots,\Phi_m\} \subset L^2(0,\pi;d\phi) V n = span { Φ 1 , … , Φ m } ⊂ L 2 ( 0 , π ; d ϕ ) P n P_n P n ( x , θ ) (x,\theta) ( x , θ ) ∥ ⋅ ∥ s \lVert\cdot\rVert_s ∥ ⋅ ∥ s H s ( D 0 ; L 2 ( d θ d ϕ ) ) H^s(D_0;L^2(d\theta\,d\phi)) H s ( D 0 ; L 2 ( d θ d ϕ )) H s ( D 0 ) H^s(D_0) H s ( D 0 ) x x x P n P_n P n D D D
The collision-scaled residual and its projection are
R n = I m n + λ D I m n − E n , C n = P n R n . (36) R_n = I_m^n + \lambda D I_m^n - E_n, \qquad C_n = P_n R_n. \tag*{(36)} R n = I m n + λ D I m n − E n , C n = P n R n . (36) These are diagnostics of the existing iteration. The continuous-coordinate version of Eq. (17) , multiplied by λ \lambda λ ∫ Φ m R n d ϕ = 0 \int\Phi_m R_n\,d\phi= 0 ∫ Φ m R n d ϕ = 0 { ψ j } \{\psi_j\} { ψ j } V n V_n V n
C n = ∑ j = 1 dim V n ψ j ∫ 0 π ψ j R n d ϕ . (37) C_n = \sum_{j=1}^{\dim V_n} \psi_j \int_0^\pi\psi_j R_n\,d\phi. \tag*{(37)} C n = j = 1 ∑ d i m V n ψ j ∫ 0 π ψ j R n d ϕ . (37) Every row uses polar contractions of the retained modes and their streaming derivatives. This avoids inferring the complete residual from a single Galerkin condition. Finite spatial-sweep errors are included in C n C_n C n
Define the residual, approximation and source-increment indicators
η n = ∥ C n ∥ 3 , (38) \eta_n = \lVert C_n\rVert_3, \tag*{(38)} η n = ∥ C n ∥ 3 , (38) ζ n = ∥ ( 1 − P n ) E n ∥ 3 + λ ∥ ( 1 − P n ) D E n ∥ 2 , (39) \zeta_n = \lVert(1-P_n)E_n\rVert_3+\lambda\lVert(1-P_n)DE_n\rVert_2, \tag*{(39)} ζ n = ∥( 1 − P n ) E n ∥ 3 + λ ∥( 1 − P n ) D E n ∥ 2 , (39) ω n = ∥ U s n + 1 − U s n ∥ 2 , a n = η n + ζ n + ω n . (40) \omega_n = \lVert U_s^{n+1}-U_s^n\rVert_2,\qquad a_n=\eta_n+\zeta_n+\omega_n. \tag*{(40)} ω n = ∥ U s n + 1 − U s n ∥ 2 , a n = η n + ζ n + ω n . (40) The quantity ζ n \zeta_n ζ n G 2 = span { 1 , sin ϕ , sin 2 ϕ } G_2=\operatorname{span}\{1,\sin\phi,\sin^2\phi\} G 2 = span { 1 , sin ϕ , sin 2 ϕ } G 2 ⊂ V n G_2\subset V_n G 2 ⊂ V n ζ n = 0 \zeta_n=0 ζ n = 0 G 2 G_2 G 2 q s q_s q s O ( ε ) O(\varepsilon) O ( ε ) G 1 = span { 1 , sin ϕ } ⊂ V n G_1=\operatorname{span}\{1,\sin\phi\}\subset V_n G 1 = span { 1 , sin ϕ } ⊂ V n ζ n = O ( ε 2 ) \zeta_n=O(\varepsilon^2) ζ n = O ( ε 2 )
Bulk reconstruction and matching
Proposition 4 (Bulk reconstruction and matching). Let λ = Kn C → 0 \lambda=\operatorname{Kn}_{\mathrm{C}}\to0 λ = Kn C → 0 λ ≍ ε \lambda\asymp\varepsilon λ ≍ ε
∥ I m n ∥ 4 + ∥ I − n ∥ 4 + ∥ U s n ∥ 4 + ∥ U s n + 1 ∥ 4 ≤ C , (41) \lVert I_m^n\rVert_4+\lVert I_-^n\rVert_4+\lVert U_s^n\rVert_4+\lVert U_s^{n+1}\rVert_4\le C, \tag*{(41)} ∥ I m n ∥ 4 + ∥ I − n ∥ 4 + ∥ U s n ∥ 4 + ∥ U s n + 1 ∥ 4 ≤ C , (41) with C C C ε \varepsilon ε
∥ I m n − E n + λ D E n ∥ 2 ≤ C ( λ 2 + η n + ζ n ) , (42) \lVert I_m^n-E_n+\lambda DE_n\rVert_2\le C(\lambda^2+\eta_n+\zeta_n), \tag*{(42)} ∥ I m n − E n + λ D E n ∥ 2 ≤ C ( λ 2 + η n + ζ n ) , (42) and its kinetic moments match the refreshed synthetic fields with
∥ K [ I m n ] − U s n + 1 ∥ 2 ≤ C ( λ 2 + a n ) . (43) \lVert K[I_m^n]-U_s^{n+1}\rVert_2\le C(\lambda^2+a_n). \tag*{(43)} ∥ K [ I m n ] − U s n + 1 ∥ 2 ≤ C ( λ 2 + a n ) . (43) Moreover,
∥ M [ I m n ] ∥ 2 ≤ C ( λ 2 + η n + ζ n ) , ∥ λ ∇ ⋅ H [ I m n ] ∥ 0 ≤ C λ ( λ 2 + η n + ζ n ) . (44) \lVert M[I_m^n]\rVert_2\le C(\lambda^2+\eta_n+\zeta_n),\qquad
\lVert\lambda\nabla\cdot H[I_m^n]\rVert_0\le C\lambda(\lambda^2+\eta_n+\zeta_n). \tag*{(44)} ∥ M [ I m n ] ∥ 2 ≤ C ( λ 2 + η n + ζ n ) , ∥ λ ∇ ⋅ H [ I m n ] ∥ 0 ≤ C λ ( λ 2 + η n + ζ n ) . (44) No bound on the number of pairs needed to make the indicators small is asserted.
Proof. Since P n I m n = I m n P_nI_m^n=I_m^n P n I m n = I m n
I m n = P n E n − λ P n D I m n + C n . (45) I_m^n=P_nE_n-\lambda P_nDI_m^n+C_n. \tag*{(45)} I m n = P n E n − λ P n D I m n + C n . (45) Write d n = ( 1 − P n ) E n d_n=(1-P_n)E_n d n = ( 1 − P n ) E n D D D
I m n − E n + λ D E n = − d n + λ ( 1 − P n ) D E n + λ P n D d n + λ 2 P n D P n D I m n − λ P n D C n + C n . (46) \begin{aligned}
I_m^n-E_n+\lambda DE_n={}&-d_n+\lambda(1-P_n)DE_n+\lambda P_nDd_n\\
&+\lambda^2P_nDP_nDI_m^n-\lambda P_nDC_n+C_n. \tag*{(46)}
\end{aligned} I m n − E n + λ D E n = − d n + λ ( 1 − P n ) D E n + λ P n D d n + λ 2 P n D P n D I m n − λ P n D C n + C n . (46) Boundedness of D : H s + 1 → H s D:H^{s+1}\to H^s D : H s + 1 → H s ∥ P n ∥ ≤ 1 \lVert P_n\rVert\le1 ∥ P n ∥ ≤ 1
The equilibrium is a P 1 P_1 P 1 D E n DE_n D E n K i j k K_{ijk} K ij k ∇ ⋅ H \nabla\cdot\mathbf{H} ∇ ⋅ H
Taking energy and flux moments of (A.9) gives
T m k i n , n = T s n − λ χ ∇ ⋅ q s n + O ( λ 2 + η n + ζ n ) , T_{m}^{\mathrm{kin},n}=T_{s}^{n}-\lambda\chi\nabla\cdot\mathbf{q}_{s}^{n}+O(\lambda^{2}+\eta_{n}+\zeta_{n}), T m kin , n = T s n − λ χ ∇ ⋅ q s n + O ( λ 2 + η n + ζ n ) , q m k i n , n = χ q s n − λ 3 ∇ T s n + O ( λ 2 + η n + ζ n ) . (47) \mathbf{q}_{m}^{\mathrm{kin},n}=\chi\mathbf{q}_{s}^{n}-\frac{\lambda}{3}\nabla T_{s}^{n}+O(\lambda^{2}+\eta_{n}+\zeta_{n}). \tag*{(47)} q m kin , n = χ q s n − 3 λ ∇ T s n + O ( λ 2 + η n + ζ n ) . (47) where the remainders are in H 2 ( D 0 ) H^{2}(D_{0}) H 2 ( D 0 ) ∇ ⋅ q s n = 0 \nabla\cdot\mathbf{q}_{s}^{n}=0 ∇ ⋅ q s n = 0
1 3 ∇ T s n + δ R q s n = λ [ 1 5 ∇ 2 q s n + ∇ ⋅ H [ I − n ] ] . (48) \frac{1}{3}\nabla T_{s}^{n}+\delta_{R}\mathbf{q}_{s}^{n}
=\lambda\left[\frac{1}{5}\nabla^{2}\mathbf{q}_{s}^{n}+\nabla\cdot\mathbf{H}[I_{-}^{n}]\right]. \tag*{(48)} 3 1 ∇ T s n + δ R q s n = λ [ 5 1 ∇ 2 q s n + ∇ ⋅ H [ I − n ] ] . (48) Using 1 − χ = λ δ R 1-\chi=\lambda\delta_{R} 1 − χ = λ δ R
q m k i n , n − q s n = − λ 2 [ 1 5 ∇ 2 q s n + ∇ ⋅ H [ I − n ] ] + O ( λ 2 + η n + ζ n ) . (49) \mathbf{q}_{m}^{\mathrm{kin},n}-\mathbf{q}_{s}^{n}
=-\lambda^{2}\left[\frac{1}{5}\nabla^{2}\mathbf{q}_{s}^{n}+\nabla\cdot\mathbf{H}[I_{-}^{n}]\right]+O(\lambda^{2}+\eta_{n}+\zeta_{n}). \tag*{(49)} q m kin , n − q s n = − λ 2 [ 5 1 ∇ 2 q s n + ∇ ⋅ H [ I − n ] ] + O ( λ 2 + η n + ζ n ) . (49) The bracket is uniformly bounded in H 2 H^{2} H 2
Diffusive limit
The diffusive scaling is K n R = ε r ˉ \mathrm{Kn}_{R}=\varepsilon\bar{r} Kn R = ε r ˉ K n N = ε n ˉ \mathrm{Kn}_{N}=\varepsilon\bar{n} Kn N = ε n ˉ r ˉ , n ˉ \bar{r},\bar{n} r ˉ , n ˉ K n C = ε c ˉ \mathrm{Kn}_{C}=\varepsilon\bar{c} Kn C = ε c ˉ c ˉ = r ˉ n ˉ / ( r ˉ + n ˉ ) \bar{c}=\bar{r}\bar{n}/(\bar{r}+\bar{n}) c ˉ = r ˉ n ˉ / ( r ˉ + n ˉ ) χ = r ˉ / ( r ˉ + n ˉ ) < 1 \chi=\bar{r}/(\bar{r}+\bar{n})<1 χ = r ˉ / ( r ˉ + n ˉ ) < 1 I ( 0 ) = T ( 0 ) / ( 4 π ) + 3 χ Ω ⋅ q ( 0 ) / ( 4 π ) I^{(0)}=T^{(0)}/(4\pi)+3\chi\boldsymbol{\Omega}\cdot\mathbf{q}^{(0)}/(4\pi) I ( 0 ) = T ( 0 ) / ( 4 π ) + 3 χ Ω ⋅ q ( 0 ) / ( 4 π ) q ( 0 ) = χ q ( 0 ) \mathbf{q}^{(0)}=\chi\mathbf{q}^{(0)} q ( 0 ) = χ q ( 0 ) q ( 0 ) = 0 \mathbf{q}^{(0)}=0 q ( 0 ) = 0
q = − K n R 3 ∇ T + O ( ε 3 ) , ∇ ⋅ q = 0. (50) \mathbf{q}=-\frac{\mathrm{Kn}_{R}}{3}\nabla T+O(\varepsilon^{3}),\qquad\nabla\cdot\mathbf{q}=0. \tag*{(50)} q = − 3 Kn R ∇ T + O ( ε 3 ) , ∇ ⋅ q = 0. (50) The limiting bulk temperature is therefore harmonic. The coefficient is K n R / 3 \mathrm{Kn}_{R}/3 Kn R /3 K n C / 3 \mathrm{Kn}_{C}/3 Kn C /3
Corollary 1 (Diffusive limit). In the diffusive scaling, let the assumptions of Proposition 4 hold and let a n = O ( ε 1 + σ A P ) a_{n}=O(\varepsilon^{1+\sigma_{\mathrm{AP}}}) a n = O ( ε 1 + σ AP ) 0 < σ A P ≤ 1 0<\sigma_{\mathrm{AP}}\le1 0 < σ AP ≤ 1 L 2 ( D 0 ) L^{2}(D_{0}) L 2 ( D 0 )
q s n + 1 = − K n R 3 ∇ T s n + 1 + O ( ε 3 ) , \mathbf{q}_{s}^{n+1}=-\frac{\mathrm{Kn}_{R}}{3}\nabla T_{s}^{n+1}+O(\varepsilon^{3}), q s n + 1 = − 3 Kn R ∇ T s n + 1 + O ( ε 3 ) , q m k i n , n = − K n R 3 ∇ T s n + 1 + O ( ε 1 + σ A P ) . (51) \mathbf{q}_{m}^{\mathrm{kin},n}=-\frac{\mathrm{Kn}_{R}}{3}\nabla T_{s}^{n+1}+O(\varepsilon^{1+\sigma_{\mathrm{AP}}}). \tag*{(51)} q m kin , n = − 3 Kn R ∇ T s n + 1 + O ( ε 1 + σ AP ) . (51) The error in the kinetic heat flux is therefore o ( ε ) o(\varepsilon) o ( ε ) q m k i n , n / ε \mathbf{q}_{m}^{\mathrm{kin},n}/\varepsilon q m kin , n / ε − r ˉ / 3 -\bar{r}/3 − r ˉ /3
Proof. The refreshed synthetic momentum equation, in H 2 H^{2} H 2 q s n + 1 = O ( ε ) \mathbf{q}_{s}^{n+1}=O(\varepsilon) q s n + 1 = O ( ε ) K n R \mathrm{Kn}_{R} Kn R λ 2 = O ( ε 1 + σ A P ) \lambda^{2}=O(\varepsilon^{1+\sigma_{\mathrm{AP}}}) λ 2 = O ( ε 1 + σ AP )
Hydrodynamic limit
The normal-dominated scaling is K n N = ε n ˉ \mathrm{Kn}_{N}=\varepsilon\bar{n} Kn N = ε n ˉ K n R \mathrm{Kn}_{R} Kn R K n C = O ( ε ) \mathrm{Kn}_{C}=O(\varepsilon) Kn C = O ( ε ) χ → 1 \chi\to1 χ → 1 [ T + 3 Ω ⋅ q ] / ( 4 π ) [T+3\boldsymbol{\Omega}\cdot\mathbf{q}]/(4\pi) [ T + 3 Ω ⋅ q ] / ( 4 π ) N i j = − ( 2 / 5 ) K n C ∂ ⟨ j q i ⟩ + O ( ε 2 ) N_{ij}=-(2/5)\mathrm{Kn}_{C}\partial_{\langle j}q_{i\rangle}+O(\varepsilon^{2}) N ij = − ( 2/5 ) Kn C ∂ ⟨ j q i ⟩ + O ( ε 2 ) (7) . At fixed K n R \mathrm{Kn}_{R} Kn R ℓ G K → 0 \ell_{\mathrm{GK}}\to0 ℓ GK → 0
K n N = ε n ˉ , δ R = ε d ˉ R , T = T ∗ + ε p G K , q = O ( 1 ) , (52) \mathrm{Kn}_{N}=\varepsilon\bar{n},\qquad\delta_{R}=\varepsilon\bar{d}_{R},\qquad T=T_{*}+\varepsilon p_{\mathrm{GK}},\qquad\mathbf{q}=O(1), \tag*{(52)} Kn N = ε n ˉ , δ R = ε d ˉ R , T = T ∗ + ε p GK , q = O ( 1 ) , (52) where T ∗ T_{*} T ∗ d ˉ R ≥ 0 \bar{d}_{R}\ge0 d ˉ R ≥ 0 K n C / ε → n ˉ \mathrm{Kn}_{C}/\varepsilon\to\bar{n} Kn C / ε → n ˉ
1 3 ∇ p G K + d ˉ R q − n ˉ 5 ∇ 2 q = 0 , ∇ ⋅ q = 0. (53) \frac{1}{3}\nabla p_{\mathrm{GK}}+\bar{d}_{R}\mathbf{q}-\frac{\bar{n}}{5}\nabla^{2}\mathbf{q}=0,\qquad\nabla\cdot\mathbf{q}=0. \tag*{(53)} 3 1 ∇ p GK + d ˉ R q − 5 n ˉ ∇ 2 q = 0 , ∇ ⋅ q = 0. (53) This distinction separates a first-order GK correction from a limiting Brinkman/Stokes problem. Corollary 2 covers both.
Corollary 2 (Hydrodynamic limit). Let the assumptions of Proposition 4 hold and let a n = O ( ε 1 + σ A P ) a_{n}=O(\varepsilon^{1+\sigma_{\mathrm{AP}}}) a n = O ( ε 1 + σ AP ) 0 < σ A P ≤ 1 0<\sigma_{\mathrm{AP}}\le1 0 < σ AP ≤ 1 (7) through its viscous order. In the distinguished scaling (A.19), if the scaled temperature p G K p_{\mathrm{GK}} p GK
Proof. In the fixed-K n R \mathrm{Kn}_{R} Kn R O ( ε 3 + ε 2 + σ A P ) O(\varepsilon^{3}+\varepsilon^{2+\sigma_{\mathrm{AP}}}) O ( ε 3 + ε 2 + σ AP ) O ( ε 1 + σ A P ) O(\varepsilon^{1+\sigma_{\mathrm{AP}}}) O ( ε 1 + σ AP ) O ( ε 1 + σ A P ) O(\varepsilon^{1+\sigma_{\mathrm{AP}}}) O ( ε 1 + σ AP ) o ( ε ) o(\varepsilon) o ( ε ) ε \varepsilon ε λ / ε → n ˉ \lambda/\varepsilon\to\bar{n} λ / ε → n ˉ O ( ε σ A P ) O(\varepsilon^{\sigma_{\mathrm{AP}}}) O ( ε σ AP )
Ballistic limit
Let δ R , δ N → 0 \delta_{R},\delta_{N}\to0 δ R , δ N → 0 δ C → 0 \delta_{C}\to0 δ C → 0 K n C → ∞ \mathrm{Kn}_{C}\to\infty Kn C → ∞ d d d δ C d \delta_{C}d δ C d sin ϕ s ⋅ ∇ I = 0 \sin\phi\,\mathbf{s}\cdot\nabla I=0 sin ϕ s ⋅ ∇ I = 0 K n C \mathrm{Kn}_{C} Kn C § 4.3 ), degrades free streaming.
Proposition 5 (Ballistic limit). Let I I I U s = ( T s , q s ) U_{s}=(T_{s},\mathbf{q}_{s}) U s = ( T s , q s ) I b I_{b} I b I f s ( x , θ ) = I b ( x b ) I_{\mathrm{fs}}(x,\theta)=I_{b}(x_{b}) I fs ( x , θ ) = I b ( x b ) ( x , θ , ϕ ) (x,\theta,\phi) ( x , θ , ϕ )
∣ I − I f s ∣ ≤ [ ∣ I b ∣ + ∥ T s ∥ ∞ + 3 χ ∥ q s ∥ ∞ 4 π ] min ( 1 , δ C d sin ϕ ) , (54) \lvert I-I_{\mathrm{fs}}\rvert\le\left[\lvert I_{b}\rvert+\frac{\lVert T_{s}\rVert_{\infty}+3\chi\lVert\mathbf{q}_{s}\rVert_{\infty}}{4\pi}\right]\min\left(1,\frac{\delta_{C}d}{\sin\phi}\right), \tag*{(54)} ∣ I − I fs ∣ ≤ [ ∣ I b ∣ + 4 π ∥ T s ∥ ∞ + 3 χ ∥ q s ∥ ∞ ] min ( 1 , sin ϕ δ C d ) , (54) and, if every ray has in-plane length at most d max d_{\max} d m a x
∥ K [ I ] − K [ I f s ] ∥ ∞ ≤ 2 π 2 δ C d max [ ∥ I b ∥ ∞ + ∥ T s ∥ ∞ + 3 χ ∥ q s ∥ ∞ 4 π ] . (55) \left\lVert\mathcal{K}[I]-\mathcal{K}[I_{\mathrm{fs}}]\right\rVert_{\infty}\leq2\pi^{2}\delta_C d_{\max}\left[\left\lVert I_b\right\rVert_{\infty}+\frac{\left\lVert T_s\right\rVert_{\infty}+3\chi\left\lVert q_s\right\rVert_{\infty}}{4\pi}\right]. \tag*{(55)} ∥ K [ I ] − K [ I fs ] ∥ ∞ ≤ 2 π 2 δ C d m a x [ ∥ I b ∥ ∞ + 4 π ∥ T s ∥ ∞ + 3 χ ∥ q s ∥ ∞ ] . (55) Proof. By Eq. (12) , I − I f s = I b ( e − δ C d / sin ϕ − 1 ) + ∫ 0 d r [ A + sin ϕ B ] e − δ C t / sin ϕ d t / sin ϕ I-I_{\mathrm{fs}}=I_b(e^{-\delta_C d/\sin\phi}-1)+\int_0^{d_r}[A+\sin\phi B]e^{-\delta_C t/\sin\phi}\,\mathrm{d}t/\sin\phi I − I fs = I b ( e − δ C d / s i n ϕ − 1 ) + ∫ 0 d r [ A + sin ϕB ] e − δ C t / s i n ϕ d t / sin ϕ A A A ( ∣ I b ∣ + ∥ T s ∥ ∞ / ( 4 π ) ) ( 1 − e − δ C d / sin ϕ ) (|I_b|+\left\lVert T_s\right\rVert_{\infty}/(4\pi))(1-e^{-\delta_C d/\sin\phi}) ( ∣ I b ∣ + ∥ T s ∥ ∞ / ( 4 π )) ( 1 − e − δ C d / s i n ϕ ) δ N = χ δ C \delta_N=\chi\delta_C δ N = χ δ C B B B 3 χ ∥ q s ∥ ∞ min ( δ C d , sin ϕ ) / ( 4 π ) 3\chi\left\lVert q_s\right\rVert_{\infty}\min(\delta_C d,\sin\phi)/(4\pi) 3 χ ∥ q s ∥ ∞ min ( δ C d , sin ϕ ) / ( 4 π ) min ( 1 , δ C d / sin ϕ ) \min(1,\delta_C d/\sin\phi) min ( 1 , δ C d / sin ϕ ) d Ω = sin ϕ d ϕ d θ \mathrm{d}\Omega=\sin\phi\,\mathrm{d}\phi\,\mathrm{d}\theta d Ω = sin ϕ d ϕ d θ ∫ 0 π min ( 1 , τ / sin ϕ ) sin ϕ d ϕ ≤ π τ \int_0^\pi\min(1,\tau/\sin\phi)\sin\phi\,\mathrm{d}\phi\leq\pi\tau ∫ 0 π min ( 1 , τ / sin ϕ ) sin ϕ d ϕ ≤ π τ
Three consequences follow. First, for every bounded source the kinetic moments approach those of free streaming at the rate δ C d max \delta_C d_{\max} δ C d m a x δ C U s \delta_C U_s δ C U s K n C \mathrm{Kn}_C Kn C O ( K n C ) O(\mathrm{Kn}_C) O ( Kn C ) 3 k K n C 3k\mathrm{Kn}_C 3 k Kn C § 4.3 ). By (A.22) that error reaches the kinetic moments with gain O ( k d max ) O(kd_{\max}) O ( k d m a x ) K n C \mathrm{Kn}_C Kn C δ C = 0 \delta_C=0 δ C = 0 (16) is proportional to sin ϕ \sin\phi sin ϕ Φ m \Phi_m Φ m (17) carries the inflow along s s s sin ϕ ≲ δ C d \sin\phi\lesssim\delta_C d sin ϕ ≲ δ C d § 3.3 .
For thermalizing walls I b I_b I b (8) ; Eq. (A.21) then shows that I I I O ( δ C d ) O(\delta_C d) O ( δ C d ) I I I d max d_{\max} d m a x 2 π 2 δ C d max 2\pi^2\delta_Cd_{\max} 2 π 2 δ C d m a x π ∫ 0 2 π min ( 2 / π , δ C d ) d θ \pi\int_0^{2\pi}\min(2/\pi,\delta_Cd)\,\mathrm{d}\theta π ∫ 0 2 π min ( 2/ π , δ C d ) d θ O ( δ C ln δ C − 1 ) O(\delta_C\ln\delta_C^{-1}) O ( δ C ln δ C − 1 ) ρ ( A G ) → 4 / 5 \rho(A_G)\to4/5 ρ ( A G ) → 4/5 § 4.3 ).
Stopping criteria, parameters and scope
For the diffusive and hydrodynamic limits, a sufficient accuracy prescription is
η n + ζ n + ω n ≤ C A P ε 1 + σ A P , 0 < σ A P ≤ 1 ; (56) \eta_n+\zeta_n+\omega_n\leq C_{\mathrm{AP}}\varepsilon^{1+\sigma_{\mathrm{AP}}},\qquad0<\sigma_{\mathrm{AP}}\leq1; \tag*{(56)} η n + ζ n + ω n ≤ C AP ε 1 + σ AP , 0 < σ AP ≤ 1 ; (56) choosing σ A P = 1 \sigma_{\mathrm{AP}}=1 σ AP = 1 O ( ε 2 ) O(\varepsilon^2) O ( ε 2 ) o ( 1 ) o(1) o ( 1 ) O ( ε ) O(\varepsilon) O ( ε ) § 4.2 do not control (A.7); their fixed 10 − 3 10^{-3} 1 0 − 3
For the parametric construction the proposition applies node by node if its bounds and indicators hold uniformly at those nodes. The parameter-averaged spatial equation (28) does not imply these pointwise bounds; a small measure-weighted residual can hide an inaccurate small-Knudsen node. Uniform certification requires pointwise control of the indicators, including the parameter-dependent polar approximation.
The estimates require spatial regularity, not only bounded amplitudes: derivatives of an unresolved layer may scale as 1 / ε 1/\varepsilon 1/ ε O ( ε 2 ) O(\varepsilon^2) O ( ε 2 ) [21 , 12 ] , but asymptotic acceleration and preservation of the limiting constitutive law are distinct properties.
Fourier form of the alternating iteration
In a periodic homogeneous medium with an error mode e i k x e^{ikx} e ik x X r = Re ( x ^ r ( θ ) e i k x ) X_r=\operatorname{Re}(\widehat{x}_r(\theta)e^{ikx}) X r = Re ( x r ( θ ) e ik x ) ( T s , q s ) = Re ( ( T ^ , q ^ ) e i k x ) (T_s,\mathbf{q}_s)=\operatorname{Re}((\widehat{T},\widehat{\mathbf{q}})e^{ikx}) ( T s , q s ) = Re (( T , q ) e ik x ) ∫ Re ( a e i k x ) Re ( b e i k x ) d x = L 2 Re ( a b ‾ ) \int\operatorname{Re}(ae^{ikx})\operatorname{Re}(be^{ikx})\,\mathrm{d}x=\frac{L}{2}\operatorname{Re}(a\overline{b}) ∫ Re ( a e ik x ) Re ( b e ik x ) d x = 2 L Re ( a b ) L / 2 L/2 L /2
α m r = Re ∫ x ^ m x ^ r ‾ d θ , β m r = Re ∫ x ^ m i k cos θ x ^ r ‾ d θ , F T , m = δ C 4 π Re ∫ x ^ m T ^ ‾ d θ , F q , m = 3 δ N 4 π Re ∫ x ^ m s ⋅ q ^ ‾ d θ . (57) \begin{aligned}
\alpha_{mr} &= \operatorname{Re}\int\widehat{x}_m\overline{\widehat{x}_r}\,\mathrm{d}\theta, &
\beta_{mr} &= \operatorname{Re}\int\widehat{x}_m\overline{ik\cos\theta\,\widehat{x}_r}\,\mathrm{d}\theta,\\
F_{T,m} &= \frac{\delta_C}{4\pi}\operatorname{Re}\int\widehat{x}_m\overline{\widehat{T}}\,\mathrm{d}\theta, &
F_{q,m} &= \frac{3\delta_N}{4\pi}\operatorname{Re}\int\widehat{x}_m\overline{\mathbf{s}\cdot\widehat{\mathbf{q}}}\,\mathrm{d}\theta. \tag*{(57)}
\end{aligned} α m r F T , m = Re ∫ x m x r d θ , = 4 π δ C Re ∫ x m T d θ , β m r F q , m = Re ∫ x m ik cos θ x r d θ , = 4 π 3 δ N Re ∫ x m s ⋅ q d θ . (57) with β m m = 0 \beta_{mm}=0 β mm = 0 μ m = 0 \mu_m=0 μ m = 0
x ^ m ( θ ) = c m T ^ + d m s ⋅ q ^ − ∑ r < m ( a m r + i k b m r cos θ ) x ^ r ( θ ) a m m + i k b m m cos θ . (58) \widehat{x}_m(\theta)=\frac{c_m\widehat{T}+d_m\mathbf{s}\cdot\widehat{\mathbf{q}}-\sum_{r<m}(a_{mr}+ikb_{mr}\cos\theta)\widehat{x}_r(\theta)}{a_{mm}+ikb_{mm}\cos\theta}. \tag*{(58)} x m ( θ ) = a mm + ik b mm cos θ c m T + d m s ⋅ q − ∑ r < m ( a m r + ik b m r cos θ ) x r ( θ ) . (58) Each new spatial factor starts constant in θ \theta θ U 0 = ( 1 , 0.5 i , 0.5 ) U^0=(1,0.5i,0.5) U 0 = ( 1 , 0.5 i , 0.5 ) γ G \gamma_G γ G n 12 n_{12} n 12 ∥ e ∥ \lVert e\rVert ∥ e ∥ (11) . For the synthetic map, the third moments of the separated state need only M 2 M_2 M 2 M 4 M_4 M 4
M x x x = ∑ r ∫ x ^ r [ cos 3 θ M 4 ( Φ r ) − 3 5 cos θ M 2 ( Φ r ) ] d θ , M_{xxx}=\sum_r\int\widehat{x}_r\left[\cos^3\theta\,M_4(\Phi_r)-\frac{3}{5}\cos\theta\,M_2(\Phi_r)\right]\,\mathrm{d}\theta, M xxx = r ∑ ∫ x r [ cos 3 θ M 4 ( Φ r ) − 5 3 cos θ M 2 ( Φ r ) ] d θ , M y x x = ∑ r ∫ x ^ r [ cos 2 θ sin θ M 4 ( Φ r ) − 1 5 sin θ M 2 ( Φ r ) ] d θ . (59) M_{yxx}=\sum_r\int\widehat{x}_r\left[\cos^2\theta\sin\theta\,M_4(\Phi_r)-\frac{1}{5}\sin\theta\,M_2(\Phi_r)\right]\,\mathrm{d}\theta. \tag*{(59)} M y xx = r ∑ ∫ x r [ cos 2 θ sin θ M 4 ( Φ r ) − 5 1 sin θ M 2 ( Φ r ) ] d θ . (59) and substituting ∇ → i k e x \nabla\to ik\mathbf{e}_x ∇ → ik e x (6) with zero data gives q ^ x = 0 \widehat{q}_x=0 q x = 0 (25) from the two momentum equations.
The full-order operators A K A_K A K A G A_G A G I e x = ( δ C T ^ + 3 δ N Ω ⋅ q ^ ) / [ 4 π ( δ C + i k sin ϕ cos θ ) ] I_{\mathrm{ex}} = (\delta_C \hat{T} + 3\delta_N \boldsymbol{\Omega} \cdot\hat{\boldsymbol{q}})/[4\pi(\delta_C + ik\sin\phi\cos\theta)] I ex = ( δ C T ^ + 3 δ N Ω ⋅ q ^ ) / [ 4 π ( δ C + ik sin ϕ cos θ )] A K A_K A K [21 ] ; its spectral radius tends to one as k K n C → 0 k\mathrm{Kn}_C \to0 k Kn C → 0 M x x x = T ^ ⋅ 1 2 ∫ − 1 1 ( ξ 3 − 3 5 ξ ) ( 1 + i k ^ ξ ) − 1 d ξ → − 4 T ^ / ( 15 i k ^ ) M_{xxx} = \hat{T} \cdot\frac{1}{2}\int_{-1}^{1}(\xi^3-\frac{3}{5}\xi)(1+i\hat{k}\xi)^{-1}\,d\xi\to-4\hat{T}/(15i\hat{k}) M xxx = T ^ ⋅ 2 1 ∫ − 1 1 ( ξ 3 − 5 3 ξ ) ( 1 + i k ^ ξ ) − 1 d ξ → − 4 T ^ / ( 15 i k ^ ) k ^ = k K n C → ∞ \hat{k}=k\mathrm{Kn}_C \to\infty k ^ = k Kn C → ∞ ξ = Ω x \xi=\Omega_x ξ = Ω x (25) then gives T ^ ′ → − 4 5 T ^ \hat{T}' \to-\frac{4}{5}\hat{T} T ^ ′ → − 5 4 T ^ ρ ( A G ) → 4 / 5 \rho(A_G)\to4/5 ρ ( A G ) → 4/5 ρ ( A K ) = O ( 1 / k ^ ) \rho(A_K)=O(1/\hat{k}) ρ ( A K ) = O ( 1/ k ^ )
Spatial sweep of pure PGD. For a fixed polar factor Φ \Phi Φ a = δ C ∫ Φ 2 d ϕ a=\delta_C\int\Phi^2\,d\phi a = δ C ∫ Φ 2 d ϕ b = ∫ sin ϕ Φ 2 d ϕ b=\int\sin\phi\,\Phi^2\,d\phi b = ∫ sin ϕ Φ 2 d ϕ c T = δ C 4 π ∫ Φ d ϕ M 1 ( Φ ) c^T=\frac{\delta_C}{4\pi}\int\Phi\,d\phi\,M_1(\Phi) c T = 4 π δ C ∫ Φ d ϕ M 1 ( Φ ) d Q = 3 δ N 4 π M 1 ( Φ ) M 2 ( Φ ) d_Q=\frac{3\delta_N}{4\pi}M_1(\Phi)M_2(\Phi) d Q = 4 π 3 δ N M 1 ( Φ ) M 2 ( Φ ) (17) –(18) with the current-mode moments on the right. One lagged sweep maps the error x ^ \hat{x} x ^
x ^ ′ ( θ ) = c T Y ^ + d Q s ( θ ) ⋅ Z ^ a + i k b cos θ , Y ^ = ∫ 0 2 π x ^ d θ , Z ^ = ∫ 0 2 π s x ^ d θ (60) \hat{x}'(\theta)=\frac{c^T\hat{Y}+d_Q\,s(\theta)\cdot\hat{\boldsymbol{Z}}}{a+ikb\cos\theta},\qquad
\hat{Y}=\int_0^{2\pi}\hat{x}\,d\theta,\qquad
\hat{\boldsymbol{Z}}=\int_0^{2\pi}s\,\hat{x}\,d\theta \tag*{(60)} x ^ ′ ( θ ) = a + ik b cos θ c T Y ^ + d Q s ( θ ) ⋅ Z ^ , Y ^ = ∫ 0 2 π x ^ d θ , Z ^ = ∫ 0 2 π s x ^ d θ (60) which is a linear map of ( Y ^ , Z ^ x , Z ^ y ) (\hat{Y},\hat{Z}_x,\hat{Z}_y) ( Y ^ , Z ^ x , Z ^ y ) ρ s w \rho_{\mathrm{sw}} ρ sw Table 1 is the spectral radius of its 3 × 3 3\times3 3 × 3 § 6 . For Φ = 1 \Phi=1 Φ = 1 a / b = π δ C / 2 a/b=\pi\delta_C/2 a / b = π δ C /2 2 π c T / a 2\pi c^T/a 2 π c T / a π d Q / a \pi d_Q/a π d Q / a 3 χ / 4 3\chi/4 3 χ /4 Φ = sin ϕ \Phi=\sin\phi Φ = sin ϕ a / b = 3 π δ C / 8 a/b=3\pi\delta_C/8 a / b = 3 π δ C /8 χ \chi χ k / δ C → 0 k/\delta_C\to0 k / δ C → 0 Φ = sin ϕ \Phi=\sin\phi Φ = sin ϕ χ → 1 \chi\to1 χ → 1 β m m = 0 \beta_{mm}=0 β mm = 0 x ^ m \hat{x}_m x ^ m 2 α m m T = α m m q 2\alpha^T_{mm}=\alpha^q_{mm} 2 α mm T = α mm q α m m q = 0 \alpha^q_{mm}=0 α mm q = 0 (19) then maps every Φ \Phi Φ sin ϕ \sin\phi sin ϕ β m m \beta_{mm} β mm
DG streaming form
Let V h V_h V h p p p K K K [40 , 41 ] , and a ± = ( a ± ∣ a ∣ ) / 2 a^\pm=(a\pm|a|)/2 a ± = ( a ± ∣ a ∣ ) /2
b h ( V , U ) = ∑ K ∫ 0 2 π [ − ( s ⋅ ∇ V , U ) K + ⟨ V , ( s ⋅ n ) + U ⟩ ∂ K + ⟨ V , ( s ⋅ n ) − U e x t ⟩ ∂ K ∩ F i ] d θ (61) b_h(V,U)=\sum_K\int_0^{2\pi}\left[-(s\cdot\nabla V,U)_K+\langle V,(s\cdot n)^+U\rangle_{\partial K}+\langle V,(s\cdot n)^-U^{\mathrm{ext}}\rangle_{\partial K\cap\mathcal{F}_i}\right]\,d\theta \tag*{(61)} b h ( V , U ) = K ∑ ∫ 0 2 π [ − ( s ⋅ ∇ V , U ) K + ⟨ V , ( s ⋅ n ) + U ⟩ ∂ K + ⟨ V , ( s ⋅ n ) − U ext ⟩ ∂ K ∩ F i ] d θ (61) where F i \mathcal{F}_i F i U e x t U^{\mathrm{ext}} U ext V = U V=U V = U
b h ( U , U ) = 1 2 ∫ 0 2 π [ ∑ F ∈ F i ∣ s ⋅ n ∣ [ [ U ] ] 2 + ∫ ∂ D w ∣ s ⋅ n ∣ U 2 ] d θ ≥ 0 , b_h(U,U)=\frac{1}{2}\int_0^{2\pi}\left[\sum_{F\in\mathcal{F}_i}|s\cdot n|[[U]]^2+\int_{\partial D_w}|s\cdot n|U^2\right]\,d\theta\ge0, b h ( U , U ) = 2 1 ∫ 0 2 π [ F ∈ F i ∑ ∣ s ⋅ n ∣ [[ U ] ] 2 + ∫ ∂ D w ∣ s ⋅ n ∣ U 2 ] d θ ≥ 0 , with ∂ D w \partial D_w ∂ D w § 4.1 . The thermal-wall load is g b ( V ) = ∫ ∂ D w ∫ 0 2 π V ( s ⋅ n ) − T b / ( 4 π ) d θ d S g_b(V)=\int_{\partial D_w}\int_0^{2\pi}V(s\cdot n)^-T_b/(4\pi)\,d\theta\,dS g b ( V ) = ∫ ∂ D w ∫ 0 2 π V ( s ⋅ n ) − T b / ( 4 π ) d θ d S β m r = b h ( X m , X r ) \beta_{mr}=b_h(X_m,X_r) β m r = b h ( X m , X r ) μ m = g b ( X m ) \mu_m=g_b(X_m) μ m = g b ( X m ) (16) ; the adiabatic load replaces T b / ( 4 π ) T_b/(4\pi) T b / ( 4 π ) F i \mathcal{F}_i F i Δ T b / ( 4 π ) \Delta T_b/(4\pi) Δ T b / ( 4 π ) θ j \theta_j θ j X r ∣ K = ∑ l x r , l K , j φ l X_r|_K=\sum_l x^{K,j}_{r,l}\varphi_l X r ∣ K = ∑ l x r , l K , j φ l (17) becomes
( a m m M K + b m m B K , j ) x m K , j = M K ( c m T s , K + d m s j ⋅ q s , K ) − T m 4 π f K , j b − b m m ∑ F C K F , j x m e x t , j (62) (a_{mm}M_K+b_{mm}B_{K,j})x^{K,j}_m=M_K(c_mT_{s,K}+d_ms_j\cdot q_{s,K})-\frac{T_m}{4\pi}f^b_{K,j}-b_{mm}\sum_F C_{KF,j}x^{\mathrm{ext},j}_m \tag*{(62)} ( a mm M K + b mm B K , j ) x m K , j = M K ( c m T s , K + d m s j ⋅ q s , K ) − 4 π T m f K , j b − b mm F ∑ C K F , j x m ext , j (62) − ∑ r = 1 m − 1 [ ( a m r M K + b m r B K , j ) x r K , j + b m r ∑ F C K F , j x r e x t , j ] , -\sum_{r=1}^{m-1}\left[(a_{mr}M_K+b_{mr}B_{K,j})x^{K,j}_r+b_{mr}\sum_F C_{KF,j}x^{\mathrm{ext},j}_r\right], − r = 1 ∑ m − 1 [ ( a m r M K + b m r B K , j ) x r K , j + b m r F ∑ C K F , j x r ext , j ] , written for PGD–GSIS (pure PGD replaces T s , K T_{s,K} T s , K q s , K q_{s,K} q s , K M K M_K M K B K , j B_{K,j} B K , j C K F , j C_{KF,j} C K F , j F F F f K , j b f^{b}_{K,j} f K , j b ⟨ φ l , ( s j ⋅ n ) − T b ⟩ ∂ K ∩ ∂ D \langle\varphi_l,(s_j\cdot n)^{-}T_b\rangle_{\partial K\cap\partial D} ⟨ φ l , ( s j ⋅ n ) − T b ⟩ ∂ K ∩ ∂ D a ˉ \bar{a} a ˉ b ˉ \bar{b} b ˉ S T , m S_{T,m} S T , m S q , m S_{q,m} S q , m
HDG discretization of the synthetic system
Mixed form. With the auxiliary tensor L L L Π \boldsymbol{\Pi} Π
Π x x = 4 3 L x x − 2 3 L y y , Π x y = Π y x = L x y + L y x , Π y y = 4 3 L y y − 2 3 L x x , (63) \Pi_{xx}=\frac{4}{3}L_{xx}-\frac{2}{3}L_{yy},\qquad\Pi_{xy}=\Pi_{yx}=L_{xy}+L_{yx},\qquad\Pi_{yy}=\frac{4}{3}L_{yy}-\frac{2}{3}L_{xx}, \tag*{(63)} Π xx = 3 4 L xx − 3 2 L y y , Π x y = Π y x = L x y + L y x , Π y y = 3 4 L y y − 3 2 L xx , (63) (6) is written on each element as
∇ ⋅ q = 0 , 1 3 ∇ T + δ R q − ∇ ⋅ Π = 0 , L − K n C 5 ∇ q = K n C H , (64) \nabla\cdot\boldsymbol{q}=0,\qquad\frac{1}{3}\nabla T+\delta_R\boldsymbol{q}-\nabla\cdot\boldsymbol{\Pi}=0,\qquad L-\frac{\mathrm{Kn}_C}{5}\nabla\boldsymbol{q}=\mathrm{Kn}_C H, \tag*{(64)} ∇ ⋅ q = 0 , 3 1 ∇ T + δ R q − ∇ ⋅ Π = 0 , L − 5 Kn C ∇ q = Kn C H , (64) where L i j L_{ij} L ij K n C 5 ∂ j q i \frac{\mathrm{Kn}_C}{5}\partial_jq_i 5 Kn C ∂ j q i H = ( H x x + 1 2 H y y , 1 2 H x y , 1 2 H x y , 1 2 H x x + H y y ) H=(H_{xx}+\frac{1}{2}H_{yy},\frac{1}{2}H_{xy},\frac{1}{2}H_{xy},\frac{1}{2}H_{xx}+H_{yy}) H = ( H xx + 2 1 H y y , 2 1 H x y , 2 1 H x y , 2 1 H xx + H y y ) ( L x x , L x y , L y x , L y y ) (L_{xx},L_{xy},L_{yx},L_{yy}) ( L xx , L x y , L y x , L y y ) ∂ z = 0 \partial_z=0 ∂ z = 0 K n C 5 [ ∇ 2 q + 1 3 ∇ ( ∇ ⋅ q ) ] \frac{\mathrm{Kn}_C}{5}[\nabla^2\boldsymbol{q}+\frac{1}{3}\nabla(\nabla\cdot\boldsymbol{q})] 5 Kn C [ ∇ 2 q + 3 1 ∇ ( ∇ ⋅ q )] H H H Π ( H ) = H \boldsymbol{\Pi}(H)=H Π ( H ) = H ∇ ⋅ Π \nabla\cdot\boldsymbol{\Pi} ∇ ⋅ Π K n C ∇ ⋅ H \mathrm{Kn}_C\nabla\cdot H Kn C ∇ ⋅ H u K = ( T , q x , q y , L x x , L x y , L y x , L y y ) ∈ [ P p ( K ) ] 7 \boldsymbol{u}_K=(T,q_x,q_y,L_{xx},L_{xy},L_{yx},L_{yy})\in[P_p(K)]^7 u K = ( T , q x , q y , L xx , L x y , L y x , L y y ) ∈ [ P p ( K ) ] 7 u ^ F = ( T ^ , q ^ x , q ^ y ) ∈ [ P p ( F ) ] 3 \hat{\boldsymbol{u}}_F=(\hat{T},\hat{q}_x,\hat{q}_y)\in[P_p(F)]^3 u ^ F = ( T ^ , q ^ x , q ^ y ) ∈ [ P p ( F ) ] 3
Local problem. With stabilization τ T = τ q = 1 \tau_T=\tau_q=1 τ T = τ q = 1 φ ∈ P p ( K ) \varphi\in P_p(K) φ ∈ P p ( K )
− ( q , ∇ φ ) K + ⟨ q ⋅ n + τ T T , φ ⟩ ∂ K = ⟨ τ T T ^ , φ ⟩ ∂ K , (65) -(\boldsymbol{q},\nabla\varphi)_K+\langle\boldsymbol{q}\cdot n+\tau_TT,\varphi\rangle_{\partial K}=\langle\tau_T\hat{T},\varphi\rangle_{\partial K}, \tag*{(65)} − ( q , ∇ φ ) K + ⟨ q ⋅ n + τ T T , φ ⟩ ∂ K = ⟨ τ T T ^ , φ ⟩ ∂ K , (65) δ R ( q i , φ ) K − 1 3 ( T , ∂ i φ ) K − ( φ , ∂ k Π i k ) K + ⟨ τ q q i , φ ⟩ ∂ K = ⟨ τ q q ^ i − 1 3 n i T ^ , φ ⟩ ∂ K , (66) \delta_R(q_i,\varphi)_K-\frac{1}{3}(T,\partial_i\varphi)_K-(\varphi,\partial_k\Pi_{ik})_K+\langle\tau_qq_i,\varphi\rangle_{\partial K}=\langle\tau_q\hat{q}_i-\frac{1}{3}n_i\hat{T},\varphi\rangle_{\partial K}, \tag*{(66)} δ R ( q i , φ ) K − 3 1 ( T , ∂ i φ ) K − ( φ , ∂ k Π ik ) K + ⟨ τ q q i , φ ⟩ ∂ K = ⟨ τ q q ^ i − 3 1 n i T ^ , φ ⟩ ∂ K , (66) ( L i j , φ ) K + K n C 5 ( q i , ∂ j φ ) K = K n C 5 ⟨ n j q ^ i , φ ⟩ ∂ K + K n C ( H i j , φ ) K , (67) (L_{ij},\varphi)_K+\frac{\mathrm{Kn}_C}{5}(q_i,\partial_j\varphi)_K=\frac{\mathrm{Kn}_C}{5}\langle n_j\hat{q}_i,\varphi\rangle_{\partial K}+\mathrm{Kn}_C(H_{ij},\varphi)_K, \tag*{(67)} ( L ij , φ ) K + 5 Kn C ( q i , ∂ j φ ) K = 5 Kn C ⟨ n j q ^ i , φ ⟩ ∂ K + Kn C ( H ij , φ ) K , (67) where H i j H_{ij} H ij A K u K = B K u ^ ∂ K + f K [ I m ] A_K\boldsymbol{u}_K=B_K\hat{\boldsymbol{u}}_{\partial K}+f_K[I_m] A K u K = B K u ^ ∂ K + f K [ I m ] A K A_K A K δ R \delta_R δ R K n C \mathrm{Kn}_C Kn C A K − 1 B K A_K^{-1}B_K A K − 1 B K
Global problem. Single-valuedness of the numerical energy flux q ⋅ n + τ T ( T − T ^ ) \boldsymbol{q}\cdot n+\tau_T(T-\hat{T}) q ⋅ n + τ T ( T − T ^ ) 1 3 T ^ n − Π n + τ q ( q − q ^ ) \frac{1}{3}\hat{T}n-\boldsymbol{\Pi}n+\tau_q(\boldsymbol{q}-\hat{\boldsymbol{q}}) 3 1 T ^ n − Π n + τ q ( q − q ^ ) u K \boldsymbol{u}_K u K K U ^ = F [ I m ] K\hat{U}=F[I_m] K U ^ = F [ I m ] [42 ] ; every synthetic refresh reassembles F F F Δ T b \Delta T_b Δ T b F F F
1 ∣ D ∣ ∫ D T d x = T ∗ , (68) \frac{1}{|D|}\int_D T\,\mathrm{d}x=T_*, \tag*{(68)} ∣ D ∣ 1 ∫ D T d x = T ∗ , (68) through a Lagrange multiplier.
Wall closure. On a physical face the conservation condition is replaced by a linearized kinetic wall law [19 , 21 ] .
Let U = ( T , q x , q y ) U=(T,q_x,q_y) U = ( T , q x , q y ) U k i n U^{\mathrm{kin}} U kin I m I_m I m ψ ( Ω ) = ( Ω ⋅ n ) ( 1 , Ω x , Ω y ) \psi(\boldsymbol{\Omega})=(\boldsymbol{\Omega}\cdot\boldsymbol{n})(1,\Omega_x,\Omega_y) ψ ( Ω ) = ( Ω ⋅ n ) ( 1 , Ω x , Ω y ) I e q ( U ) = a ( Ω ) ⋅ U I^{\mathrm{eq}}(U)=\boldsymbol{a}(\boldsymbol{\Omega})\cdot U I eq ( U ) = a ( Ω ) ⋅ U a = ( 1 , 3 χ Ω x , 3 χ Ω y ) / ( 4 π ) \boldsymbol{a}=(1,3\chi\Omega_x,3\chi\Omega_y)/(4\pi) a = ( 1 , 3 χ Ω x , 3 χ Ω y ) / ( 4 π ) J F + = ∫ Ω ⋅ n > 0 ψ a T d Ω J_F^+=\int_{\boldsymbol{\Omega}\cdot\boldsymbol{n}>0}\psi\boldsymbol{a}^{\mathsf{T}}\,\mathrm{d}\Omega J F + = ∫ Ω ⋅ n > 0 ψ a T d Ω J F = J F + J_F=J_F^+ J F = J F +
On a diffuse adiabatic wall the incoming level (8) is linear in U U U J F − = ( ∫ Ω ⋅ n < 0 ψ d Ω ) e + T / ( − D ) J_F^-=(\int_{\boldsymbol{\Omega}\cdot\boldsymbol{n}<0}\psi\,\mathrm{d}\Omega)e_+^{\mathsf{T}}/(-D) J F − = ( ∫ Ω ⋅ n < 0 ψ d Ω ) e + T / ( − D ) e + T e_+^{\mathsf{T}} e + T J F + J_F^+ J F + D = ∫ Ω ⋅ n < 0 Ω ⋅ n d Ω D=\int_{\boldsymbol{\Omega}\cdot\boldsymbol{n}<0}\boldsymbol{\Omega}\cdot\boldsymbol{n}\,\mathrm{d}\Omega D = ∫ Ω ⋅ n < 0 Ω ⋅ n d Ω
F ^ ( U ) = F k i n [ I m ] + J F ( U − U k i n ) , F k i n [ I m ] = ∫ ψ ( Ω ) I m u p d Ω , (69) \widehat{\mathcal{F}}(U)=\mathcal{F}^{\mathrm{kin}}[I_m]+J_F(U-U^{\mathrm{kin}}),\qquad\mathcal{F}^{\mathrm{kin}}[I_m]=\int\psi(\boldsymbol{\Omega})I_m^{\mathrm{up}}\,\mathrm{d}\Omega, \tag*{(69)} F ( U ) = F kin [ I m ] + J F ( U − U kin ) , F kin [ I m ] = ∫ ψ ( Ω ) I m up d Ω , (69) with I m u p I_m^{\mathrm{up}} I m up U = U k i n U=U^{\mathrm{kin}} U = U kin I m I_m I m J F J_F J F H \mathbf{H} H F k i n \mathcal{F}^{\mathrm{kin}} F kin I m I_m I m T s y n T^{\mathrm{syn}} T syn q s y n \boldsymbol{q}^{\mathrm{syn}} q syn
