Snapshot from this documentation build. SHA-256: 8da043b1ae0326ae29c96d74ba042503c418d60f5df42cb8dd910d76d2633d59
1 # Divergence-preserving terminal-regime induction for early-PMF evolution.
2 # The pressure response is supplied by the caller; this module advances B with
3 # SSPRK2 using a centered curl plus LLF-equivalent resistivity. On a periodic
4 # uniform grid, div(curl(.))=0 and div(laplacian(B))=laplacian(div(B)) for the
5 # matching centered operators, so an initially solenoidal field stays so without
6 # a projection FFT.
7
8 export terminal_ct_induction!, terminal_magnetic_force!,
9 terminal_magnetic_force_split!, terminal_velocity_max!
10
11 @inline _tct_lin3(i, j, k, N) = begin
12 ii = mod(i - 1, N) + 1
13 jj = mod(j - 1, N) + 1
14 kk = mod(k - 1, N) + 1
15 ((kk - 1) * N + (jj - 1)) * N + ii
16 end
17
18 @inline function _tct_magnetic_force_at(Bx, By, Bz, i, j, k, N, inv2dx)
19 T = typeof(inv2dx)
20 cxp = _tct_lin3(i + 1, j, k, N); cxm = _tct_lin3(i - 1, j, k, N)
21 cyp = _tct_lin3(i, j + 1, k, N); cym = _tct_lin3(i, j - 1, k, N)
22 czp = _tct_lin3(i, j, k + 1, N); czm = _tct_lin3(i, j, k - 1, N)
23
24 bxp = Bx[cxp]; byp = By[cxp]; bzp = Bz[cxp]
25 bxm = Bx[cxm]; bym = By[cxm]; bzm = Bz[cxm]
26 b2p = bxp*bxp + byp*byp + bzp*bzp
27 b2m = bxm*bxm + bym*bym + bzm*bzm
28 dtxx = (bxp*bxp - T(0.5)*b2p) - (bxm*bxm - T(0.5)*b2m)
29 dtyx = byp*bxp - bym*bxm
30 dtzx = bzp*bxp - bzm*bxm
31
32 bxp = Bx[cyp]; byp = By[cyp]; bzp = Bz[cyp]
33 bxm = Bx[cym]; bym = By[cym]; bzm = Bz[cym]
34 b2p = bxp*bxp + byp*byp + bzp*bzp
35 b2m = bxm*bxm + bym*bym + bzm*bzm
36 dtxy = bxp*byp - bxm*bym
37 dtyy = (byp*byp - T(0.5)*b2p) - (bym*bym - T(0.5)*b2m)
38 dtzy = bzp*byp - bzm*bym
39
40 bxp = Bx[czp]; byp = By[czp]; bzp = Bz[czp]
41 bxm = Bx[czm]; bym = By[czm]; bzm = Bz[czm]
42 b2p = bxp*bxp + byp*byp + bzp*bzp
43 b2m = bxm*bxm + bym*bym + bzm*bzm
44 dtxz = bxp*bzp - bxm*bzm
45 dtyz = byp*bzp - bym*bzm
46 dtzz = (bzp*bzp - T(0.5)*b2p) - (bzm*bzm - T(0.5)*b2m)
47
48 return (dtxx + dtxy + dtxz) * inv2dx,
49 (dtyx + dtyy + dtyz) * inv2dx,
50 (dtzx + dtzy + dtzz) * inv2dx
51 end
52
53 @inline function _tct_magnetic_force_split_at(Bx, By, Bz, i, j, k, N, inv2dx)
54 T = typeof(inv2dx)
55 fx, fy, fz = _tct_magnetic_force_at(Bx, By, Bz, i, j, k, N, inv2dx)
56 cxp = _tct_lin3(i + 1, j, k, N); cxm = _tct_lin3(i - 1, j, k, N)
57 cyp = _tct_lin3(i, j + 1, k, N); cym = _tct_lin3(i, j - 1, k, N)
58 czp = _tct_lin3(i, j, k + 1, N); czm = _tct_lin3(i, j, k - 1, N)
59 b2xp = Bx[cxp]*Bx[cxp] + By[cxp]*By[cxp] + Bz[cxp]*Bz[cxp]
60 b2xm = Bx[cxm]*Bx[cxm] + By[cxm]*By[cxm] + Bz[cxm]*Bz[cxm]
61 b2yp = Bx[cyp]*Bx[cyp] + By[cyp]*By[cyp] + Bz[cyp]*Bz[cyp]
62 b2ym = Bx[cym]*Bx[cym] + By[cym]*By[cym] + Bz[cym]*Bz[cym]
63 b2zp = Bx[czp]*Bx[czp] + By[czp]*By[czp] + Bz[czp]*Bz[czp]
64 b2zm = Bx[czm]*Bx[czm] + By[czm]*By[czm] + Bz[czm]*Bz[czm]
65 px = -T(0.5) * (b2xp - b2xm) * inv2dx
66 py = -T(0.5) * (b2yp - b2ym) * inv2dx
67 pz = -T(0.5) * (b2zp - b2zm) * inv2dx
68 return px, py, pz, fx - px, fy - py, fz - pz
69 end
70
71 @kernel function _tct_magnetic_force_k!(fx, fy, fz,
72 @Const(Bx), @Const(By), @Const(Bz),
73 N::Int, inv2dx)
74 c = @index(Global)
75 @inbounds begin
76 i = (c - 1) % N + 1
77 j = ((c - 1) ÷ N) % N + 1
78 k = (c - 1) ÷ (N * N) + 1
79 x, y, z = _tct_magnetic_force_at(Bx, By, Bz, i, j, k, N, inv2dx)
80 fx[c] = x; fy[c] = y; fz[c] = z
81 end
82 end
83
84 @kernel function _tct_magnetic_force_split_k!(px, py, pz, tx, ty, tz,
85 @Const(Bx), @Const(By), @Const(Bz),
86 N::Int, inv2dx)
87 c = @index(Global)
88 @inbounds begin
89 i = (c - 1) % N + 1
90 j = ((c - 1) ÷ N) % N + 1
91 k = (c - 1) ÷ (N * N) + 1
92 pxi, pyi, pzi, txi, tyi, tzi =
93 _tct_magnetic_force_split_at(Bx, By, Bz, i, j, k, N, inv2dx)
94 px[c] = pxi; py[c] = pyi; pz[c] = pzi
95 tx[c] = txi; ty[c] = tyi; tz[c] = tzi
96 end
97 end
98
99 function terminal_magnetic_force!(s::MHDState, fx, fy, fz;
100 B=(s.U[6], s.U[7], s.U[8]))
101 N = s.dims[1]
102 all(==(N), s.dims) || error("terminal magnetic force requires a cubic grid")
103 _tct_magnetic_force_k!(s.be)(fx, fy, fz, B..., N, inv(typeof(s.dx)(2) * s.dx);
104 ndrange=ncells(s))
105 return fx, fy, fz
106 end
107
108 """
109 terminal_magnetic_force_split!(s, px, py, pz, tx, ty, tz; B=...)
110
111 Compute the centered conservative Maxwell-stress force as magnetic pressure
112 `-grad(B^2/2)` plus the remaining tension contribution. The remainder uses the
113 same discrete total-force stencil as [`terminal_magnetic_force!`](@ref), so the
114 two terms reconstruct the force used by the solver. For a solenoidal field the
115 tension contribution approaches `(B dot grad)B`.
116 """
117 function terminal_magnetic_force_split!(s::MHDState, px, py, pz, tx, ty, tz;
118 B=(s.U[6], s.U[7], s.U[8]))
119 N = s.dims[1]
120 all(==(N), s.dims) || error("terminal magnetic force requires a cubic grid")
121 _tct_magnetic_force_split_k!(s.be)(px, py, pz, tx, ty, tz, B..., N,
122 inv(typeof(s.dx)(2) * s.dx);
123 ndrange=ncells(s))
124 return px, py, pz, tx, ty, tz
125 end
126
127 @inline function _tct_velocity_at(rho, fx, fy, fz, i, j, k, N, inv2dx,
128 gamma_drag, pressure_coeff, smallr, vcap)
129 c = _tct_lin3(i, j, k, N)
130 cxp = _tct_lin3(i + 1, j, k, N); cxm = _tct_lin3(i - 1, j, k, N)
131 cyp = _tct_lin3(i, j + 1, k, N); cym = _tct_lin3(i, j - 1, k, N)
132 czp = _tct_lin3(i, j, k + 1, N); czm = _tct_lin3(i, j, k - 1, N)
133 r = max(rho[c], smallr)
134 gx = pressure_coeff * (rho[cxp] - rho[cxm]) * inv2dx
135 gy = pressure_coeff * (rho[cyp] - rho[cym]) * inv2dx
136 gz = pressure_coeff * (rho[czp] - rho[czm]) * inv2dx
137 invrg = one(gamma_drag) / (r * max(gamma_drag, eps(gamma_drag)))
138 vx = (fx[c] - gx) * invrg
139 vy = (fy[c] - gy) * invrg
140 vz = (fz[c] - gz) * invrg
141 v2 = vx*vx + vy*vy + vz*vz
142 if v2 > vcap*vcap
143 q = vcap / sqrt(v2)
144 vx *= q; vy *= q; vz *= q
145 end
146 return vx, vy, vz
147 end
148
149 @inline function _tct_emf_at(rho, fx, fy, fz, Bx, By, Bz,
150 i, j, k, N, inv2dx, gamma_drag,
151 pressure_coeff, smallr, vcap)
152 c = _tct_lin3(i, j, k, N)
153 vx, vy, vz = _tct_velocity_at(rho, fx, fy, fz, i, j, k, N, inv2dx,
154 gamma_drag, pressure_coeff, smallr, vcap)
155 bx = Bx[c]; by = By[c]; bz = Bz[c]
156 return vy*bz - vz*by, vz*bx - vx*bz, vx*by - vy*bx
157 end
158
159 @inline function _tct_lap_at(A, i, j, k, N, invdx2)
160 c = _tct_lin3(i, j, k, N)
161 return (A[_tct_lin3(i + 1, j, k, N)] + A[_tct_lin3(i - 1, j, k, N)] +
162 A[_tct_lin3(i, j + 1, k, N)] + A[_tct_lin3(i, j - 1, k, N)] +
163 A[_tct_lin3(i, j, k + 1, N)] + A[_tct_lin3(i, j, k - 1, N)] -
164 typeof(invdx2)(6) * A[c]) * invdx2
165 end
166
167 @inline function _tct_bilap_at(A, i, j, k, N, invdx2)
168 T = typeof(invdx2)
169 c = A[_tct_lin3(i, j, k, N)]
170 axis1 = A[_tct_lin3(i + 1, j, k, N)] + A[_tct_lin3(i - 1, j, k, N)] +
171 A[_tct_lin3(i, j + 1, k, N)] + A[_tct_lin3(i, j - 1, k, N)] +
172 A[_tct_lin3(i, j, k + 1, N)] + A[_tct_lin3(i, j, k - 1, N)]
173 axis2 = A[_tct_lin3(i + 2, j, k, N)] + A[_tct_lin3(i - 2, j, k, N)] +
174 A[_tct_lin3(i, j + 2, k, N)] + A[_tct_lin3(i, j - 2, k, N)] +
175 A[_tct_lin3(i, j, k + 2, N)] + A[_tct_lin3(i, j, k - 2, N)]
176 diag = A[_tct_lin3(i + 1, j + 1, k, N)] + A[_tct_lin3(i + 1, j - 1, k, N)] +
177 A[_tct_lin3(i - 1, j + 1, k, N)] + A[_tct_lin3(i - 1, j - 1, k, N)] +
178 A[_tct_lin3(i + 1, j, k + 1, N)] + A[_tct_lin3(i + 1, j, k - 1, N)] +
179 A[_tct_lin3(i - 1, j, k + 1, N)] + A[_tct_lin3(i - 1, j, k - 1, N)] +
180 A[_tct_lin3(i, j + 1, k + 1, N)] + A[_tct_lin3(i, j + 1, k - 1, N)] +
181 A[_tct_lin3(i, j - 1, k + 1, N)] + A[_tct_lin3(i, j - 1, k - 1, N)]
182 return (T(42)*c - T(12)*axis1 + axis2 + T(2)*diag) * invdx2 * invdx2
183 end
184
185 @inline function _tct_rhs_at(rho, fx, fy, fz, Bx, By, Bz,
186 i, j, k, N, inv2dx, invdx2, eta, dissipation_order::Int,
187 gamma_drag, pressure_coeff, smallr, vcap)
188 exyp = _tct_emf_at(rho, fx, fy, fz, Bx, By, Bz, i, j + 1, k, N, inv2dx,
189 gamma_drag, pressure_coeff, smallr, vcap)
190 exym = _tct_emf_at(rho, fx, fy, fz, Bx, By, Bz, i, j - 1, k, N, inv2dx,
191 gamma_drag, pressure_coeff, smallr, vcap)
192 exzp = _tct_emf_at(rho, fx, fy, fz, Bx, By, Bz, i, j, k + 1, N, inv2dx,
193 gamma_drag, pressure_coeff, smallr, vcap)
194 exzm = _tct_emf_at(rho, fx, fy, fz, Bx, By, Bz, i, j, k - 1, N, inv2dx,
195 gamma_drag, pressure_coeff, smallr, vcap)
196 exxp = _tct_emf_at(rho, fx, fy, fz, Bx, By, Bz, i + 1, j, k, N, inv2dx,
197 gamma_drag, pressure_coeff, smallr, vcap)
198 exxm = _tct_emf_at(rho, fx, fy, fz, Bx, By, Bz, i - 1, j, k, N, inv2dx,
199 gamma_drag, pressure_coeff, smallr, vcap)
200 dx = ((exyp[3] - exym[3]) - (exzp[2] - exzm[2])) * inv2dx
201 dy = ((exzp[1] - exzm[1]) - (exxp[3] - exxm[3])) * inv2dx
202 dz = ((exxp[2] - exxm[2]) - (exyp[1] - exym[1])) * inv2dx
203 if dissipation_order == 4
204 return dx - eta * _tct_bilap_at(Bx, i, j, k, N, invdx2),
205 dy - eta * _tct_bilap_at(By, i, j, k, N, invdx2),
206 dz - eta * _tct_bilap_at(Bz, i, j, k, N, invdx2)
207 end
208 return dx + eta * _tct_lap_at(Bx, i, j, k, N, invdx2),
209 dy + eta * _tct_lap_at(By, i, j, k, N, invdx2),
210 dz + eta * _tct_lap_at(Bz, i, j, k, N, invdx2)
211 end
212
213 @kernel function _tct_speed_k!(speed, @Const(rho), @Const(fx), @Const(fy), @Const(fz),
214 N::Int, inv2dx, gamma_drag, pressure_coeff, smallr, vcap)
215 c = @index(Global)
216 @inbounds begin
217 i = (c - 1) % N + 1
218 j = ((c - 1) ÷ N) % N + 1
219 k = (c - 1) ÷ (N*N) + 1
220 vx, vy, vz = _tct_velocity_at(rho, fx, fy, fz, i, j, k, N, inv2dx,
221 gamma_drag, pressure_coeff, smallr, vcap)
222 speed[c] = sqrt(vx*vx + vy*vy + vz*vz)
223 end
224 end
225
226 function terminal_velocity_max!(s::MHDState{T}; gamma_drag::Real,
227 pressure_coeff::Real, vcap::Real=Inf) where {T}
228 N = s.dims[1]
229 fx, fy, fz = s.scratch[2], s.scratch[3], s.scratch[4]
230 speed = s.scratch[1]
231 inv2dx = inv(T(2) * s.dx)
232 terminal_magnetic_force!(s, fx, fy, fz)
233 _tct_speed_k!(s.be)(speed, s.U[1], fx, fy, fz, N, inv2dx,
234 T(gamma_drag), T(pressure_coeff), s.smallr, T(vcap);
235 ndrange=ncells(s))
236 KA.synchronize(s.be)
237 return Float64(maximum(speed))
238 end
239
240 @kernel function _tct_predict_k!(pBx, pBy, pBz,
241 @Const(rho), @Const(fx), @Const(fy), @Const(fz),
242 @Const(Bx), @Const(By), @Const(Bz),
243 N::Int, dt, inv2dx, invdx2, eta, dissipation_order::Int,
244 gamma_drag, pressure_coeff, smallr, vcap)
245 c = @index(Global)
246 @inbounds begin
247 i = (c - 1) % N + 1
248 j = ((c - 1) ÷ N) % N + 1
249 k = (c - 1) ÷ (N*N) + 1
250 x, y, z = _tct_rhs_at(rho, fx, fy, fz, Bx, By, Bz, i, j, k, N,
251 inv2dx, invdx2, eta, dissipation_order,
252 gamma_drag, pressure_coeff, smallr, vcap)
253 pBx[c] = Bx[c] + dt*x
254 pBy[c] = By[c] + dt*y
255 pBz[c] = Bz[c] + dt*z
256 end
257 end
258
259 @kernel function _tct_correct_k!(Bx, By, Bz, E,
260 @Const(rho), @Const(fx), @Const(fy), @Const(fz),
261 @Const(pBx), @Const(pBy), @Const(pBz),
262 N::Int, dt, inv2dx, invdx2, eta, dissipation_order::Int,
263 gamma_drag, pressure_coeff, smallr, smallE, vcap)
264 c = @index(Global)
265 @inbounds begin
266 T = eltype(E)
267 i = (c - 1) % N + 1
268 j = ((c - 1) ÷ N) % N + 1
269 k = (c - 1) ÷ (N*N) + 1
270 x, y, z = _tct_rhs_at(rho, fx, fy, fz, pBx, pBy, pBz, i, j, k, N,
271 inv2dx, invdx2, eta, dissipation_order,
272 gamma_drag, pressure_coeff, smallr, vcap)
273 ox = Bx[c]; oy = By[c]; oz = Bz[c]
274 nx = T(0.5) * (ox + pBx[c] + dt*x)
275 ny = T(0.5) * (oy + pBy[c] + dt*y)
276 nz = T(0.5) * (oz + pBz[c] + dt*z)
277 eb0 = T(0.5) * (ox*ox + oy*oy + oz*oz)
278 eb1 = T(0.5) * (nx*nx + ny*ny + nz*nz)
279 Bx[c] = nx; By[c] = ny; Bz[c] = nz
280 E[c] = max(E[c] + eb1 - eb0, smallE)
281 end
282 end
283
284 # ── Face-staggered (constrained-transport) terminal induction ────────────────────────────
285 # For `mhd_solver=:ct` states the same SSPRK2 scheme advances the FACE fields. The cell EMF
286 # e = v x B is computed exactly as above (cell-centred v and B), averaged to the cell EDGES
287 # (each edge EMF is the mean of the four cells sharing the edge), and the faces advance by the
288 # discrete curl of the edge EMFs. The dissipation is the face-component Laplacian (order 2) or
289 # bilaplacian (order 4) with the identical stencils; both commute with the face divergence
290 # on a periodic grid, so div b stays at round-off. Cell-centred B is the face average after every
291 # stage and the total energy follows the magnetic-energy change as in `_tct_correct_k!`.
292 # Workspace (free outside the Godunov step): ct.ec* (cell EMF), ct.e* (edge EMF), ct.mf*
293 # (predicted faces).
294
295 @kernel function _tctc_cell_emf_k!(ecx, ecy, ecz,
296 @Const(rho), @Const(fx), @Const(fy), @Const(fz),
297 @Const(Bx), @Const(By), @Const(Bz),
298 N::Int, inv2dx, gamma_drag, pressure_coeff, smallr, vcap)
299 c = @index(Global)
300 @inbounds begin
301 i = (c - 1) % N + 1
302 j = ((c - 1) ÷ N) % N + 1
303 k = (c - 1) ÷ (N*N) + 1
304 x, y, z = _tct_emf_at(rho, fx, fy, fz, Bx, By, Bz, i, j, k, N, inv2dx,
305 gamma_drag, pressure_coeff, smallr, vcap)
306 ecx[c] = x; ecy[c] = y; ecz[c] = z
307 end
308 end
309
310 # Cell EMF v x B from prescribed cell velocity arrays (radiation terminal closure).
311 @kernel function _tctc_cell_emf_v_k!(ecx, ecy, ecz, @Const(vx), @Const(vy), @Const(vz),
312 @Const(Bx), @Const(By), @Const(Bz), N::Int)
313 c = @index(Global)
314 @inbounds begin
315 i = (c - 1) % N + 1
316 j = ((c - 1) ÷ N) % N + 1
317 k = (c - 1) ÷ (N*N) + 1
318 x, y, z = _rtct_emf_at(vx, vy, vz, Bx, By, Bz, i, j, k, N)
319 ecx[c] = x; ecy[c] = y; ecz[c] = z
320 end
321 end
322
323 # Edge EMF = mean of the four cells sharing the edge. Lower-face convention of ct.jl:
324 # ez at (i-1/2,j-1/2,k), ex at (i,j-1/2,k-1/2), ey at (i-1/2,j,k-1/2).
325 @kernel function _tctc_edge_k!(ex, ey, ez, @Const(ecx), @Const(ecy), @Const(ecz), N::Int)
326 c = @index(Global)
327 @inbounds begin
328 T = eltype(ex)
329 i = (c - 1) % N + 1
330 j = ((c - 1) ÷ N) % N + 1
331 k = (c - 1) ÷ (N*N) + 1
332 c0 = c
333 cim = _tct_lin3(i - 1, j, k, N); cjm = _tct_lin3(i, j - 1, k, N)
334 ckm = _tct_lin3(i, j, k - 1, N)
335 cimjm = _tct_lin3(i - 1, j - 1, k, N)
336 cjmkm = _tct_lin3(i, j - 1, k - 1, N)
337 ckmim = _tct_lin3(i - 1, j, k - 1, N)
338 q = T(0.25)
339 ez[c] = q * (ecz[c0] + ecz[cim] + ecz[cjm] + ecz[cimjm])
340 ex[c] = q * (ecx[c0] + ecx[cjm] + ecx[ckm] + ecx[cjmkm])
341 ey[c] = q * (ecy[c0] + ecy[cim] + ecy[ckm] + ecy[ckmim])
342 end
343 end
344
345 # d/dt of the faces: curl of the edge EMFs (e = v x B) plus face dissipation of `A*`.
346 @inline function _tctc_rate_at(ex, ey, ez, Ax, Ay, Az, i, j, k, N, invdx, invdx2,
347 eta, dissipation_order::Int)
348 c = _tct_lin3(i, j, k, N)
349 ip = _tct_lin3(i + 1, j, k, N); jp = _tct_lin3(i, j + 1, k, N)
350 kp = _tct_lin3(i, j, k + 1, N)
351 dx = ((ez[jp] - ez[c]) - (ey[kp] - ey[c])) * invdx
352 dy = ((ex[kp] - ex[c]) - (ez[ip] - ez[c])) * invdx
353 dz = ((ey[ip] - ey[c]) - (ex[jp] - ex[c])) * invdx
354 if dissipation_order == 4
355 return dx - eta * _tct_bilap_at(Ax, i, j, k, N, invdx2),
356 dy - eta * _tct_bilap_at(Ay, i, j, k, N, invdx2),
357 dz - eta * _tct_bilap_at(Az, i, j, k, N, invdx2)
358 end
359 return dx + eta * _tct_lap_at(Ax, i, j, k, N, invdx2),
360 dy + eta * _tct_lap_at(Ay, i, j, k, N, invdx2),
361 dz + eta * _tct_lap_at(Az, i, j, k, N, invdx2)
362 end
363
364 @kernel function _tctc_predict_k!(pbx, pby, pbz, @Const(bxf), @Const(byf), @Const(bzf),
365 @Const(ex), @Const(ey), @Const(ez),
366 N::Int, dt, invdx, invdx2, eta, dissipation_order::Int)
367 c = @index(Global)
368 @inbounds begin
369 i = (c - 1) % N + 1
370 j = ((c - 1) ÷ N) % N + 1
371 k = (c - 1) ÷ (N*N) + 1
372 x, y, z = _tctc_rate_at(ex, ey, ez, bxf, byf, bzf, i, j, k, N, invdx, invdx2,
373 eta, dissipation_order)
374 pbx[c] = bxf[c] + dt*x
375 pby[c] = byf[c] + dt*y
376 pbz[c] = bzf[c] + dt*z
377 end
378 end
379
380 # SSPRK2 corrector on the faces; edge EMFs belong to the predicted state. `COMP` selects
381 # Kahan-compensated accumulation (mirrors `_rtct_correct_k!`'s compensated branch).
382 @kernel function _tctc_correct_k!(bxf, byf, bzf, rbx, rby, rbz,
383 @Const(pbx), @Const(pby), @Const(pbz),
384 @Const(ex), @Const(ey), @Const(ez),
385 N::Int, dt, invdx, invdx2, eta, dissipation_order::Int,
386 compensated::Val{COMP}) where {COMP}
387 c = @index(Global)
388 @inbounds begin
389 T = eltype(bxf)
390 i = (c - 1) % N + 1
391 j = ((c - 1) ÷ N) % N + 1
392 k = (c - 1) ÷ (N*N) + 1
393 x, y, z = _tctc_rate_at(ex, ey, ez, pbx, pby, pbz, i, j, k, N, invdx, invdx2,
394 eta, dissipation_order)
395 ox = bxf[c]; oy = byf[c]; oz = bzf[c]
396 if COMP
397 bxf[c] = _compensated_update!(rbx, c, ox, T(0.5) * ((pbx[c] - ox) + dt*x), Val(true))
398 byf[c] = _compensated_update!(rby, c, oy, T(0.5) * ((pby[c] - oy) + dt*y), Val(true))
399 bzf[c] = _compensated_update!(rbz, c, oz, T(0.5) * ((pbz[c] - oz) + dt*z), Val(true))
400 else
401 bxf[c] = T(0.5) * (ox + pbx[c] + dt*x)
402 byf[c] = T(0.5) * (oy + pby[c] + dt*y)
403 bzf[c] = T(0.5) * (oz + pbz[c] + dt*z)
404 end
405 end
406 end
407
408 # Cell B <- face average; E += magnetic-energy change (floor smallE), as `_tct_correct_k!`.
409 @kernel function _tctc_cell_energy_k!(Bx, By, Bz, E, @Const(bxf), @Const(byf), @Const(bzf),
410 N::Int, smallE)
411 c = @index(Global)
412 @inbounds begin
413 T = eltype(E)
414 i = (c - 1) % N + 1
415 j = ((c - 1) ÷ N) % N + 1
416 k = (c - 1) ÷ (N*N) + 1
417 nx = T(0.5) * (bxf[c] + bxf[_tct_lin3(i + 1, j, k, N)])
418 ny = T(0.5) * (byf[c] + byf[_tct_lin3(i, j + 1, k, N)])
419 nz = T(0.5) * (bzf[c] + bzf[_tct_lin3(i, j, k + 1, N)])
420 ox = Bx[c]; oy = By[c]; oz = Bz[c]
421 eb0 = T(0.5) * (ox*ox + oy*oy + oz*oz)
422 eb1 = T(0.5) * (nx*nx + ny*ny + nz*nz)
423 Bx[c] = nx; By[c] = ny; Bz[c] = nz
424 E[c] = max(E[c] + eb1 - eb0, smallE)
425 end
426 end
427
428 # Edge EMFs from the cell EMFs held in ct.ec*.
429 function _tctc_edges!(s::MHDState, N::Int)
430 ct = s.ct
431 _tctc_edge_k!(s.be)(ct.ex, ct.ey, ct.ez, ct.ecx, ct.ecy, ct.ecz, N; ndrange=ncells(s))
432 return nothing
433 end
434
435 # Predictor: faces -> ct.mf* (stage state), cell B of the predictor -> scratch[6:8].
436 function _tctc_predict!(s::MHDState{T}, N, dt, invdx, invdx2, eta, order) where {T}
437 ct = s.ct
438 _tctc_predict_k!(s.be)(ct.mfx, ct.mfy, ct.mfz, ct.bxf, ct.byf, ct.bzf,
439 ct.ex, ct.ey, ct.ez, N, dt, invdx, invdx2, eta, order;
440 ndrange=ncells(s))
441 _ct_bcc_k!(s.be, 256)(s.scratch[6], s.scratch[7], s.scratch[8],
442 ct.mfx, ct.mfy, ct.mfz, N, N, N; ndrange=ncells(s))
443 return nothing
444 end
445
446 # Corrector: faces advance, then cell B / E follow.
447 function _tctc_correct!(s::MHDState{T}, N, dt, invdx, invdx2, eta, order, comp::Val) where {T}
448 ct = s.ct
449 _tctc_correct_k!(s.be)(ct.bxf, ct.byf, ct.bzf, ct.rbx, ct.rby, ct.rbz,
450 ct.mfx, ct.mfy, ct.mfz, ct.ex, ct.ey, ct.ez,
451 N, dt, invdx, invdx2, eta, order, comp; ndrange=ncells(s))
452 _tctc_cell_energy_k!(s.be)(s.U[6], s.U[7], s.U[8], s.U[5], ct.bxf, ct.byf, ct.bzf,
453 N, T(1e-30); ndrange=ncells(s))
454 return nothing
455 end
456
457 # One forward-Euler face update from cell EMFs already stored in ct.ec* (used by the
458 # first-order velocity-step terminal paths of the production runner): edge average, faces +=
459 # dt curl(e) (Kahan-compensated through ct.rb*), then cell B <- face average with E += (|B_face|^2 - |B_cell|^2)/2, replacing the
460 # provisional cell-centred B (and its energy) that the caller's cell kernel wrote.
461 @kernel function _tctc_euler_k!(bxf, byf, bzf, rbx, rby, rbz,
462 @Const(ex), @Const(ey), @Const(ez), N::Int, dt, invdx)
463 c = @index(Global)
464 @inbounds begin
465 i = (c - 1) % N + 1
466 j = ((c - 1) ÷ N) % N + 1
467 k = (c - 1) ÷ (N*N) + 1
468 ip = _tct_lin3(i + 1, j, k, N); jp = _tct_lin3(i, j + 1, k, N)
469 kp = _tct_lin3(i, j, k + 1, N)
470 bxf[c] = _compensated_update!(rbx, c, bxf[c], dt * ((ez[jp] - ez[c]) - (ey[kp] - ey[c])) * invdx, Val(true))
471 byf[c] = _compensated_update!(rby, c, byf[c], dt * ((ex[kp] - ex[c]) - (ez[ip] - ez[c])) * invdx, Val(true))
472 bzf[c] = _compensated_update!(rbz, c, bzf[c], dt * ((ey[ip] - ey[c]) - (ex[jp] - ex[c])) * invdx, Val(true))
473 end
474 end
475
476 function terminal_ct_euler_from_cell_emf!(s::MHDState{T}, dt::Real) where {T}
477 uses_ct(s) || error("terminal_ct_euler_from_cell_emf! needs a CT state")
478 N = s.dims[1]; ct = s.ct
479 _tctc_edges!(s, N)
480 _tctc_euler_k!(s.be)(ct.bxf, ct.byf, ct.bzf, ct.rbx, ct.rby, ct.rbz, ct.ex, ct.ey, ct.ez, N, T(dt), inv(s.dx);
481 ndrange=ncells(s))
482 _ct_cell_from_faces_k!(s.be, 256)(s.U[5], s.U[6], s.U[7], s.U[8], s.U[9],
483 ct.bxf, ct.byf, ct.bzf, N, N, N; ndrange=ncells(s))
484 KA.synchronize(s.be)
485 return nothing
486 end
487
488 """
489 terminal_ct_induction!(s, dt; gamma_drag, pressure_coeff, cfl=0.4,
490 dissipation=0.05, dissipation_order=4,
491 vcap=Inf, max_subcycles=256)
492
493 Advance only the terminal-regime induction equation. Density and pressure are
494 held at the caller-provided macro-step state in `s`; Lorentz force and velocity
495 are recomputed at both SSPRK2 stages. Returns `(nsubcycles, vmax, eta)`.
496
497 For `mhd_solver=:ct` states the faces advance (see the face-staggered section above); with
498 `compensate=true` (default, CT only) each face update is Kahan-compensated through the persistent
499 face remainders `s.ct.rbx/rby/rbz`, which keeps the Float32 face divergence at round-off over
500 thousands of sub-ulp updates. GLM states ignore `compensate`.
501 """
502 function terminal_ct_induction!(s::MHDState{T}, dt::Real;
503 gamma_drag::Real, pressure_coeff::Real,
504 cfl::Real=0.4, dissipation::Real=0.05,
505 dissipation_order::Int=4,
506 vcap::Real=Inf, max_subcycles::Int=256,
507 compensate::Bool=true) where {T}
508 N = s.dims[1]
509 all(==(N), s.dims) || error("terminal CT requires a cubic grid")
510 0 < cfl <= 1 || error("terminal CT cfl must be in (0,1]")
511 0 <= dissipation <= 1 || error("terminal CT dissipation must be in [0,1]")
512 dissipation_order in (2, 4) || error("terminal CT dissipation_order must be 2 or 4")
513 fx, fy, fz = s.scratch[2], s.scratch[3], s.scratch[4]
514 inv2dx = inv(T(2) * s.dx)
515 invdx2 = inv(s.dx * s.dx)
516 vmax = terminal_velocity_max!(s; gamma_drag=gamma_drag,
517 pressure_coeff=pressure_coeff, vcap=vcap)
518 nsub = max(1, ceil(Int, Float64(dt) * vmax /
519 max(Float64(cfl) * Float64(s.dx), eps(Float64))))
520 nsub <= max_subcycles || error("terminal CT needs $nsub subcycles (limit $max_subcycles)")
521 dts = T(dt / nsub)
522 eta = dissipation_order == 4 ?
523 T(dissipation / 12) * T(vmax) * s.dx^3 :
524 T(dissipation) * T(vmax) * s.dx
525 if uses_ct(s)
526 _terminal_ct_face_subcycles!(s, nsub, dts, eta, dissipation_order;
527 gamma_drag, pressure_coeff, vcap, compensate)
528 KA.synchronize(s.be)
529 return nsub, vmax, Float64(eta)
530 end
531 for _ in 1:nsub
532 terminal_magnetic_force!(s, fx, fy, fz)
533 _tct_predict_k!(s.be)(s.scratch[6], s.scratch[7], s.scratch[8],
534 s.U[1], fx, fy, fz, s.U[6], s.U[7], s.U[8],
535 N, dts, inv2dx, invdx2, eta, dissipation_order,
536 T(gamma_drag), T(pressure_coeff), s.smallr, T(vcap);
537 ndrange=ncells(s))
538 _tct_magnetic_force_k!(s.be)(fx, fy, fz,
539 s.scratch[6], s.scratch[7], s.scratch[8],
540 N, inv2dx; ndrange=ncells(s))
541 _tct_correct_k!(s.be)(s.U[6], s.U[7], s.U[8], s.U[5],
542 s.U[1], fx, fy, fz,
543 s.scratch[6], s.scratch[7], s.scratch[8],
544 N, dts, inv2dx, invdx2, eta, dissipation_order,
545 T(gamma_drag), T(pressure_coeff), s.smallr,
546 T(1e-30), T(vcap); ndrange=ncells(s))
547 end
548 KA.synchronize(s.be)
549 return nsub, vmax, Float64(eta)
550 end
551
552 # CT variant of the subcycle loop (see the face-staggered section above).
553 function _terminal_ct_face_subcycles!(s::MHDState{T}, nsub, dts, eta, order;
554 gamma_drag, pressure_coeff, vcap, compensate) where {T}
555 N = s.dims[1]; ct = s.ct
556 fx, fy, fz = s.scratch[2], s.scratch[3], s.scratch[4]
557 inv2dx = inv(T(2) * s.dx); invdx = inv(s.dx); invdx2 = inv(s.dx * s.dx)
558 gd = T(gamma_drag); pc = T(pressure_coeff); vc = T(vcap)
559 for _ in 1:nsub
560 terminal_magnetic_force!(s, fx, fy, fz)
561 _tctc_cell_emf_k!(s.be)(ct.ecx, ct.ecy, ct.ecz, s.U[1], fx, fy, fz,
562 s.U[6], s.U[7], s.U[8], N, inv2dx, gd, pc, s.smallr, vc;
563 ndrange=ncells(s))
564 _tctc_edges!(s, N)
565 _tctc_predict!(s, N, dts, invdx, invdx2, eta, order)
566 _tct_magnetic_force_k!(s.be)(fx, fy, fz, s.scratch[6], s.scratch[7], s.scratch[8],
567 N, inv2dx; ndrange=ncells(s))
568 _tctc_cell_emf_k!(s.be)(ct.ecx, ct.ecy, ct.ecz, s.U[1], fx, fy, fz,
569 s.scratch[6], s.scratch[7], s.scratch[8],
570 N, inv2dx, gd, pc, s.smallr, vc; ndrange=ncells(s))
571 _tctc_edges!(s, N)
572 _tctc_correct!(s, N, dts, invdx, invdx2, eta, order, Val(compensate))
573 end
574 return nothing
575 end