Snapshot from this documentation build. SHA-256: 54daec92822e01c3ac2aa7399f475743a774951dc0375e4a3f0641d2719a006c
1 # ── MHDState: 9 SoA conserved fields + a ping-pong scratch set + run parameters ─
2 # Fields are flat column-major length-Nx·Ny·Nz device vectors; cell (i,j,k) lives at
3 # idx = (k-1)*Nx*Ny + (j-1)*Nx + i (1-based).
4 # Boundaries are PERIODIC for now (the smooth-wave convergence gate needs nothing
5 # else); `bc` is carried for the outflow/reflecting extension. The conserved order is
6 # U=(ρ,ρvx,ρvy,ρvz,E,Bx,By,Bz,ψ); ψ is the GLM cleaning potential (damped, not
7 # conserved). γ is the adiabatic index; the floors/Riemann-switch params feed
8 # `riemann`. `be` is the KA backend the fields live on.
9
10 export MHDState, allocate_state, nvars, ncells, linidx, fields_to_host
11 export CTState, uses_ct, mhd_solver
12
13 const NVAR = 9
14
15 """
16 CTState{A}
17
18 Face-staggered constrained-transport (CT) storage attached to an [`MHDState`](@ref)
19 allocated with `mhd_solver=:ct` (see `ct.jl` for the scheme). All fields are flat
20 length-`ncells` device vectors indexed like the cell fields, cell `c=(i,j,k)`:
21
22 * `bxf, byf, bzf` — the evolved face fields. `bxf[c]` is `Bx` on the LOWER-x face
23 `(i-1/2,j,k)` of cell `c`; likewise `byf` at `(i,j-1/2,k)`, `bzf` at `(i,j,k-1/2)`.
24 * `mfx, mfy, mfz` — per-step mass flux through the lower x/y/z face of each cell.
25 * `fxby, fxbz` — lower-x-face magnetic fluxes `F_x[By]`, `F_x[Bz]`;
26 `fybx, fybz` — lower-y-face `F_y[Bx]`, `F_y[Bz]`; `fzbx, fzby` — lower-z-face
27 `F_z[Bx]`, `F_z[By]`.
28 * `ecx, ecy, ecz` — cell-centred reference EMF `E = -v x B`.
29 * `ex, ey, ez` — edge EMFs: `ez[c]` at `(i-1/2,j-1/2,k)`, `ex[c]` at
30 `(i,j-1/2,k-1/2)`, `ey[c]` at `(i-1/2,j,k-1/2)`.
31
32 * `rbx, rby, rbz` — Kahan remainders of the face updates (slots 19-21). They are used
33 only when `ct_update!` is asked to compensate (the radiation-coupled steps with
34 `MHD_RADIATION_MAGNETIC_COMPENSATION`, like GLM's cell-B remainders) and stay zero
35 otherwise.
36
37 All 21 are contiguous views, in this order, of ONE device buffer `packed` of length
38 `21*ncells`: the fused Godunov kernel takes that single buffer (slot `q` of cell `c` at
39 `(q-1)*ncells + c`), because Metal allows a kernel 31 buffer arguments and counts every
40 scalar argument as one: with 18 separate CT arrays it failed with 38/31.
41
42 (The struct lives in `state.jl` because `MHDState` has a field of this type.)
43 """
44 struct CTState{A,V}
45 packed::A
46 bxf::V; byf::V; bzf::V
47 mfx::V; mfy::V; mfz::V
48 fxby::V; fxbz::V; fybx::V; fybz::V; fzbx::V; fzby::V
49 ecx::V; ecy::V; ecz::V
50 ex::V; ey::V; ez::V
51 rbx::V; rby::V; rbz::V
52 end
53
54 function CTState(packed::AbstractVector, nc::Int)
55 length(packed) == 21nc || error("CT packed buffer has $(length(packed)) entries, expected $(21nc)")
56 v = ntuple(q -> view(packed, (q-1)*nc+1:q*nc), 21)
57 CTState{typeof(packed),eltype(v)}(packed, v...)
58 end
59
60 "`ct_arrays(ct)` — the 21 CT field views (18 fields + 3 face remainders) as a tuple, in field order."
61 ct_arrays(c::CTState) = ntuple(i -> getfield(c, i + 1), 21)
62
63 mutable struct MHDState{T,A,B}
64 U::NTuple{NVAR,A} # conserved
65 scratch::NTuple{NVAR,A} # ping-pong target for the integrator
66 dims::NTuple{3,Int} # (Nx, Ny, Nz)
67 dx::T
68 γ::T
69 smallr::T
70 pfl::T
71 llf_dmin::T
72 llf_pmin::T
73 use_hlld::Bool
74 riemann::Symbol
75 bcs::NTuple{3,Symbol} # per-axis: :periodic / :outflow / :reflecting
76 recon::Symbol # :plm (MonCen), :plm_minmod, or :ppm (local CW84)
77 be::B
78 ct::Union{Nothing,CTState{A}} # (CTState{A,V} for any view type V) face-staggered CT storage (`mhd_solver=:ct`), else nothing
79 end
80
81 # GLM-only construction (the historic 14-field form).
82 MHDState{T,A,B}(U, scratch, dims, dx, γ, smallr, pfl, llf_dmin, llf_pmin, use_hlld,
83 riemann, bcs, recon, be) where {T,A,B} =
84 MHDState{T,A,B}(U, scratch, dims, dx, γ, smallr, pfl, llf_dmin, llf_pmin, use_hlld,
85 riemann, bcs, recon, be, nothing)
86
87 "`uses_ct(s)` — true when `s` carries face-staggered CT storage (`mhd_solver=:ct`)."
88 uses_ct(s::MHDState) = s.ct !== nothing
89
90 "`mhd_solver(s)` — `:ct` (constrained transport) or `:glm` (Dedner GLM cleaning, the default)."
91 mhd_solver(s::MHDState) = uses_ct(s) ? :ct : :glm
92
93 nvars(::MHDState) = NVAR
94 ncells(s::MHDState) = prod(s.dims)
95 @inline linidx(Nx::Int, Ny::Int, i::Int, j::Int, k::Int) = ((k-1)*Ny + (j-1))*Nx + i
96
97 "`bc_codes(s)` → the per-axis BC integer codes (cx,cy,cz) for the kernels."
98 bc_codes(s::MHDState) = (bc_code(s.bcs[1]), bc_code(s.bcs[2]), bc_code(s.bcs[3]))
99
100 "`recon_code_of(s)` → reconstruction code (0=MonCen PLM, 1=PPM, 2=minmod PLM)."
101 recon_code_of(s::MHDState) = recon_code(s.recon)
102
103 "`riemann_code_of(s)` → the compile-time Riemann code (:hlld, :hll, :llf)."
104 riemann_code_of(s::MHDState) = riemann_code(s.riemann)
105
106 "`all_periodic(s)` → true when every axis is periodic (enables the cube's fast load path)."
107 all_periodic(s::MHDState) = all(==(:periodic), s.bcs)
108
109 """
110 allocate_state(be, T, dims; dx, gamma=5/3, smallr=1e-6, pfl=1e-7,
111 llf_dmin=1e-4, llf_pmin=0, use_hlld=true, riemann=nothing,
112 bc=:periodic, bcs=(bc,bc,bc), recon=:plm, mhd_solver=:glm)
113
114 Allocate an `MHDState` of element type `T` and shape `dims=(Nx,Ny,Nz)` on backend
115 `be` (zero-filled; fill it with a problem initialiser). `dx` is the uniform cell
116 width (domain assumed `[0,L]^3` with L = Nx·dx). `bc` sets all three axes; `bcs`
117 overrides per-axis (e.g. `bcs=(:outflow,:periodic,:periodic)` for a 1-D shock tube).
118
119 `mhd_solver` selects the divergence treatment: `:glm` (default; Dedner cleaning with
120 the ψ field) or `:ct` (face-staggered constrained transport, see `ct.jl`; allocates
121 3 face + 15 workspace + 3 face-remainder arrays and requires all-periodic BCs and the `:ref` integrator).
122 """
123 function allocate_state(be, ::Type{T}, dims::NTuple{3,Int};
124 dx::Real, gamma::Real = 5//3,
125 smallr::Real = 1e-6, pfl::Real = 1e-7,
126 llf_dmin::Real = 1e-4, llf_pmin::Real = 0,
127 use_hlld::Bool = true, riemann = nothing, bc::Symbol = :periodic,
128 bcs::NTuple{3,Symbol} = (bc, bc, bc), recon::Symbol = :plm,
129 mhd_solver::Symbol = :glm) where {T}
130 mhd_solver in (:glm, :ct) ||
131 error("unknown mhd_solver :$mhd_solver (have :glm, :ct)")
132 mhd_solver === :ct && !all(==(:periodic), bcs) &&
133 error("mhd_solver=:ct requires all-periodic boundaries (bcs=$bcs); constrained transport is periodic-only")
134 nc = prod(dims)
135 U = ntuple(_ -> device_zeros(be, T, (nc,)), NVAR)
136 sc = ntuple(_ -> device_zeros(be, T, (nc,)), NVAR)
137 rsym = riemann === nothing ? (use_hlld ? :hlld : :llf) : Symbol(riemann)
138 ct = mhd_solver === :ct ? CTState(device_zeros(be, T, (21nc,)), nc) : nothing
139 MHDState{T,eltype(U),typeof(be)}(U, sc, dims, T(dx), T(gamma), T(smallr), T(pfl),
140 T(llf_dmin), T(llf_pmin), use_hlld, rsym, bcs, recon, be, ct)
141 end
142
143 "`fields_to_host(s)` → an `NTuple{9,Array}` host copy of the conserved fields."
144 fields_to_host(s::MHDState) = ntuple(v -> to_host(s.U[v]), NVAR)