Snapshot from this documentation build. SHA-256: d167f8c4038ef3fff3ae527ba3a7cecf7a714475f84299af251a144a466397e2
1 """
2 MHDKernels
3
4 A KernelAbstractions.jl ideal-MHD solver with **Dedner mixed-GLM divergence
5 cleaning** — MUSCL-Hancock + PLM(MonCen) reconstruction + HLLD (LLF fallback),
6 written **once** and run on the CPU (the f64 convergence/conservation oracle) and
7 on CUDA / Metal GPUs (f32). The overdamped PMF path is separate: it advances B
8 with a curl-form constrained-transport update and does not require GLM cleaning
9 or a projection FFT while terminal mode is active.
10
11 Design contract (mirrors `PPMKernels`):
12
13 * **One source, two+ devices, f32-first.** Every compute kernel is a
14 precision-generic `@kernel` parameterised on the element type `T`, but the
15 DEFAULT precision is `Float32` on EVERY backend — the CPU runs f32 too, so
16 CPU↔GPU is an apples-to-apples comparison. `Float64` is available (pass it to
17 `allocate_state`) only when a physical reason needs the dynamic range. The
18 convergence gate runs in f32 (finite-amplitude exact Alfvén wave keeps the
19 discretisation error above the f32 round-off floor).
20 * **Backend by name.** `backend(:cpu)` always works; `backend(:cuda)` /
21 `backend(:metal)` resolve after `using CUDA` / `using Metal` loads the matching
22 package extension. Allocation/transfer go through [`device_zeros`](@ref) /
23 [`to_device`](@ref) / [`to_host`](@ref).
24 * **Conserved order** `U = (ρ, ρvx, ρvy, ρvz, E, Bx, By, Bz, ψ)`;
25 **primitive order** `W = (ρ, vx, vy, vz, p, Bx, By, Bz, ψ)`.
26
27 Optional **constrained transport** (`allocate_state(...; mhd_solver=:ct)`): a
28 face-staggered, Gardiner-Stone-2005 upwinded-EMF alternative to GLM for the `:ref`
29 integrator on periodic grids (see `ct.jl`); GLM remains the default.
30
31 Two integrators, cross-checked: a portable per-cell **reference** kernel
32 (`step_ref!`, runs on every backend incl. CPU-f64) and a shared-memory **cube**
33 kernel (`step_cube!`, GPU throughput path). Both implement the identical scheme.
34 """
35 module MHDKernels
36
37 using FFTW
38 using KernelAbstractions
39 using Random
40 using StaticArrays
41 const KA = KernelAbstractions
42
43 export backend, has_backend, device_zeros, to_device, to_host
44
45 # ── backend registry ─────────────────────────────────────────────────────────
46 const _BACKENDS = Dict{Symbol,Any}(:cpu => CPU())
47 register_backend!(name::Symbol, be) = (_BACKENDS[name] = be)
48 has_backend(name::Symbol) = haskey(_BACKENDS, name)
49
50 """
51 backend(name::Symbol = :cpu)
52
53 The KernelAbstractions backend registered under `name`. `:cpu` is always
54 available; `:cuda`/`:metal` require `using CUDA`/`using Metal` first.
55 """
56 function backend(name::Symbol = :cpu)
57 return get(_BACKENDS, name) do
58 error("MHD backend :$name is not available. " *
59 (name === :cuda ? "Run `using CUDA` first." :
60 name === :metal ? "Run `using Metal` first (Apple Silicon only)." :
61 "Known backends: $(collect(keys(_BACKENDS)))."))
62 end
63 end
64
65 # ── device array helpers (specialised by the CUDA/Metal extensions) ──────────
66 device_zeros(::CPU, ::Type{T}, dims::Dims) where {T} = zeros(T, dims)
67
68 function to_device(be, a::AbstractArray, ::Type{T} = eltype(a)) where {T}
69 d = device_zeros(be, T, size(a))
70 copyto!(d, convert(Array{T}, a))
71 return d
72 end
73
74 "`to_host(a)` — a plain host `Array` copy of a device array; synchronizes first."
75 function to_host(a::AbstractArray)
76 KA.synchronize(KA.get_backend(a))
77 return Array(a)
78 end
79
80 include("glm_physics.jl") # precision-generic @fastmath GLM-MHD device functions
81 include("boundary.jl") # per-axis BCs (periodic/outflow/reflecting), ghost-free
82 include("state.jl") # MHDState (9 SoA fields) + cons/prim host helpers
83 include("passive_transport.jl") # allocation-free bounded chemistry-species remap
84 include("integrator_ref.jl") # portable per-cell reference integrator (CPU + GPU)
85 include("ct.jl") # face-staggered constrained transport (alternative to GLM)
86 include("integrator_cube.jl") # fused shared-memory cube integrator (GPU throughput)
87 include("timestep.jl") # CFL Δt reduction
88 include("diagnostics.jl") # max|∇·B|, conserved totals
89 include("problems.jl") # initial conditions (smooth waves, Brio-Wu, Orszag-Tang)
90 include("streaming.jl") # coherent boosts and particle-position relative velocities
91 include("pmf.jl") # divergence-free primordial magnetic field initializers
92 include("terminal_ct.jl") # divergence-preserving terminal-regime induction
93 include("terminal_pressure.jl") # pressure-implicit terminal response + CT coupling
94 include("gravity_coupling.jl") # time-centered gas source for particle gravity
95 include("driver.jl") # step! / evolve! (ch + ψ-damping schedule)
96 include("drag.jl") # exact linear drag and Strang-coupled full-MHD step
97 include("radiation_coupling.jl") # spectral photon-baryon response inside full MHD
98 include("recombination_transport.jl") # nonlocal Lyman-alpha chemistry response
99 include("radiation_terminal_ct.jl") # B-only overdamped radiation response + CT
100
101 end # module