Snapshot from this documentation build. SHA-256: 9d87274f54b75006dfe79518839a81ac136bbf34367fe2b289010c551b82b112
1 # Peculiar baryon minus DM velocities, sampled at the actual particle positions.
2 export boost_gas_velocity!, sample_baryon_dm_velocity!
3
4 @kernel function _boost_gas_velocity_k!(rho, mx, my, mz, E, ux, uy, uz)
5 p = @index(Global)
6 @inbounds begin
7 r = rho[p]
8 E[p] += ux*mx[p] + uy*my[p] + uz*mz[p] + r*(ux^2+uy^2+uz^2)/2
9 mx[p] += r*ux
10 my[p] += r*uy
11 mz[p] += r*uz
12 end
13 end
14
15 "Add a uniform code velocity while preserving gas internal and magnetic energy."
16 function boost_gas_velocity!(s::MHDState{T}, velocity) where T
17 length(velocity) == 3 && all(isfinite, velocity) ||
18 throw(ArgumentError("boost requires three finite components"))
19 _boost_gas_velocity_k!(s.be)(s.U[1],s.U[2],s.U[3],s.U[4],s.U[5],
20 T.(velocity)...; ndrange=ncells(s))
21 s
22 end
23
24 @kernel function _sample_baryon_dm_velocity_k!(outx, outy, outz,
25 @Const(rho), @Const(mx), @Const(my), @Const(mz),
26 @Const(px), @Const(py), @Const(pz), @Const(vx), @Const(vy), @Const(vz),
27 N::Int, lattice::Bool, gas_unit, dm_unit, smallr)
28 p = @index(Global)
29 @inbounds begin
30 T = eltype(outx)
31 q = p-1
32 # Work with offsets separately to retain tiny lattice displacements.
33 bx = lattice ? q % N : 0
34 by = lattice ? (q ÷ N) % N : 0
35 bz = lattice ? q ÷ (N*N) : 0
36 ux = T(N)*px[p] - (lattice ? zero(T) : T(0.5))
37 uy = T(N)*py[p] - (lattice ? zero(T) : T(0.5))
38 uz = T(N)*pz[p] - (lattice ? zero(T) : T(0.5))
39 ix = floor(Int,ux); iy = floor(Int,uy); iz = floor(Int,uz)
40 fx = ux-T(ix); fy = uy-T(iy); fz = uz-T(iz)
41 gx = zero(T); gy = zero(T); gz = zero(T)
42 for k in 0:1, j in 0:1, i in 0:1
43 c = mod(bx+ix+i,N) + N*mod(by+iy+j,N) + N*N*mod(bz+iz+k,N) + 1
44 w = (i==0 ? 1-fx : fx)*(j==0 ? 1-fy : fy)*(k==0 ? 1-fz : fz)
45 r = max(rho[c],smallr)
46 gx += w*mx[c]/r; gy += w*my[c]/r; gz += w*mz[c]/r
47 end
48 outx[p] = gas_unit*gx-dm_unit*vx[p]
49 outy[p] = gas_unit*gy-dm_unit*vy[p]
50 outz[p] = gas_unit*gz-dm_unit*vz[p]
51 end
52 end
53
54 """Sample v_b - v_c at DM positions using periodic CIC of gas velocity.
55
56 Outputs have the units supplied by gas_unit and dm_unit (normally cm/s).
57 Equal-mass particle averages are DM-mass weighted, not Eulerian volume averages;
58 these arrays must not be FFTed on the original lattice after particles move.
59 """
60 function sample_baryon_dm_velocity!(out, s, positions, velocities;
61 lattice_displacements::Bool=true, gas_unit::Real, dm_unit::Real)
62 N = s.dims[1]
63 s.dims == (N,N,N) || throw(ArgumentError("streaming sampler requires a cube"))
64 all(x->isfinite(x) && x>0,(gas_unit,dm_unit)) ||
65 throw(ArgumentError("velocity units must be finite and positive"))
66 length(out)==length(positions)==length(velocities)==3 ||
67 throw(ArgumentError("require three vector components"))
68 all(a->length(a)==ncells(s),(out...,positions...,velocities...)) ||
69 throw(ArgumentError("sampler requires one particle per original lattice cell"))
70 T = eltype(s.U[1])
71 _sample_baryon_dm_velocity_k!(s.be)(out...,s.U[1],s.U[2],s.U[3],s.U[4],
72 positions...,velocities...,N,lattice_displacements,T(gas_unit),T(dm_unit),s.smallr;
73 ndrange=ncells(s))
74 out
75 end