Snapshot from this documentation build. SHA-256: f3923416fb55b5a65c8a47dbc405f8549cea44f616c9c9255a5c5eda614655dd
1 # Shared CPU and real-device checks; no device scalar indexing.
2 function streaming_checks(name, T)
3 M = PMFGLM.MHDKernels
4 K = KernelAbstractions
5 n = 8
6 be = M.backend(name)
7 s = M.allocate_state(be,T,(n,n,n); dx=1/n,gamma=5/3,recon=:plm,riemann=:hll)
8 r = T[1+0.2cospi(2*(i-0.5)/n) for i in 1:n,j in 1:n,k in 1:n]
9 for (d,field) in enumerate((r,zero.(r),zero.(r),zero.(r),fill(T(2),size(r)),
10 zero.(r),zero.(r),zero.(r),zero.(r)))
11 copyto!(s.U[d],vec(field))
12 end
13 v = (T(0.2),T(-0.3),T(0.4))
14 M.boost_gas_velocity!(s,v)
15 K.synchronize(be)
16 fields = M.fields_to_host(s)
17 @test maximum(abs,fields[5] .- (2 .+ vec(r).*sum(abs2,v)/2)) < 2e-6
18 for d in 1:3
19 @test maximum(abs,fields[d+1] ./ vec(r) .- v[d]) < 1e-6
20 end
21 positions = ntuple(_->M.device_zeros(be,T,(n^3,)),3)
22 velocities = ntuple(d->M.to_device(be,fill(v[d]/T(2),n^3)),3)
23 out = ntuple(_->M.device_zeros(be,T,(n^3,)),3)
24 # A common physical boost must cancel, despite different code units.
25 M.sample_baryon_dm_velocity!(out,s,positions,velocities;gas_unit=3,dm_unit=6)
26 K.synchronize(be)
27 @test maximum(maximum(abs,M.to_host(x)) for x in out) < 1e-6
28 # Probe an x-dependent gas velocity at wrapped positions, not particle labels.
29 speed = T[cospi(2*(i-0.5)/n) for i in 1:n,j in 1:n,k in 1:n]
30 copyto!(s.U[2],vec(r.*speed))
31 for d in 1:3
32 copyto!(positions[d],fill(T(d==1 ? -1.25/n : 0),n^3))
33 end
34 M.sample_baryon_dm_velocity!(out,s,positions,velocities;gas_unit=3,dm_unit=6)
35 K.synchronize(be)
36 expected = T[3*(0.25speed[mod1(i-2,n),j,k]+0.75speed[mod1(i-1,n),j,k])-3v[1]
37 for i in 1:n,j in 1:n,k in 1:n]
38 @test maximum(abs,M.to_host(out[1]).-vec(expected)) < 3e-6
39 # Absolute-position convention gives the same interpolation.
40 absolute = ntuple(d->T[(d==1 ? i-0.5-1.25 : d==2 ? j-0.5 : k-0.5)/n
41 for i in 1:n,j in 1:n,k in 1:n],3)
42 for d in 1:3
43 copyto!(positions[d],vec(absolute[d]))
44 end
45 M.sample_baryon_dm_velocity!(out,s,positions,velocities;
46 gas_unit=3,dm_unit=6,lattice_displacements=false)
47 K.synchronize(be)
48 @test maximum(abs,M.to_host(out[1]).-vec(expected)) < 3e-6
49
50 for axis in 1:3
51 stats = M.init_pmf_batchelor!(s;brms=0.1,p0=1,seed=14,kmax=4pi,
52 invariant_axis=axis)
53 K.synchronize(be)
54 initial = M.fields_to_host(s)
55 @test stats.invariant_axis == axis
56 for d in 6:8
57 field = reshape(initial[d],n,n,n)
58 @test maximum(abs,field.-circshift(field,ntuple(i->i==axis ? 1 : 0,3))) < 2e-7
59 @test maximum(abs,field) > 0.001
60 end
61 M.boost_gas_velocity!(s,(0.1,0.2,0.3))
62 M.step!(s,1e-4;ch=1,integrator=name === :cpu ? :ref : :cube)
63 K.synchronize(be)
64 for field in M.fields_to_host(s)
65 cube = reshape(field,n,n,n)
66 @test maximum(abs,cube.-circshift(cube,ntuple(i->i==axis ? 1 : 0,3))) < 2e-6
67 end
68 end
69 # Reading photon maps must preserve both the live Fourier modes and future
70 # dynamics. FFTW c2r may destroy input even for an out-of-place transform.
71 states = ntuple(2) do _
72 state = M.allocate_state(be,T,(n,n,n);dx=1/n,gamma=5/3,recon=:plm,riemann=:hll)
73 fill!(state.U[1],one(T));fill!(state.U[5],T(2))
74 M.boost_gas_velocity!(state,(0.125,0.0625,0.03125))
75 state
76 end
77 workspaces = map(M.allocate_radiation_workspace,states)
78 photons = map(M.allocate_photon_moment_state,states)
79 for d in 1:2
80 M.initialize_photon_moments!(photons[d],workspaces[d],states[d])
81 end
82 closure = M.RadiationClosure(inertia_ratio=.23,mean_free_path=70.,
83 photon_speed=3e4,gas_sound_speed=1.5,response_model=:moments_exponential)
84 for _ in 1:2
85 before = M.to_host(photons[1].hat_batch)
86 M.photon_fields_to_workspace!(workspaces[1],photons[1])
87 @test M.to_host(photons[1].hat_batch) == before
88 real = M.to_host(workspaces[1].real_batch)
89 @test maximum(abs,real[:,:,:,2].-T(.125)) < 2e-6
90 for d in 1:2
91 M.step_radiation_godunov!(states[d],workspaces[d],1e-4;
92 drag_impulse=.05,closure=closure,ch=2,
93 integrator=name === :cpu ? :ref : :cube,photon_state=photons[d])
94 end
95 end
96 for (left,right) in zip(M.fields_to_host(states[1]),M.fields_to_host(states[2]))
97 @test maximum(abs,left.-right) < 2e-6
98 end
99 @test maximum(abs,M.to_host(photons[1].hat_batch).-M.to_host(photons[2].hat_batch)) < 2e-6
100
101 end