Snapshot from this documentation build. SHA-256: 33711f727aa1982d6933d34cb2aa44adc84851b115706726c8e1b98755d2e071
1 using Test
2
3 @testset "nonlocal Lyman-alpha recombination transport" begin
4 N = 16
5 backend_name = Symbol(get(ENV, "MHD_TEST_BACKEND", "cpu"))
6 be = backend(backend_name)
7 s = allocate_state(be, Float32, (N, N, N); dx=1 / N)
8 fh = 0.76f0
9 xe_mean = 0.1f0
10 density_amplitude = 0.01f0
11 ionization_amplitude = 0.005f0
12 velocity_amplitude = 0.02f0
13 rho_host = zeros(Float32, N^3)
14 mx_host = zeros(Float32, N^3)
15 hii_host = zeros(Float32, N^3)
16 @inbounds for k in 1:N, j in 1:N, i in 1:N
17 c = i + N * ((j - 1) + N * (k - 1))
18 phase = 2f0 * Float32(pi) * (i - 1) / N
19 rho = 1f0 + density_amplitude * cos(phase)
20 rho_host[c] = rho
21 mx_host[c] = rho * velocity_amplitude * sin(phase)
22 hii_host[c] = fh * rho * (xe_mean + ionization_amplitude * cos(phase))
23 end
24 copyto!(s.U[1], rho_host)
25 copyto!(s.U[2], mx_host)
26 HII = to_device(be, hii_host, Float32)
27 H2I = device_zeros(be, Float32, (N^3,))
28 workspace = allocate_radiation_workspace(s)
29 redshifts = expm1.(range(log1p(500.0), log1p(2000.0), length=3))
30 wavenumbers = exp.(range(log(1.0), log(100.0), length=3))
31 density_rate = fill(-2.0e-12, 3, 3)
32 ionization_rate = fill(-3.0e-12, 3, 3)
33 velocity_rate = fill(4.0e-12, 3, 3)
34 table = NonlocalRecombinationTable(be, redshifts, wavenumbers,
35 density_rate, ionization_rate, velocity_rate; precision=Float32)
36 dt_seconds = 1.0e8
37 applied = apply_nonlocal_recombination!(HII, H2I, s, workspace, table;
38 redshift=redshifts[2], xe_mean=xe_mean, dt_seconds=dt_seconds,
39 fundamental_k_mpc=10.0, theta_over_aH_conversion=1.0,
40 hydrogen_mass_fraction=fh)
41 @test applied
42 q = -3.0e-12 * dt_seconds
43 phi = expm1(q) / -3.0e-12
44 theta_amplitude = velocity_amplitude * sinpi(2 / N) / (1 / N)
45 expected_amplitude = exp(q) * ionization_amplitude + phi * (
46 -2.0e-12 * density_amplitude + 4.0e-12 * theta_amplitude)
47 xe = to_host(HII) ./ (fh .* rho_host)
48 expected = Float32[
49 xe_mean + expected_amplitude * cospi(2f0 * (i - 1) / N)
50 for k in 1:N for j in 1:N for i in 1:N
51 ]
52 @test maximum(abs, xe .- expected) < 8f-7
53 @test abs(sum(xe) / length(xe) - xe_mean) < 2f-7
54
55 # The high-ionization representation stores neutral H in the same slot.
56 # It must yield the identical physical electron response without forming
57 # 1 - x_HI in Float32 until the diagnostic conversion below.
58 hi_host = fh .* rho_host .- hii_host
59 HI = to_device(be, hi_host, Float32)
60 workspace_neutral = allocate_radiation_workspace(s)
61 @test apply_nonlocal_recombination!(HI, H2I, s, workspace_neutral, table;
62 redshift=redshifts[2], xe_mean=xe_mean, dt_seconds=dt_seconds,
63 fundamental_k_mpc=10.0, theta_over_aH_conversion=1.0,
64 hydrogen_mass_fraction=fh, neutral_hydrogen_storage=true)
65 xe_neutral = 1f0 .- to_host(HI) ./ (fh .* rho_host)
66 @test maximum(abs, xe_neutral .- xe) < 8f-7
67 @test maximum(abs, xe_neutral .- expected) < 8f-7
68
69 # Centered storage carries x_HII-<x_HII> directly. Verify that a mode much
70 # smaller than an ULP of a nearly ionized background survives the FFT
71 # response instead of disappearing when added to that background.
72 centered_reference = Float32(1 - 2.0^-20)
73 centered_amplitude = 3f-8
74 centered_host = Float32[
75 centered_amplitude * cospi(2f0 * (i - 1) / N)
76 for k in 1:N for j in 1:N for i in 1:N
77 ]
78 HII_centered = to_device(be, centered_host, Float32)
79 negligible_density_rate = fill(-1.0e-30, 3, 3)
80 negligible_velocity_rate = fill(1.0e-30, 3, 3)
81 centered_table = NonlocalRecombinationTable(
82 be, redshifts, wavenumbers, negligible_density_rate, ionization_rate,
83 negligible_velocity_rate;
84 precision=Float32,
85 )
86 workspace_centered = allocate_radiation_workspace(s)
87 @test apply_nonlocal_recombination!(
88 HII_centered, H2I, s, workspace_centered, centered_table;
89 redshift=redshifts[2], xe_mean=centered_reference,
90 dt_seconds=dt_seconds, fundamental_k_mpc=10.0,
91 theta_over_aH_conversion=0.0, hydrogen_mass_fraction=fh,
92 centered_hydrogen_storage=true,
93 hydrogen_reference=centered_reference,
94 hydrogen_complement_reference=Float32(1 - Float64(centered_reference)),
95 )
96 centered_result = to_host(HII_centered)
97 centered_expected = exp(Float32(q)) .* centered_host
98 @test maximum(abs, centered_result .- centered_expected) < 2f-10
99 @test maximum(abs, centered_result) > 2f-8
100
101 # A production table also carries the k=0 Jacobian already advanced by
102 # the nonlinear local network. The Fourier correction must replace that
103 # finite-step local map with the total map, rather than Lie-splitting two
104 # stiff exponentials.
105 local_density_rate = -2.0e-12
106 local_ionization_rate = -4.0e-12
107 differential_density_rate = fill(-3.0e-12, 3, 3)
108 differential_ionization_rate = fill(-6.0e-12, 3, 3)
109 local_response = zeros(3, 4)
110 local_response[:, 1] .= local_density_rate
111 local_response[:, 2] .= local_ionization_rate
112 composed_table = NonlocalRecombinationTable(
113 be, redshifts, wavenumbers, differential_density_rate,
114 differential_ionization_rate, negligible_velocity_rate;
115 local_response_rate=local_response, precision=Float32,
116 )
117 composed_dt = 1.0e11
118 local_phi = expm1(local_ionization_rate * composed_dt) /
119 local_ionization_rate
120 local_amplitude = local_phi * local_density_rate * density_amplitude
121 composed_initial = Float32[
122 local_amplitude * cospi(2f0 * (i - 1) / N)
123 for k in 1:N for j in 1:N for i in 1:N
124 ]
125 HII_composed = to_device(be, composed_initial, Float32)
126 workspace_composed = allocate_radiation_workspace(s)
127 composed_reference = 0.5f0
128 @test apply_nonlocal_recombination!(
129 HII_composed, H2I, s, workspace_composed, composed_table;
130 redshift=redshifts[2], xe_mean=composed_reference,
131 dt_seconds=composed_dt, fundamental_k_mpc=10.0,
132 theta_over_aH_conversion=0.0, hydrogen_mass_fraction=fh,
133 centered_hydrogen_storage=true,
134 hydrogen_reference=composed_reference,
135 hydrogen_complement_reference=1f0 - composed_reference,
136 )
137 total_ionization_rate = local_ionization_rate - 6.0e-12
138 total_density_rate = local_density_rate - 3.0e-12
139 total_phi = expm1(total_ionization_rate * composed_dt) /
140 total_ionization_rate
141 expected_composed_amplitude = total_phi * total_density_rate * density_amplitude
142 expected_composed = Float32[
143 expected_composed_amplitude * cospi(2f0 * (i - 1) / N)
144 for k in 1:N for j in 1:N for i in 1:N
145 ]
146 @test maximum(abs, to_host(HII_composed) .- expected_composed) < 8f-7
147
148 before = to_host(HII)
149 @test !apply_nonlocal_recombination!(HII, H2I, s, workspace, table;
150 redshift=3000.0, xe_mean=xe_mean, dt_seconds=dt_seconds,
151 fundamental_k_mpc=10.0, theta_over_aH_conversion=0.0,
152 hydrogen_mass_fraction=fh)
153 @test to_host(HII) == before
154
155 # Below the first tabulated k, the differential response must vanish as k^2.
156 hii_lowk_host = similar(hii_host)
157 @inbounds for k in 1:N, j in 1:N, i in 1:N
158 c = i + N * ((j - 1) + N * (k - 1))
159 phase = 2f0 * Float32(pi) * (i - 1) / N
160 hii_lowk_host[c] = fh * rho_host[c] *
161 (xe_mean + ionization_amplitude * cos(phase))
162 end
163 HII_lowk = to_device(be, hii_lowk_host, Float32)
164 low_scale = (0.1 / first(wavenumbers))^2
165 @test apply_nonlocal_recombination!(HII_lowk, H2I, s, workspace, table;
166 redshift=redshifts[2], xe_mean=xe_mean, dt_seconds=dt_seconds,
167 fundamental_k_mpc=0.1, theta_over_aH_conversion=1.0,
168 hydrogen_mass_fraction=fh)
169 q_low = low_scale * -3.0e-12 * dt_seconds
170 phi_low = expm1(q_low) / (low_scale * -3.0e-12)
171 expected_low_amplitude = exp(q_low) * ionization_amplitude + phi_low * low_scale * (
172 -2.0e-12 * density_amplitude + 4.0e-12 * theta_amplitude)
173 xe_low = to_host(HII_lowk) ./ (fh .* rho_host)
174 expected_low = Float32[
175 xe_mean + expected_low_amplitude * cospi(2f0 * (i - 1) / N)
176 for k in 1:N for j in 1:N for i in 1:N
177 ]
178 @test maximum(abs, xe_low .- expected_low) < 8f-7
179
180 # The fourth real-field slot carries delta(T_b)/T_b through the same batched
181 # transform, so temperature response adds no FFT or full-grid workspace.
182 s_temperature = allocate_state(be, Float32, (N, N, N); dx=1 / N)
183 temperature_reference = 3000f0
184 temperature_amplitude = 0.003f0
185 velocity_unit2 = 1.0f12
186 temperature_rate = fill(5.0e-12, 3, 3)
187 table_temperature = NonlocalRecombinationTable(
188 be, redshifts, wavenumbers, density_rate, ionization_rate, velocity_rate;
189 temperature_rate=temperature_rate, precision=Float32,
190 )
191 rho_temperature = ones(Float32, N^3)
192 hii_temperature = fill(fh * xe_mean, N^3)
193 energy_temperature = similar(rho_temperature)
194 particles_per_mass_h = fh + fh * xe_mean + (1f0 - fh) / 4f0
195 @inbounds for k in 1:N, j in 1:N, i in 1:N
196 c = i + N * ((j - 1) + N * (k - 1))
197 phase = 2f0 * Float32(pi) * (i - 1) / N
198 temperature = temperature_reference *
199 (1f0 + temperature_amplitude * cos(phase))
200 specific_e = temperature * 1.380649f-16 * particles_per_mass_h /
201 ((Float32(s_temperature.γ) - 1f0) * 1.6735575f-24 *
202 velocity_unit2)
203 energy_temperature[c] = specific_e
204 end
205 copyto!(s_temperature.U[1], rho_temperature)
206 copyto!(s_temperature.U[5], energy_temperature)
207 HII_temperature = to_device(be, hii_temperature, Float32)
208 H2I_temperature = device_zeros(be, Float32, (N^3,))
209 workspace_temperature = allocate_radiation_workspace(s_temperature)
210 @test apply_nonlocal_recombination!(
211 HII_temperature, H2I_temperature, s_temperature, workspace_temperature,
212 table_temperature; redshift=redshifts[2], xe_mean=xe_mean,
213 dt_seconds=dt_seconds, fundamental_k_mpc=10.0,
214 theta_over_aH_conversion=0.0, hydrogen_mass_fraction=fh,
215 velocity_unit2=velocity_unit2,
216 temperature_reference=temperature_reference,
217 )
218 phi_temperature = expm1(q) / -3.0e-12
219 expected_temperature_amplitude = phi_temperature * 5.0e-12 *
220 temperature_amplitude
221 xe_temperature = to_host(HII_temperature) ./ fh
222 expected_temperature = Float32[
223 xe_mean + expected_temperature_amplitude * cospi(2f0 * (i - 1) / N)
224 for k in 1:N for j in 1:N for i in 1:N
225 ]
226 @test maximum(abs, xe_temperature .- expected_temperature) < 8f-7
227 @test table_temperature.has_temperature
228
229 mktemp() do path, io
230 println(io, join(("redshift", "k_mpc", "delta_density_rate_s",
231 "delta_ionization_rate_s",
232 "delta_velocity_divergence_rate_s"), '\t'))
233 for z in redshifts, k in wavenumbers
234 println(io, join((z, k, -2e-12, -3e-12, 4e-12), '\t'))
235 end
236 close(io)
237 loaded = load_nonlocal_recombination_table(path, be; precision=Float32)
238 @test loaded.nz == 3
239 @test loaded.nk == 3
240 @test isapprox(loaded.zlog_min, Float32(log1p(first(redshifts))); rtol=2f-6)
241 @test isapprox(loaded.logk_max, Float32(log(last(wavenumbers))); rtol=2f-6)
242 @test !loaded.has_temperature
243 end
244
245 mktemp() do path, io
246 println(io, join(("redshift", "k_mpc", "delta_density_rate_s",
247 "delta_ionization_rate_s",
248 "delta_velocity_divergence_rate_s",
249 "delta_temperature_rate_s",
250 "local_velocity_divergence_rate_s"), '\t'))
251 for z in redshifts, k in wavenumbers
252 println(io, join((z, k, -2e-12, -3e-12, 4e-12, 5e-12, -1e-12), '\t'))
253 end
254 close(io)
255 loaded = load_nonlocal_recombination_table(path, be; precision=Float32)
256 @test loaded.has_temperature
257 @test all(isapprox.(to_host(loaded.local_velocity_divergence_rate),
258 -1f-12; rtol=2f-6))
259 end
260
261 @test_throws ErrorException NonlocalRecombinationTable(
262 be, redshifts, wavenumbers, density_rate, -ionization_rate,
263 velocity_rate; precision=Float32)
264 end