Snapshot from this documentation build. SHA-256: 2fdd553ab93a20b005066f010af41ce27a126afd0978aedfa512450720a53e71
1 # ── Face-staggered constrained transport (CT), an alternative to Dedner GLM ────
2 #
3 # SCHEME. The hydrodynamic part is shared with the GLM path: cell-centred Godunov,
4 # PLM(MonCen) with the full slope, a 3-direction MUSCL-Hancock predictor, and the
5 # selected Riemann solver (HLLD/HLL/LLF). What differs is the induction equation:
6 #
7 # * B lives on the FACES of the grid (`CTState.bxf/byf/bzf`; `bxf[c]` is Bx on the
8 # lower-x face (i-1/2,j,k) of cell c, etc.). The normal field at a face is that one
9 # continuous value: both Riemann states are overwritten with it, there is no
10 # Riemann problem for B_n, and no ψ (`riemann_ct`).
11 # * The Godunov transverse-B fluxes of the faces are the face EMFs. The edge EMF is
12 # built from the four surrounding face EMFs plus upwind (contact-mode) derivative
13 # corrections selected by the sign of the face mass flux, from the cell-centred
14 # reference EMF E = -v x B (v of the predicted half step, B at the start of the
15 # step) — the "CT-contact" scheme of Gardiner & Stone 2005 (`_ct_gs05`).
16 # * Faces are advanced by the discrete curl of the edge EMFs (`_ct_faceB_k!`), which
17 # keeps the face divergence at round-off for all time; cell-centred B (used by the
18 # hydro) is the average of the two faces (`_ct_bcc_k!`). Total energy is advanced by
19 # the Godunov energy flux alone (no energy correction).
20 #
21 # PROVENANCE. A KernelAbstractions, precision-generic port of the fixed CT prototype on
22 # the Solaceonx fork of GLMMHDTurb (`ct_mhd.jl`, branch `sweep-fixes`): commit 9d20e52
23 # (full MonCen slope + 3-direction Hancock predictor) and commit 7256715 (Gardiner &
24 # Stone 2005 upwinded edge EMFs with a cell-centred reference EMF). Gardiner & Stone
25 # 2005, J. Comput. Phys. 205, 509. Deviations from the fork: the reference EMF velocity
26 # includes the same drag/radiation source modifications the face states receive
27 # (`_facep_sources`); generic (Nx,Ny,Nz) grids; pmf-glm's pressure/density floors.
28 #
29 # MEMORY. 3 face arrays + 3 face-remainder arrays + 15 workspace arrays (6 face-flux arrays, 3 mass-flux arrays,
30 # 3 cell EMFs, 3 edge EMFs) of N^3 on top of GLM's 18 (9 state + 9 scratch). The ψ slot
31 # stays allocated but is identically zero.
32 #
33 # LIMITS. Periodic boundaries only; `:ref` integrator only (the shared-memory cube
34 # kernels error for CT states); the terminal-regime induction paths (which overwrite
35 # cell-centred B) error for CT states.
36
37 export max_divb_face, ct_faces_from_cells!, ct_faces_from_potential!
38
39
40 function _require_glm(s::MHDState, what::AbstractString)
41 uses_ct(s) && error("$what overwrites cell-centred B and is not available for " *
42 "mhd_solver=:ct states (constrained transport keeps B on faces)")
43 nothing
44 end
45
46 function _require_ct_integrator(s::MHDState, what::AbstractString)
47 uses_ct(s) && error("$what: CT is implemented on the :ref integrator only " *
48 "(mhd_solver=:ct states cannot use the :cube integrators)")
49 nothing
50 end
51
52 # ── Gardiner & Stone (2005) CT-contact edge EMF ──────────────────────────────
53 # The edge sits at (a-1/2, b-1/2) in the plane of the two directions a, b across it.
54 # Ea0, Eam are the EMFs on the a-faces (a-1/2, b) and (a-1/2, b-1); Eb0, Ebm those on
55 # the b-faces (a, b-1/2) and (a-1, b-1/2); mA*, mB* the mass fluxes through those
56 # faces; Er** the cell-centred EMFs of cells (a+da, b+db), da, db in {0,-1}.
57 @inline function _ct_gs05(Ea0::T,Eam::T,Eb0::T,Ebm::T,mA0::T,mAm::T,mB0::T,mBm::T,
58 Er00::T,Erm0::T,Er0m::T,Ermm::T) where {T}
59 h = T(0.5)
60 up(m,lower,upper) = ifelse(m>zero(T),lower,ifelse(m<zero(T),upper,h*(lower+upper)))
61 A1 = up(mA0,Erm0-Ebm,Er00-Eb0)
62 A3 = up(mAm,Ebm-Ermm,Eb0-Er0m)
63 B1 = up(mB0,Er0m-Eam,Er00-Ea0)
64 B3 = up(mBm,Eam-Ermm,Ea0-Erm0)
65 T(0.25)*((Ea0+Eam+Eb0+Ebm)+(A3-A1)+(B3-B1))
66 end
67
68 # Edge EMFs from the lower-face Godunov magnetic fluxes (E_z on x-faces = -F_x[By],
69 # on y-faces = +F_y[Bx]; E_x on y-faces = -F_y[Bz], on z-faces = +F_z[By];
70 # E_y on z-faces = -F_z[Bx], on x-faces = +F_x[Bz]).
71 @kernel function _ct_emf_k!(ex,ey,ez,@Const(mfx),@Const(mfy),@Const(mfz),
72 @Const(fxby),@Const(fxbz),@Const(fybx),@Const(fybz),@Const(fzbx),@Const(fzby),
73 @Const(ecx),@Const(ecy),@Const(ecz),Nx::Int,Ny::Int,Nz::Int)
74 c=@index(Global,Linear)
75 @inbounds begin
76 i=(c-1)%Nx+1; j=((c-1)÷Nx)%Ny+1; k=(c-1)÷(Nx*Ny)+1
77 L(a,b,d)=_periodic_corr_index(Nx,Ny,Nz,a,b,d)
78 c0=c; cim=L(i-1,j,k); cjm=L(i,j-1,k); ckm=L(i,j,k-1)
79 cimjm=L(i-1,j-1,k); cjmkm=L(i,j-1,k-1); ckmim=L(i-1,j,k-1)
80 # ez at (i-1/2,j-1/2,k): a=x, b=y
81 ez[c]=_ct_gs05(-fxby[c0],-fxby[cjm],fybx[c0],fybx[cim],
82 mfx[c0],mfx[cjm],mfy[c0],mfy[cim],
83 ecz[c0],ecz[cim],ecz[cjm],ecz[cimjm])
84 # ex at (i,j-1/2,k-1/2): a=y, b=z
85 ex[c]=_ct_gs05(-fybz[c0],-fybz[ckm],fzby[c0],fzby[cjm],
86 mfy[c0],mfy[ckm],mfz[c0],mfz[cjm],
87 ecx[c0],ecx[cjm],ecx[ckm],ecx[cjmkm])
88 # ey at (i-1/2,j,k-1/2): a=z, b=x
89 ey[c]=_ct_gs05(-fzbx[c0],-fzbx[cim],fxbz[c0],fxbz[ckm],
90 mfz[c0],mfz[cim],mfx[c0],mfx[ckm],
91 ecy[c0],ecy[ckm],ecy[cim],ecy[ckmim])
92 end
93 end
94
95 # Face update by the discrete curl of the edge EMFs.
96 @kernel function _ct_faceB_k!(bxf,byf,bzf,@Const(ex),@Const(ey),@Const(ez),
97 Nx::Int,Ny::Int,Nz::Int,dtdx::T) where {T}
98 c=@index(Global,Linear)
99 @inbounds begin
100 i=(c-1)%Nx+1; j=((c-1)÷Nx)%Ny+1; k=(c-1)÷(Nx*Ny)+1
101 L(a,b,d)=_periodic_corr_index(Nx,Ny,Nz,a,b,d)
102 ip=L(i+1,j,k); jp=L(i,j+1,k); kp=L(i,j,k+1)
103 bxf[c]-=dtdx*((ez[jp]-ez[c])-(ey[kp]-ey[c]))
104 byf[c]-=dtdx*((ex[kp]-ex[c])-(ez[ip]-ez[c]))
105 bzf[c]-=dtdx*((ey[ip]-ey[c])-(ex[jp]-ex[c]))
106 end
107 end
108
109 # The same curl update with Kahan (compensated) summation of every face, mirroring GLM's
110 # `_compensated_update!` for cell B: the part of each increment that the Float32 sum
111 # `b + inc` cannot represent is carried in a per-face remainder (`rbx/rby/rbz`) and added
112 # to the next increment instead of being lost. With a strong guide field the per-step
113 # increments `dtdx*curl(E)` of a weak perturbation are comparable to half an ulp of `b`,
114 # and plain summation drops them face by face; since the dropped parts differ from face to
115 # face the stored faces' divergence then random-walks away from round-off. Here the exact
116 # sum `stored + remainder` follows the curl update, so the divergence of the stored faces
117 # stays bounded by the remainders (at most half an ulp of `b` each).
118 @kernel function _ct_faceB_comp_k!(bxf,byf,bzf,rbx,rby,rbz,@Const(ex),@Const(ey),@Const(ez),
119 Nx::Int,Ny::Int,Nz::Int,dtdx::T) where {T}
120 c=@index(Global,Linear)
121 @inbounds begin
122 i=(c-1)%Nx+1; j=((c-1)÷Nx)%Ny+1; k=(c-1)÷(Nx*Ny)+1
123 L(a,b,d)=_periodic_corr_index(Nx,Ny,Nz,a,b,d)
124 ip=L(i+1,j,k); jp=L(i,j+1,k); kp=L(i,j,k+1)
125 bxf[c]=_compensated_update!(rbx,c,bxf[c],-(dtdx*((ez[jp]-ez[c])-(ey[kp]-ey[c]))),Val(true))
126 byf[c]=_compensated_update!(rby,c,byf[c],-(dtdx*((ex[kp]-ex[c])-(ez[ip]-ez[c]))),Val(true))
127 bzf[c]=_compensated_update!(rbz,c,bzf[c],-(dtdx*((ey[ip]-ey[c])-(ex[jp]-ex[c]))),Val(true))
128 end
129 end
130
131 # Cell-centred B = average of the two faces.
132 @kernel function _ct_bcc_k!(o6,o7,o8,@Const(bxf),@Const(byf),@Const(bzf),
133 Nx::Int,Ny::Int,Nz::Int)
134 c=@index(Global,Linear)
135 @inbounds begin
136 i=(c-1)%Nx+1; j=((c-1)÷Nx)%Ny+1; k=(c-1)÷(Nx*Ny)+1
137 L(a,b,d)=_periodic_corr_index(Nx,Ny,Nz,a,b,d)
138 T=eltype(o6)
139 o6[c]=T(0.5)*(bxf[c]+bxf[L(i+1,j,k)])
140 o7[c]=T(0.5)*(byf[c]+byf[L(i,j+1,k)])
141 o8[c]=T(0.5)*(bzf[c]+bzf[L(i,j,k+1)])
142 end
143 end
144
145 """
146 ct_update!(s, dtdx; compensate=false)
147
148 Second half of a CT step, run after the Godunov kernel and the `U`/`scratch` swap: edge
149 EMFs (Gardiner & Stone 2005 CT-contact) from the workspace fluxes, face update by the
150 discrete curl, and cell-centred B (`s.U[6:8]`) as the face average. `s.scratch` keeps the
151 full old state.
152
153 `compensate=true` sums every face update with Kahan compensation, carrying the remainders
154 in `s.ct.rbx/rby/rbz` (see `_ct_faceB_comp_k!`); the cell-centred B is still the plain
155 average of the stored faces. The radiation-coupled steps request it exactly when GLM
156 compensates cell B (`track_magnetic`).
157 """
158 function ct_update!(s::MHDState{T}, dtdx::Real; compensate::Bool=false) where {T}
159 ct=s.ct; ct === nothing && error("ct_update! needs a CT state")
160 Nx,Ny,Nz=s.dims; nc=ncells(s)
161 _ct_emf_k!(s.be,256)(ct.ex,ct.ey,ct.ez,ct.mfx,ct.mfy,ct.mfz,ct.fxby,ct.fxbz,
162 ct.fybx,ct.fybz,ct.fzbx,ct.fzby,ct.ecx,ct.ecy,ct.ecz,Nx,Ny,Nz;ndrange=nc)
163 if compensate
164 _ct_faceB_comp_k!(s.be,256)(ct.bxf,ct.byf,ct.bzf,ct.rbx,ct.rby,ct.rbz,
165 ct.ex,ct.ey,ct.ez,Nx,Ny,Nz,T(dtdx);ndrange=nc)
166 else
167 _ct_faceB_k!(s.be,256)(ct.bxf,ct.byf,ct.bzf,ct.ex,ct.ey,ct.ez,Nx,Ny,Nz,T(dtdx);
168 ndrange=nc)
169 end
170 _ct_bcc_k!(s.be,256)(s.U[6],s.U[7],s.U[8],ct.bxf,ct.byf,ct.bzf,Nx,Ny,Nz;ndrange=nc)
171 KA.synchronize(s.be)
172 s
173 end
174
175 # ── Diagnostics ──────────────────────────────────────────────────────────────
176 @kernel function _ct_divb_face_k!(d,@Const(bxf),@Const(byf),@Const(bzf),
177 Nx::Int,Ny::Int,Nz::Int,invdx::T) where {T}
178 c=@index(Global,Linear)
179 @inbounds begin
180 i=(c-1)%Nx+1; j=((c-1)÷Nx)%Ny+1; k=(c-1)÷(Nx*Ny)+1
181 L(a,b,e)=_periodic_corr_index(Nx,Ny,Nz,a,b,e)
182 d[c]=abs(((bxf[L(i+1,j,k)]-bxf[c])+(byf[L(i,j+1,k)]-byf[c])+
183 (bzf[L(i,j,k+1)]-bzf[c]))*invdx)
184 end
185 end
186
187 """
188 max_divb_face(s) -> Float64
189
190 Maximum over cells of the exact face divergence `|Σ_d (b_d[upper] - b_d[lower])/dx|` of a
191 CT state — the constraint the CT solver preserves to round-off. (The cell-centred
192 [`max_divb`](@ref) is a different, centred-difference quantity.)
193 """
194 function max_divb_face(s::MHDState{T}) where {T}
195 uses_ct(s) || error("max_divb_face needs a CT state (mhd_solver=:ct)")
196 Nx,Ny,Nz=s.dims; ct=s.ct
197 d=device_zeros(s.be,T,(ncells(s),))
198 _ct_divb_face_k!(s.be,256)(d,ct.bxf,ct.byf,ct.bzf,Nx,Ny,Nz,one(T)/s.dx;ndrange=ncells(s))
199 KA.synchronize(s.be)
200 Float64(maximum(d))
201 end
202
203 # ── Initialization ───────────────────────────────────────────────────────────
204 # cell B <- face average; E <- E + (|B_new|^2 - |B_old|^2)/2 (thermal pressure fixed); psi <- 0
205 @kernel function _ct_cell_from_faces_k!(U5,U6,U7,U8,U9,@Const(bxf),@Const(byf),@Const(bzf),
206 Nx::Int,Ny::Int,Nz::Int)
207 c=@index(Global,Linear)
208 @inbounds begin
209 T=eltype(U5)
210 i=(c-1)%Nx+1; j=((c-1)÷Nx)%Ny+1; k=(c-1)÷(Nx*Ny)+1
211 L(a,b,e)=_periodic_corr_index(Nx,Ny,Nz,a,b,e)
212 bx=T(0.5)*(bxf[c]+bxf[L(i+1,j,k)])
213 by=T(0.5)*(byf[c]+byf[L(i,j+1,k)])
214 bz=T(0.5)*(bzf[c]+bzf[L(i,j,k+1)])
215 ox=U6[c]; oy=U7[c]; oz=U8[c]
216 U5[c]+=T(0.5)*((bx*bx+by*by+bz*bz)-(ox*ox+oy*oy+oz*oz))
217 U6[c]=bx; U7[c]=by; U8[c]=bz; U9[c]=zero(T)
218 end
219 end
220
221 "Zero the CT face remainders (no-op for GLM states)."
222 function _ct_reset_face_residual!(s::MHDState)
223 ct=s.ct
224 ct === nothing || foreach(r -> fill!(r,zero(eltype(r))),(ct.rbx,ct.rby,ct.rbz))
225 nothing
226 end
227
228 function _ct_install_faces!(s::MHDState{T}, bxf, byf, bzf) where {T}
229 ct=s.ct; Nx,Ny,Nz=s.dims
230 _ct_reset_face_residual!(s)
231 copyto!(ct.bxf,bxf); copyto!(ct.byf,byf); copyto!(ct.bzf,bzf)
232 _ct_cell_from_faces_k!(s.be,256)(s.U[5],s.U[6],s.U[7],s.U[8],s.U[9],
233 ct.bxf,ct.byf,ct.bzf,Nx,Ny,Nz;ndrange=ncells(s))
234 KA.synchronize(s.be)
235 s
236 end
237
238 # Half-cell phase θ = k dx/2 of rfft/fft index m (0-based) on an N-point axis; Nyquist
239 # (even N, m = N/2) is reported separately because cos θ = 0 there.
240 @inline function _ct_half_phase(m::Int,N::Int)
241 nyq = iseven(N) && m == N÷2
242 mm = m > N/2 ? m-N : m
243 (nyq ? 0.0 : π*mm/N, nyq)
244 end
245
246 """
247 ct_faces_from_cells!(s::MHDState) -> NamedTuple
248
249 Build a discretely divergence-free face field from the cell-centred `B = s.U[6:8]` of a CT
250 state (host Float64 FFTW), changing the cell field as little as possible.
251
252 The CT cell field is the average of its two faces, so with `θ_d = k_d dx/2` a face mode
253 `f_d = b_d exp(-iθ_d)/cos θ_d` averages back to exactly `b_d`, and its face divergence
254 `Σ_d (exp(2iθ_d)-1) f_d/dx = (2i/dx) Σ_d tan θ_d b_d`. The cell field is therefore
255 projected orthogonally onto `t·b = 0` with `t_d = tan θ_d` (the smallest L2 change of the
256 cell field that admits divergence-free faces), and the faces are built from the result.
257 On a Nyquist plane of axis d (`cos θ_d = 0`) the face average of component d vanishes for
258 any face field, so that component is dropped there. Both operations are
259 Hermitian-symmetric, so the faces are real and divergence-free to round-off. (A plain
260 half-cell shift of the cell field followed by a face projection instead filters every
261 component by `cos θ_d` on the way back to the cell average: 7% of the magnetic energy of a
262 Batchelor field cut at 8 cells per wavelength.)
263
264 The centred-curl and continuum PMF fields are solenoidal under different symbols
265 (`sin 2θ`, `θ`) than `tan θ`, so the projection changes them by O((k dx)^2) near the
266 cutoff. Then `s.U[6:8]` become the face averages, `s.U[5]` is adjusted so the thermal
267 pressure is unchanged (`E += ½|B_new|² - ½|B_old|²`), and `s.U[9] = 0`.
268
269 Returns `(cell_b_rel_l2_change, magnetic_energy_ratio, max_divb_face_dx_over_brms)`.
270 """
271 function ct_faces_from_cells!(s::MHDState{T}) where {T}
272 uses_ct(s) || error("ct_faces_from_cells! needs a CT state (mhd_solver=:ct)")
273 Nx,Ny,Nz=s.dims; dx=Float64(s.dx)
274 Bold=ntuple(d->reshape(Float64.(to_host(s.U[5+d])),Nx,Ny,Nz),3)
275 Bh=ntuple(d->rfft(Bold[d]),3)
276 Fh=ntuple(d->similar(Bh[d]),3)
277 nxr=Nx÷2+1
278 px=[_ct_half_phase(m,Nx) for m in 0:nxr-1]
279 py=[_ct_half_phase(m,Ny) for m in 0:Ny-1]
280 pz=[_ct_half_phase(m,Nz) for m in 0:Nz-1]
281 for kk in 1:Nz, jj in 1:Ny, ii in 1:nxr
282 (θx,nx),(θy,ny),(θz,nz)=px[ii],py[jj],pz[kk]
283 bx=nx ? zero(ComplexF64) : Bh[1][ii,jj,kk]
284 by=ny ? zero(ComplexF64) : Bh[2][ii,jj,kk]
285 bz=nz ? zero(ComplexF64) : Bh[3][ii,jj,kk]
286 tx=tan(θx); ty=tan(θy); tz=tan(θz)
287 t2=tx*tx+ty*ty+tz*tz
288 if t2>0
289 dot=(tx*bx+ty*by+tz*bz)/t2
290 bx-=tx*dot; by-=ty*dot; bz-=tz*dot
291 end
292 Bh[1][ii,jj,kk]=bx; Bh[2][ii,jj,kk]=by; Bh[3][ii,jj,kk]=bz
293 Fh[1][ii,jj,kk]=bx*cis(-θx)/cos(θx)
294 Fh[2][ii,jj,kk]=by*cis(-θy)/cos(θy)
295 Fh[3][ii,jj,kk]=bz*cis(-θz)/cos(θz)
296 end
297 faces=ntuple(d->vec(irfft(Fh[d],Nx)),3)
298 _ct_install_faces!(s,to_device(s.be,faces[1],T),to_device(s.be,faces[2],T),
299 to_device(s.be,faces[3],T))
300 Bnew=ntuple(d->reshape(Float64.(to_host(s.U[5+d])),Nx,Ny,Nz),3)
301 num=sum(sum(abs2,Bnew[d].-Bold[d]) for d in 1:3)
302 den=sum(sum(abs2,Bold[d]) for d in 1:3)
303 enew=sum(sum(abs2,Bnew[d]) for d in 1:3)
304 brms=sqrt(enew/length(Bnew[1]))
305 (cell_b_rel_l2_change = den>0 ? sqrt(num/den) : 0.0,
306 magnetic_energy_ratio = den>0 ? enew/den : 1.0,
307 max_divb_face_dx_over_brms = brms>0 ? max_divb_face(s)*dx/brms : 0.0)
308 end
309
310 @kernel function _ct_curl_k!(bxf,byf,bzf,@Const(Ax),@Const(Ay),@Const(Az),
311 Nx::Int,Ny::Int,Nz::Int,invdx::T) where {T}
312 c=@index(Global,Linear)
313 @inbounds begin
314 i=(c-1)%Nx+1; j=((c-1)÷Nx)%Ny+1; k=(c-1)÷(Nx*Ny)+1
315 L(a,b,e)=_periodic_corr_index(Nx,Ny,Nz,a,b,e)
316 bxf[c]=((Az[L(i,j+1,k)]-Az[c])-(Ay[L(i,j,k+1)]-Ay[c]))*invdx
317 byf[c]=((Ax[L(i,j,k+1)]-Ax[c])-(Az[L(i+1,j,k)]-Az[c]))*invdx
318 bzf[c]=((Ay[L(i+1,j,k)]-Ay[c])-(Ax[L(i,j+1,k)]-Ax[c]))*invdx
319 end
320 end
321
322 """
323 ct_faces_from_potential!(s::MHDState, Ax, Ay, Az)
324
325 Set the CT face field to the discrete curl of an EDGE-centred vector potential (arrays of
326 size `s.dims` or length `ncells`; cell `c=(i,j,k)` holds `Ax` at `(i,j-1/2,k-1/2)`, `Ay` at
327 `(i-1/2,j,k-1/2)`, `Az` at `(i-1/2,j-1/2,k)`, with the same lower-face convention as the
328 face fields). The face divergence is zero to round-off by construction. Then, as in
329 [`ct_faces_from_cells!`](@ref), cell B becomes the face average, `E` is adjusted so the
330 thermal pressure is unchanged, and `ψ = 0`.
331 """
332 function ct_faces_from_potential!(s::MHDState{T}, Ax, Ay, Az) where {T}
333 uses_ct(s) || error("ct_faces_from_potential! needs a CT state (mhd_solver=:ct)")
334 Nx,Ny,Nz=s.dims
335 for A in (Ax,Ay,Az)
336 length(A)==ncells(s) || error("vector-potential array has $(length(A)) entries, expected $(ncells(s))")
337 end
338 dev(A)=vec(to_device(s.be,A isa Array ? A : Array(A),T))
339 dAx=dev(Ax); dAy=dev(Ay); dAz=dev(Az)
340 nc=ncells(s)
341 bx=device_zeros(s.be,T,(nc,)); by=device_zeros(s.be,T,(nc,)); bz=device_zeros(s.be,T,(nc,))
342 _ct_curl_k!(s.be,256)(bx,by,bz,dAx,dAy,dAz,Nx,Ny,Nz,one(T)/s.dx;ndrange=nc)
343 KA.synchronize(s.be)
344 _ct_install_faces!(s,bx,by,bz)
345 end