Snapshot from this documentation build. SHA-256: a822aa016d33bc6a989c2548e959bed4e005b819b5a0edd65178221e59b764ad
1 using LinearAlgebra
2 using Test
3
4 @testset "spectral radiation coupling" begin
5 @testset "robust PLM limiter" begin
6 L=ntuple(_->0f0,9)
7 M=ntuple(_->1f0,9)
8 R=ntuple(_->3f0,9)
9 moncen=MHDKernels.prim_slope(L,M,R)
10 limited=MHDKernels.prim_slope_minmod(L,M,R)
11 @test all(limited .== 0.5f0)
12 @test all(limited .<= moncen)
13
14 turning=ntuple(_->0f0,9)
15 @test all(MHDKernels.prim_slope_minmod(L,M,turning) .== 0f0)
16 @test MHDKernels.recon_code(:plm_minmod) ==
17 MHDKernels.RECON_PLM_MINMOD
18 end
19
20 @testset "authoritative density perturbation" begin
21 delta=zeros(Float32,1)
22 rho=1f0
23 for _ in 1:32
24 rho=MHDKernels._cube_density_update!(delta,1,rho,1f-8,Val(true))
25 end
26 @test delta[1] ≈ 3.2f-7 rtol=2f-6
27 @test rho == 1f0
28
29 s=allocate_state(backend(:cpu),Float32,(4,4,4);dx=0.25)
30 fill!(s.U[1],1f0)
31 ws=allocate_radiation_workspace(s)
32 ws.density_perturbation[1]=4f-8
33 reconstruct_radiation_density!(ws,s)
34 @test ws.density_perturbation[1] == 4f-8
35 @test s.U[1][1] == 1f0
36
37 fill!(ws.density_perturbation,0f0)
38 fill!(ws.fft_real,1f-8)
39 foreach(x->fill!(x,0f0),ws.correction)
40 for _ in 1:32
41 MHDKernels._radiation_apply_correction!(s,ws,ws.correction)
42 end
43 @test ws.density_perturbation[1] ≈ 3.2f-7 rtol=2f-6
44 @test s.U[1][1] == 1f0+Float32(3.2f-7)
45
46 fill!(ws.fft_real,4f-8)
47 MHDKernels._radiation_apply_correction!(s,ws,ws.correction;
48 density_absolute=true)
49 @test ws.density_perturbation[1] == 4f-8
50 @test s.U[1][1] == 1f0
51 end
52
53 @testset "authoritative density half spectrum" begin
54 N=8
55 s=allocate_state(backend(:cpu),Float32,(N,N,N);dx=1/N)
56 fill!(s.U[1],1f0)
57 ws=allocate_radiation_workspace(s)
58 @test !ws.density_hat_valid[]
59 @inbounds for k in 1:N,j in 1:N,i in 1:N
60 ws.density_perturbation[i,j,k]=
61 1f-2*cospi(2f0*Float32(i-1)/Float32(N))
62 end
63 initial=copy(MHDKernels._radiation_ensure_density_hat!(ws,s))
64 @test ws.density_hat_valid[]
65
66 # Once valid, imported/compatibility real-space round-off cannot replace
67 # the authoritative spectral state without an explicit refresh.
68 fill!(ws.density_perturbation,0f0)
69 MHDKernels._radiation_ensure_density_hat!(ws,s)
70 @test ws.authoritative_density_hat==initial
71
72 copyto!(ws.density_hat,2f0.*initial)
73 MHDKernels._radiation_commit_density_hat!(ws,s;
74 additive=false,density_blend=0.25f0)
75 @test ws.authoritative_density_hat ≈ 1.25f0.*initial rtol=2f-6
76
77 copyto!(ws.density_hat,initial)
78 MHDKernels._radiation_commit_density_hat!(ws,s;
79 additive=true,density_blend=0.5f0)
80 @test ws.authoritative_density_hat ≈ 1.75f0.*initial rtol=2f-6
81
82 fill!(s.U[1],1.01f0)
83 MHDKernels._radiation_refresh_density!(ws,s)
84 @test !ws.density_hat_valid[]
85 refreshed=MHDKernels._radiation_ensure_density_hat!(ws,s)
86 @test abs(refreshed[1]-ComplexF32(0.01f0*N^3)) < 2f-5*N^3
87 end
88
89 @testset "positivity blend scales the coupled correction" begin
90 s=allocate_state(backend(:cpu),Float32,(4,4,4);dx=0.25)
91 fill!(s.U[1],1f0); fill!(s.U[2],0f0); fill!(s.U[3],0f0)
92 fill!(s.U[4],0f0); fill!(s.U[5],1f0)
93 foreach(x->fill!(x,0f0),s.U[6:9])
94 ws=allocate_radiation_workspace(s)
95 fill!(ws.density_perturbation,0f0)
96 fill!(ws.fft_real,-0.4f0)
97 fill!(ws.correction[1],2f0)
98 fill!(ws.correction[2],0f0)
99 fill!(ws.correction[3],0f0)
100
101 MHDKernels._radiation_apply_correction!(s,ws,ws.correction;
102 density_blend=0.25f0)
103
104 @test all(s.U[1].==0.9f0)
105 @test all(ws.density_perturbation.==-0.1f0)
106 @test all(s.U[2]./s.U[1].==0.5f0)
107 end
108
109 @testset "matched asymptotic coefficients" begin
110 rc=RadiationClosure(inertia_ratio=0.2f0,mean_free_path=0.01f0,
111 photon_speed=20f0,gas_sound_speed=1f0)
112 low=radiation_mode_coefficients(1f-3,7f0,rc)
113 high=radiation_mode_coefficients(1f8,7f0,rc)
114 @test low.streaming_weight < 1f-9
115 @test low.inertia_response ≈ 1f0/6f0 rtol=2f-6
116 @test low.radiation_sound_speed2 ≈ 400f0/3.6f0 rtol=2f-6
117 @test high.streaming_weight > 0.999999
118 @test high.streaming_weight == 1
119 @test high.inertia_response ≈ 1 rtol=2f-6
120 @test high.gamma_transverse == 7
121 @test high.radiation_sound_speed2 == 0
122 end
123
124 @testset "angular transport response" begin
125 q1=radiation_angular_response(1.0)
126 @test q1.monopole ≈ π/4 rtol=2e-7
127 @test q1.longitudinal ≈ 3*(1-π/4) rtol=2e-7
128 @test q1.transverse ≈ (3π/4-q1.longitudinal)/2 rtol=2e-7
129
130 R=0.2f0; lambda=0.01f0; c=20f0
131 rc=RadiationClosure(inertia_ratio=R,mean_free_path=lambda,
132 photon_speed=c,gas_sound_speed=1f0,
133 response_model=:angular)
134 gamma=c/(R*lambda)
135 tight=radiation_mode_coefficients(1f0,gamma,rc)
136 expected=c*lambda/(5f0*(1f0+R))
137 @test tight.response_model === :angular
138 @test tight.gamma_transverse ≈ expected rtol=3f-3
139 @test tight.gamma_longitudinal/tight.gamma_transverse ≈ 4f0/3f0 rtol=2f-4
140 @test tight.radiation_sound_speed2 ≈ c*c/(3f0*(1f0+R)) rtol=2f-3
141
142 free=radiation_mode_coefficients(1f8,7f0,rc)
143 @test free.streaming_weight > 0.998f0
144 @test free.inertia_response > 0.999f0
145 @test free.gamma_transverse ≈ 7f0 rtol=3f-3
146 @test free.gamma_longitudinal ≈ 7f0 rtol=3f-3
147 @test free.radiation_sound_speed2 < 1f-8
148 end
149
150 @testset "stiff acoustic propagator retains the slow mode in Float32" begin
151 g=1f8; c2=1f0; k2=1f8; dt=0.25f0
152 B32,S32=MHDKernels._radiation_acoustic_coefficients(g,c2,k2,dt)
153 B64,S64=MHDKernels._radiation_acoustic_coefficients(
154 Float64(g),Float64(c2),Float64(k2),Float64(dt))
155 @test isfinite(B32) && isfinite(S32)
156 @test B32 < 0 && S32 > 0
157 @test B32 ≈ Float32(B64) rtol=2f-6
158 @test S32 ≈ Float32(S64) rtol=2f-6
159 end
160
161 @testset "photon exponential split has an exact zero-step identity" begin
162 ub=ComplexF32(0.01f0,0.02f0)
163 ug=ComplexF32(-0.03f0,0.04f0)
164 force=ComplexF32(0.2f0,-0.1f0)
165 ubn,ugn=MHDKernels._photon_exact_drag_affine(
166 ub,ug,force,0f0,13f0,3f0,1500f0)
167 @test ubn == ub
168 @test ugn == ug
169
170 db=ComplexF32(-0.015f0,0.004f0)
171 dbn,ubn=MHDKernels._photon_exact_acoustic(
172 db,ub,0f0,Float32(2pi),1f0,3f0)
173 @test dbn == db
174 @test ubn == ub
175 end
176
177 @testset "exact exponential photon-moment propagator" begin
178 function matrix_reference(x,force,dt,G,nu,stress,k,cs2,cph2)
179 ik=im*k
180 operator=ComplexF64[0 -ik 0 0 0;
181 -ik*cs2 -G 0 G 1;
182 0 0 0 -(4/3)*ik 0;
183 0 nu -ik*cph2/4 -(nu+stress) 0;
184 0 0 0 0 0]
185 (exp(dt*operator)*ComplexF64[x...;force])[1:4]
186 end
187 G=37.0; nu=11.0; stress=2.5; k=8.0; cs2=0.7; cph2=19.0
188 force=0.13-0.07im
189 x0=ComplexF64[0.02+0.03im,-0.04+0.01im,0.01-0.02im,0.03+0.04im]
190 dt=0.08
191 ref=matrix_reference(x0,force,dt,G,nu,stress,k,cs2,cph2)
192 exact=collect(MHDKernels._photon_split_longitudinal(
193 Tuple(x0)...,force,dt,G,nu,stress,k,cs2,cph2))
194 @test exact ≈ ref rtol=3e-12 atol=3e-13
195
196 x32=ComplexF32.(x0); force32=ComplexF32(force)
197 exact32=collect(MHDKernels._photon_split_longitudinal(
198 Tuple(x32)...,force32,Float32(dt),Float32(G),Float32(nu),
199 Float32(stress),Float32(k),Float32(cs2),Float32(cph2)))
200 @test exact32 ≈ ComplexF32.(ref) rtol=2f-5 atol=2f-6
201
202 N=16
203 s=allocate_state(backend(:cpu),Float32,(N,N,N);dx=1/N)
204 ws=allocate_radiation_workspace(s)
205 rc=RadiationClosure(inertia_ratio=Float32(nu/G),mean_free_path=0.03f0,
206 photon_speed=sqrt(Float32(cph2)),gas_sound_speed=sqrt(Float32(cs2)),
207 bridge_q2=5f0,longitudinal_viscosity=4f0/3f0,
208 nyquist_guard_fraction=1f0,
209 response_model=:moments_exponential)
210 kfac=Float32(2pi)
211 mode2=Int32(3)
212 kcache=kfac*sqrt(Float32(mode2))
213 _,sl=MHDKernels._photon_stress_rates(kcache^2,Float32(nu),rc)
214 MHDKernels._radiation_ensure_longitudinal_propagator!(ws,s,
215 Float32(dt),Float32(G),rc,kfac)
216 cached=collect(MHDKernels._photon_apply_longitudinal_propagator(
217 ws.longitudinal_propagator,mode2,Tuple(x32)...,force32,
218 rc.photon_speed))
219 direct=collect(MHDKernels._photon_split_longitudinal(
220 Tuple(x32)...,force32,Float32(dt),Float32(G),Float32(nu),sl,
221 kcache,Float32(cs2),Float32(cph2)))
222 @test cached ≈ direct rtol=3f-6 atol=3f-7
223
224 # The affine drag solve must retain both the zero/slow mode and forcing
225 # when Gamma*dt is far outside the range of an explicit update.
226 ub,ug=MHDKernels._photon_exact_drag_affine(1f0+0f0im,-2f0+0f0im,
227 0.25f0+0f0im,1f0,1f8,2f7,0f0)
228 conserved=(1f8*ug+2f7*ub)/(1.2f8)
229 expected=(1f8*(-2f0)+2f7*1f0)/(1.2f8)+0.25f0*(2f7/1.2f8)
230 @test isfinite(real(ub)) && isfinite(real(ug))
231 @test real(ub) ≈ expected rtol=3f-6
232 @test real(ug) ≈ expected rtol=3f-6
233 @test conserved ≈ expected rtol=3f-6
234 end
235
236 @testset "RFFT self-conjugate planes remain Hermitian" begin
237 N=16
238 s=allocate_state(backend(:cpu),Float32,(N,N,N);dx=1/N)
239 ws=allocate_radiation_workspace(s)
240 real=randn(Float32,N,N,N,4)
241 copyto!(ws.real_batch,real)
242 MHDKernels._radiation_forward_batch!(ws.hat_batch,ws.real_batch)
243 rc=MHDKernels.RadiationClosure{Float32}(2f-4,1.5f-4,3.7f4,58f0,
244 5f0,4f0/3f0,15f0/16f0)
245 MHDKernels._radiation_homogeneous_hat_k!(s.be)(
246 ws.velocity_hat...,ws.density_hat,Int32(N),Float32(2pi),8.4f-8,
247 1.3f12,1f0,rc,Val(true);ndrange=length(ws.density_hat))
248 KernelAbstractions.synchronize(s.be)
249 hat=Array(ws.hat_batch)
250 for ii in (1,N÷2+1)
251 mismatch=0f0
252 scale=0f0
253 @inbounds for d in 1:4,k in 1:N,j in 1:N
254 jp=mod(-(j-1),N)+1
255 kp=mod(-(k-1),N)+1
256 mismatch=max(mismatch,abs(hat[ii,j,k,d]-conj(hat[ii,jp,kp,d])))
257 scale=max(scale,abs(hat[ii,j,k,d]))
258 end
259 @test mismatch <= 2f-5*max(scale,1f0)
260 end
261 end
262
263 @testset "pseudospectral correction has a Nyquist guard" begin
264 N=16
265 s=allocate_state(backend(:cpu),Float32,(N,N,N);dx=1/N)
266 ws=allocate_radiation_workspace(s)
267 copyto!(ws.real_batch,randn(Float32,N,N,N,4))
268 MHDKernels._radiation_forward_batch!(ws.hat_batch,ws.real_batch)
269 rc=MHDKernels.RadiationClosure{Float32}(2f-4,1.5f-4,3.7f4,58f0,
270 5f0,4f0/3f0)
271 MHDKernels._radiation_homogeneous_hat_k!(s.be)(
272 ws.velocity_hat...,ws.density_hat,Int32(N),Float32(2pi),8.4f-8,
273 1.3f12,1f0,rc,Val(true);ndrange=length(ws.density_hat))
274 KernelAbstractions.synchronize(s.be)
275 hat=Array(ws.hat_batch)
276 @inbounds for k in 1:N,j in 1:N,i in 1:N÷2+1
277 ii=Int32(i-1); jj=Int32(j-1); kk=Int32(k-1)
278 MHDKernels._radiation_mode_dealiased(
279 ii,jj,kk,Int32(N),rc.nyquist_guard_fraction) && continue
280 @test all(iszero,hat[i,j,k,:])
281 end
282 end
283
284 @testset "transverse force uses exact damped response" begin
285 N=16; s=allocate_state(backend(:cpu),Float32,(N,N,N);dx=1/N)
286 ws=allocate_radiation_workspace(s)
287 f=ntuple(_->zeros(Float32,N,N,N),3)
288 @inbounds for k in 1:N,j in 1:N,i in 1:N
289 f[2][i,j,k]=sinpi(2f0*Float32(i-1)/Float32(N))
290 end
291 MHDKernels._radiation_forward_triplet!(ws,f)
292 rc=MHDKernels.RadiationClosure{Float32}(0.2,0.01,20,1,5,4/3)
293 dt=0.03f0; gamma=7f0; k=Float32(2pi)
294 MHDKernels._radiation_source_hat_k!(s.be)(ws.velocity_hat...,ws.density_hat,
295 Int32(N),k,dt,gamma,1f0,rc,Val(false);
296 ndrange=length(ws.velocity_hat[1]))
297 MHDKernels._radiation_inverse_triplet!(f,ws)
298 mode=radiation_mode_coefficients(k,gamma,rc)
299 expected=mode.inertia_response*dt*MHDKernels._drag_phi1(mode.gamma_transverse*dt)-
300 dt*MHDKernels._drag_phi1(gamma*dt)
301 measured=2f0*sum(f[2][i,1,1]*sinpi(2f0*Float32(i-1)/Float32(N)) for i in 1:N)/N
302 @test measured ≈ expected rtol=3f-5 atol=2f-7
303 @test maximum(abs,f[1]) < 2f-7
304 @test maximum(abs,f[3]) < 2f-7
305 end
306
307 @testset "radiation-loaded longitudinal EOS response" begin
308 N=16; s=allocate_state(backend(:cpu),Float32,(N,N,N);dx=1/N)
309 amp=1f-3
310 fill!(s.U[1],1f0); fill!(s.U[5],2f0)
311 @inbounds for k in 1:N,j in 1:N,i in 1:N
312 s.U[1][linidx(N,N,i,j,k)]=1f0+amp*cospi(2f0*Float32(i-1)/Float32(N))
313 end
314 ws=allocate_radiation_workspace(s)
315 rc=MHDKernels.RadiationClosure{Float32}(0.2,0.01,20,1,5,4/3)
316 dt=0.004f0; gamma=7f0
317 MHDKernels._radiation_transform_velocity_density!(ws,s,s.U,dt,gamma,rc)
318 vx=ws.correction[1]
319 measured=2f0*sum(vx[i,1,1]*sinpi(2f0*Float32(i-1)/Float32(N)) for i in 1:N)/N
320 @test measured > 0
321 @test maximum(abs,ws.correction[2]) < 2f-7
322 @test maximum(abs,ws.correction[3]) < 2f-7
323 end
324
325 @testset "full step follows exact longitudinal mode above radiation CFL" begin
326 N=64; gamma=5f0/3f0; amp=2f-4; p0=1f0/gamma
327 s=allocate_state(backend(:cpu),Float32,(N,N,N);dx=1/N,
328 gamma=gamma,riemann=:hll)
329 @inbounds for k in 1:N,j in 1:N,i in 1:N
330 phase=2f0*Float32(pi)*Float32(i-1)/Float32(N)
331 rho=1f0+amp*cos(phase)
332 p=p0*rho^gamma
333 c=linidx(N,N,i,j,k)
334 s.U[1][c]=rho; s.U[2][c]=0; s.U[3][c]=0; s.U[4][c]=0
335 s.U[5][c]=p/(gamma-1f0)
336 s.U[6][c]=0; s.U[7][c]=0; s.U[8][c]=0; s.U[9][c]=0
337 end
338 ws=allocate_radiation_workspace(s)
339 rc=MHDKernels.RadiationClosure{Float32}(0.2,0.001,20,1,5,4/3)
340 dt=0.004f0; drag_gamma=7f0; k1=Float32(2pi)
341 step_radiation_godunov!(s,ws,dt;drag_impulse=drag_gamma*dt,
342 closure=rc,ch=1f0,integrator=:ref)
343 rho=reshape(Array(s.U[1]),N,N,N)
344 delta=Array(ws.density_perturbation)
345 vx=reshape(Array(s.U[2]),N,N,N)./rho
346 rhoamp=2f0*sum(delta[i,1,1]*cospi(2f0*Float32(i-1)/N) for i in 1:N)/N
347 vxamp=2f0*sum(vx[i,1,1]*sinpi(2f0*Float32(i-1)/N) for i in 1:N)/N
348 mode=radiation_mode_coefficients(k1,drag_gamma,rc)
349 c2=1f0+mode.radiation_sound_speed2
350 B,S=MHDKernels._radiation_acoustic_coefficients(
351 mode.gamma_longitudinal,c2,k1*k1,dt)
352 exact_rho=amp*(B+mode.gamma_longitudinal*S)
353 exact_v=amp*k1*c2*S
354 @info "radiation longitudinal mode" rhoamp exact_rho vxamp exact_v
355 @test rhoamp ≈ exact_rho rtol=3f-2 atol=2f-6
356 @test vxamp ≈ exact_v rtol=4f-2 atol=2f-6
357 @test abs(Float64(sum(delta)) / N^3) < 1e-6
358 @test minimum(rho) > 0
359 end
360
361 @testset "Lorentz-forced mode follows exact radiation response" begin
362 N=64; gamma=5f0/3f0; p0=1f0/gamma; bg=0.5f0; db=0.05f0
363 s=allocate_state(backend(:cpu),Float32,(N,N,N);dx=1/N,
364 gamma=gamma,riemann=:hll)
365 @inbounds for k in 1:N,j in 1:N,i in 1:N
366 phase=2f0*Float32(pi)*Float32(i-1)/Float32(N)
367 by=bg+db*cos(phase); c=linidx(N,N,i,j,k)
368 s.U[1][c]=1; s.U[2][c]=0; s.U[3][c]=0; s.U[4][c]=0
369 s.U[5][c]=p0/(gamma-1f0)+0.5f0*by*by
370 s.U[6][c]=0; s.U[7][c]=by; s.U[8][c]=0; s.U[9][c]=0
371 end
372 ws=allocate_radiation_workspace(s)
373 rc=MHDKernels.RadiationClosure{Float32}(0.2,0.001,20,1,5,4/3)
374 dt=0.004f0; drag_gamma=7f0; k1=Float32(2pi)
375 step_radiation_godunov!(s,ws,dt;drag_impulse=drag_gamma*dt,
376 closure=rc,ch=1f0,integrator=:ref)
377 rho=reshape(Array(s.U[1]),N,N,N); vx=reshape(Array(s.U[2]),N,N,N)./rho
378 rhoamp=2f0*sum((rho[i,1,1]-1f0)*cospi(2f0*Float32(i-1)/N) for i in 1:N)/N
379 vxamp=2f0*sum(vx[i,1,1]*sinpi(2f0*Float32(i-1)/N) for i in 1:N)/N
380 mode=radiation_mode_coefficients(k1,drag_gamma,rc)
381 c2=1f0+mode.radiation_sound_speed2
382 B,S=MHDKernels._radiation_acoustic_coefficients(
383 mode.gamma_longitudinal,c2,k1*k1,dt)
384 J=MHDKernels._radiation_integral_S(
385 B,S,mode.gamma_longitudinal,c2*k1*k1,dt)
386 accel=bg*db*sin(k1/N)*N
387 exact_rho=-mode.inertia_response*k1*J*accel
388 exact_v=mode.inertia_response*S*accel
389 @info "radiation Lorentz-forced mode" rhoamp exact_rho vxamp exact_v
390 @test rhoamp ≈ exact_rho rtol=8f-2 atol=2f-7
391 @test vxamp ≈ exact_v rtol=8f-2 atol=2f-6
392 @test abs(Float64(sum(rho)) / N^3 - 1) < 1e-6
393 @test minimum(rho) > 0
394 end
395
396 @testset "free-streaming limit recovers local Godunov drag" begin
397 N=8
398 a=allocate_state(backend(:cpu),Float32,(N,N,N);dx=1/N,riemann=:hll)
399 b=allocate_state(backend(:cpu),Float32,(N,N,N);dx=1/N,riemann=:hll)
400 init_turb_field!(a)
401 for v in 1:9 copyto!(b.U[v],a.U[v]) end
402 dt,smax=compute_dt(a;cfl=0.03)
403 step_drag_godunov!(a,dt;drag_impulse=0.7f0,ch=smax,integrator=:ref)
404 rc=RadiationClosure(inertia_ratio=0.2f0,mean_free_path=1f6,
405 photon_speed=20f0,gas_sound_speed=1f0)
406 ws=allocate_radiation_workspace(b)
407 step_radiation_godunov!(b,ws,dt;drag_impulse=0.7f0,closure=rc,
408 ch=smax,integrator=:ref)
409 ha=fields_to_host(a); hb=fields_to_host(b)
410 @test maximum(maximum(abs,hb[v].-ha[v]) for v in 1:9) < 3f-6
411 end
412
413
414 @testset "device energy ledger classifies the coupled step" begin
415 N=8
416 s=allocate_state(backend(:cpu),Float32,(N,N,N);dx=1/N,riemann=:hll)
417 init_turb_field!(s)
418 ws=allocate_radiation_workspace(s)
419 before=radiation_energy_components(ws,s)
420 hbefore=fields_to_host(s)
421 @test before.total ≈ before.kinetic+before.internal+before.magnetic rtol=2e-6
422 rc=RadiationClosure(inertia_ratio=0.2f0,mean_free_path=0.01f0,
423 photon_speed=20f0,gas_sound_speed=1f0,
424 response_model=:angular)
425 ledger=RadiationEnergyLedger()
426 dt,smax=compute_dt(s;cfl=0.02)
427 step_radiation_godunov!(s,ws,dt;drag_impulse=0.3f0,closure=rc,
428 ch=smax,integrator=:ref,energy_ledger=ledger)
429 after=radiation_energy_components(ws,s)
430 hafter=fields_to_host(s)
431 function host_components(h)
432 rho=Float64.(h[1]); mx=Float64.(h[2]); my=Float64.(h[3]); mz=Float64.(h[4])
433 bx=Float64.(h[6]); by=Float64.(h[7]); bz=Float64.(h[8])
434 ek=sum(0.5 .* (mx.^2 .+ my.^2 .+ mz.^2) ./ rho)
435 eb=sum(0.5 .* (bx.^2 .+ by.^2 .+ bz.^2))
436 total=sum(Float64.(h[5]))
437 (total=total,kinetic=ek,internal=total-ek-eb,magnetic=eb)
438 end
439 hb=host_components(hbefore); ha=host_components(hafter)
440 @test ledger.samples == 1
441 @test ledger.total_change ≈ ha.total-hb.total rtol=2e-5 atol=2e-6
442 @test ledger.kinetic_change ≈ ha.kinetic-hb.kinetic rtol=2e-5 atol=2e-6
443 @test ledger.internal_change ≈ ha.internal-hb.internal rtol=1e-4 atol=2e-6
444 @test ledger.magnetic_change ≈ ha.magnetic-hb.magnetic rtol=2e-5 atol=2e-6
445 @test all(isfinite,(ledger.total_change,ledger.kinetic_change,
446 ledger.internal_change,ledger.magnetic_change))
447 @test ledger.component_residual_abs >= abs(ledger.component_residual)
448 @test ledger.component_magnitude_abs >= ledger.component_residual_abs
449 component_scale=abs(ledger.kinetic_change)+abs(ledger.internal_change)+
450 abs(ledger.magnetic_change)
451 @test abs(ledger.component_residual)/component_scale < 5e-4
452 @test ledger.component_residual_abs/ledger.component_magnitude_abs < 1e-6
453 end
454
455 @testset "evolved photon moments close momentum and angular stress" begin
456 R=0.2f0; G=7f0; nu=R*G; dt=0.4f0
457 ub0=ComplexF32(0.7,-0.2); ug0=ComplexF32(-0.1,0.3)
458 ub1,ug1=MHDKernels._photon_sdirk_transverse(ub0,ug0,0f0,dt,G,nu,0f0)
459 @test ub1+ug1/R ≈ ub0+ug0/R rtol=3f-6
460
461 rc=RadiationClosure(inertia_ratio=R,mean_free_path=1f-4,
462 photon_speed=20f0,gas_sound_speed=1f0,
463 bridge_q2=5f0,longitudinal_viscosity=4f0/3f0,
464 nyquist_guard_fraction=1f0,
465 response_model=:moments)
466 k=3f0
467 st,sl=MHDKernels._photon_stress_rates(k*k,nu,rc)
468 @test st ≈ nu*(k*rc.mean_free_path)^2/5 rtol=5f-3
469 @test sl/st ≈ 4f0/3f0 rtol=2f-4
470 @test radiation_mode_coefficients(k,G,rc).response_model === :moments
471
472 streaming=RadiationClosure(inertia_ratio=R,mean_free_path=1f4,
473 photon_speed=20f0,gas_sound_speed=1f0,
474 bridge_q2=5f0,longitudinal_viscosity=4f0/3f0,
475 nyquist_guard_fraction=1f0,
476 response_model=:moments)
477 stfree,_=MHDKernels._photon_stress_rates(k*k,nu,streaming)
478 @test stfree > 0
479 @test isfinite(stfree)
480
481 # Eq. 57: an adiabatic tight photon-baryon mode propagates with
482 # c/sqrt(3(1+R)). Integrate the closed moments directly so this gate is
483 # independent of finite-volume spatial error.
484 Gtight=1f4; nutight=R*Gtight; lambda=20f0/nutight
485 tight=RadiationClosure(inertia_ratio=R,mean_free_path=lambda,
486 photon_speed=20f0,gas_sound_speed=0f0,bridge_q2=5f0,
487 longitudinal_viscosity=4f0/3f0,nyquist_guard_fraction=1f0,
488 response_model=:moments)
489 _,sltight=MHDKernels._photon_stress_rates(1f0,nutight,tight)
490 db=ComplexF32(1,0); ub=ComplexF32(0); dg=ComplexF32(4f0/3f0); ug=ComplexF32(0)
491 dti=1f-4
492 for _ in 1:1000
493 db,ub,dg,ug=MHDKernels._photon_sdirk_longitudinal(
494 db,ub,dg,ug,0f0,dti,Gtight,nutight,sltight,1f0,0f0,400f0)
495 end
496 ceq=20f0/sqrt(3f0*(1f0+R))
497 expected=cos(ceq*0.1f0)*exp(-0.5f0*sltight/(1f0+R)*0.1f0)
498 @test real(db) ≈ expected rtol=8f-3
499 @test abs(imag(db)) < 2f-5
500
501 N=16; sn=allocate_state(backend(:cpu),Float32,(N,N,N);dx=1/N)
502 wn=allocate_radiation_workspace(sn); pn=allocate_photon_moment_state(sn)
503 fill!(wn.hat_batch,0); fill!(pn.hat_batch,0)
504 pn.velocity_hat[1][N÷2+1,1,1]=ComplexF32(1)
505 MHDKernels._photon_moment_homogeneous_hat_k!(sn.be)(
506 wn.velocity_hat...,wn.density_hat,wn.longitudinal_propagator,
507 pn.velocity_hat...,pn.delta_hat,
508 Int32(N),Float32(2pi),0.1f0,7f0,1f0,tight,Val(true),Val(true);
509 ndrange=length(wn.density_hat))
510 KernelAbstractions.synchronize(sn.be)
511 @test pn.velocity_hat[1][N÷2+1,1,1]==ComplexF32(1)
512 @test all(iszero,wn.hat_batch)
513 end
514
515 @testset "photon moment state advances with full Godunov MHD" begin
516 N=16; gamma=5f0/3f0; amp=1f-4
517 s=allocate_state(backend(:cpu),Float32,(N,N,N);dx=1/N,
518 gamma=gamma,riemann=:hll)
519 @inbounds for k in 1:N,j in 1:N,i in 1:N
520 phase=2f0*Float32(pi)*Float32(i-1)/Float32(N)
521 rho=1f0+amp*cos(phase); vx=amp*sin(phase); c=linidx(N,N,i,j,k)
522 s.U[1][c]=rho; s.U[2][c]=rho*vx; s.U[3][c]=0; s.U[4][c]=0
523 s.U[5][c]=1f0/(gamma-1f0)+0.5f0*rho*vx*vx
524 s.U[6][c]=0; s.U[7][c]=0; s.U[8][c]=0; s.U[9][c]=0
525 end
526 ws=allocate_radiation_workspace(s)
527 @test length(ws.packed_psi_correction)==9*N^3
528 @test all(r -> all(iszero,r),ws.magnetic_residual)
529 photons=allocate_photon_moment_state(s)
530 initialize_photon_moments!(photons,ws,s)
531 rc=RadiationClosure(inertia_ratio=0.2f0,mean_free_path=0.001f0,
532 photon_speed=20f0,gas_sound_speed=1f0,
533 response_model=:moments)
534 e0=photon_moment_energy(photons,ws,s,rc)
535 p0=copy(Array(photons.hat_batch))
536 ledger=RadiationEnergyLedger()
537 step_radiation_godunov!(s,ws,2f-4;drag_impulse=7f0*2f-4,
538 closure=rc,ch=1f0,integrator=:ref,
539 photon_state=photons,energy_ledger=ledger)
540 e1=photon_moment_energy(photons,ws,s,rc)
541 h=fields_to_host(s)
542 @test all(all(isfinite,x) for x in h)
543 @test minimum(h[1])>0
544 @test isfinite(e0) && e0>0
545 @test isfinite(e1) && e1>0
546 @test maximum(abs,Array(photons.hat_batch).-p0)>0
547 @test ledger.samples==1
548 @test ledger.photon_change ≈ e1-e0 rtol=2f-5 atol=1f-8
549 @test ledger.resolved_coupled_change ≈
550 ledger.total_change+ledger.photon_change rtol=2f-6 atol=1f-9
551 @test ledger.component_residual_abs >= abs(ledger.component_residual)
552 @test ledger.component_magnitude_abs >= ledger.component_residual_abs
553
554 missing=allocate_state(backend(:cpu),Float32,(8,8,8);dx=1/8)
555 mw=allocate_radiation_workspace(missing)
556 @test_throws ErrorException step_radiation_godunov!(missing,mw,1f-4;
557 drag_impulse=1f-3,closure=rc,ch=1f0,integrator=:ref)
558 end
559
560 @testset "gated free-streaming moment path recovers spectral update" begin
561 N=8
562 full=allocate_state(backend(:cpu),Float32,(N,N,N);dx=1/N,riemann=:hll)
563 fast=allocate_state(backend(:cpu),Float32,(N,N,N);dx=1/N,riemann=:hll)
564 init_turb_field!(full)
565 for v in 1:9
566 copyto!(fast.U[v],full.U[v])
567 end
568 wf=allocate_radiation_workspace(full)
569 wa=allocate_radiation_workspace(fast)
570 pf=allocate_photon_moment_state(full)
571 pa=allocate_photon_moment_state(fast)
572 initialize_photon_moments!(pf,wf,full)
573 copyto!(pa.hat_batch,pf.hat_batch)
574 rc=RadiationClosure(inertia_ratio=0.3f0,mean_free_path=1f4,
575 photon_speed=20f0,gas_sound_speed=1f0,
576 response_model=:moments)
577 dt,smax=compute_dt(full;cfl=0.01)
578 impulse=0.2f0
579 bound=radiation_free_streaming_bound(2f0*Float32(pi),impulse/dt,rc)
580 @test bound < 1f-3
581 tf=RadiationStepTiming()
582 ta=RadiationStepTiming()
583 la=RadiationEnergyLedger()
584 step_radiation_godunov!(full,wf,dt;drag_impulse=impulse,closure=rc,
585 ch=smax,integrator=:ref,photon_state=pf,
586 timing=tf)
587 step_radiation_godunov!(fast,wa,dt;drag_impulse=impulse,closure=rc,
588 ch=smax,integrator=:ref,photon_state=pa,
589 timing=ta,energy_ledger=la,
590 free_streaming_tolerance=1f-3)
591 hf=fields_to_host(full)
592 ha=fields_to_host(fast)
593 gas_error=maximum(maximum(abs,ha[v].-hf[v]) for v in 1:9)
594 photon_scale=max(maximum(abs,Array(pf.hat_batch)),1f-20)
595 photon_error=maximum(abs,Array(pa.hat_batch).-Array(pf.hat_batch))/
596 photon_scale
597 @test gas_error < 4f-5
598 @test photon_error < 2f-3
599 @test tf.spectral_steps==1 && tf.fast_steps==0
600 @test ta.fast_steps==1 && ta.spectral_steps==0
601 @test la.samples==1
602 @test la.component_residual_abs <
603 2f-5*max(abs(la.total_change),1f-8)
604
605 coupled=RadiationClosure(inertia_ratio=0.3f0,mean_free_path=0.01f0,
606 photon_speed=20f0,gas_sound_speed=1f0,
607 response_model=:moments)
608 @test radiation_free_streaming_bound(2f0*Float32(pi),
609 impulse/dt,coupled)>1f-3
610 end
611
612 @testset "checkpoint fallback rejects without accepting the trial" begin
613 N=8
614 s=allocate_state(backend(:cpu),Float32,(N,N,N);dx=1/N,riemann=:hll)
615 init_turb_field!(s)
616 ws=allocate_radiation_workspace(s)
617 before=fields_to_host(s)
618 for v in 1:9
619 copyto!(s.scratch[v],s.U[v])
620 end
621 fill!(s.U[1],-1f-4)
622 err=try
623 MHDKernels._checkpoint_fallback_rejection!(s,ws,true;
624 dt=1f-7,debug_context="unit-test",trial_rho_min=-1f-4,
625 trial_wave_speed=0f0,wave_limit=1f0)
626 nothing
627 catch caught
628 caught
629 end
630 @test err isa ErrorException
631 @test occursin("RADIATION_CHECKPOINT_CFL_FALLBACK",sprint(showerror,err))
632 after=fields_to_host(s)
633 @test all(after[v]==before[v] for v in 1:9)
634 expected_delta=reshape(before[1],s.dims).-ws.rho_mean
635 @test Array(ws.density_perturbation)==expected_delta
636 end
637 end