Snapshot from this documentation build. SHA-256: 262546b8320d0b364bd48b36db8f0498b4b4418312e4fbecd736c9ada0abf53d
1 # ── GLM-MHD device functions (precision-generic, @fastmath) ───────────────────
2 # Ideal MHD + Dedner mixed-GLM divergence cleaning. Conserved U and primitive W:
3 # U = (ρ, ρvx, ρvy, ρvz, E, Bx, By, Bz, ψ)
4 # W = (ρ, vx, vy, vz, p, Bx, By, Bz, ψ)
5 # Reconstruction MUSCL-Hancock + PLM(MonCen); Riemann HLLD (Miyoshi & Kusano 2005)
6 # with an LLF (Rusanov) fallback for low ρ/p; the (Bn,ψ) pair is cleaned per face
7 # (glm_pair). All functions are pure, allocation-free, `@inline`, and parameterised
8 # on `T` so they compile to CPU (T=Float64) and GPU (T=Float32). `@fastmath` lowers
9 # the trace's many min/max/abs to hardware intrinsics (the dominant perf factor).
10
11 @fastmath @inline function cons2prim(c::NTuple{9,T}, γ::T, smallr::T, pfl::T) where {T}
12 rho = max(c[1], smallr); ir = one(T)/rho
13 vx=c[2]*ir; vy=c[3]*ir; vz=c[4]*ir; bx=c[6]; by=c[7]; bz=c[8]
14 ekin = T(0.5)*(c[2]*c[2]+c[3]*c[3]+c[4]*c[4])*ir
15 emag = T(0.5)*(bx*bx+by*by+bz*bz)
16 p = max((γ-one(T))*(c[5]-ekin-emag), pfl)
17 (rho,vx,vy,vz,p,bx,by,bz,c[9])
18 end
19 @inline function prim2cons(q::NTuple{9,T}, γ::T) where {T}
20 rho=q[1];vx=q[2];vy=q[3];vz=q[4];p=q[5];bx=q[6];by=q[7];bz=q[8]
21 E = p/(γ-one(T))+T(0.5)*rho*(vx*vx+vy*vy+vz*vz)+T(0.5)*(bx*bx+by*by+bz*bz)
22 (rho,rho*vx,rho*vy,rho*vz,E,bx,by,bz,q[9])
23 end
24 @fastmath @inline function fast_speed(q::NTuple{9,T}, γ::T, bn::T) where {T}
25 rho=q[1]; c2=γ*q[5]/rho; b2=(q[6]*q[6]+q[7]*q[7]+q[8]*q[8])/rho; d2=T(0.5)*(b2+c2)
26 sqrt(d2 + sqrt(max(d2*d2 - c2*bn*bn/rho, zero(T))))
27 end
28 @inline function phys_flux_x(q::NTuple{9,T}, γ::T) where {T}
29 rho=q[1];vx=q[2];vy=q[3];vz=q[4];p=q[5];bx=q[6];by=q[7];bz=q[8]
30 b2=bx*bx+by*by+bz*bz; ptot=p+T(0.5)*b2
31 E=p/(γ-one(T))+T(0.5)*rho*(vx*vx+vy*vy+vz*vz)+T(0.5)*b2; vb=vx*bx+vy*by+vz*bz
32 (rho*vx, rho*vx*vx+ptot-bx*bx, rho*vx*vy-bx*by, rho*vx*vz-bx*bz,
33 (E+ptot)*vx-bx*vb, zero(T), vx*by-vy*bx, vx*bz-vz*bx, zero(T))
34 end
35 # Rotate so direction `dir` is the x-normal (cyclic v & B permutation).
36 @inline function rot_to(q::NTuple{9,T}, dir::Int) where {T}
37 dir==1 ? q :
38 dir==2 ? (q[1],q[3],q[4],q[2],q[5],q[7],q[8],q[6],q[9]) :
39 (q[1],q[4],q[2],q[3],q[5],q[8],q[6],q[7],q[9])
40 end
41 @inline function rot_flux_from(f::NTuple{9,T}, dir::Int) where {T}
42 dir==1 ? f :
43 dir==2 ? (f[1],f[4],f[2],f[3],f[5],f[8],f[6],f[7],f[9]) :
44 (f[1],f[3],f[4],f[2],f[5],f[7],f[8],f[6],f[9])
45 end
46 # MonCen slope limiter (slope_type=2); @fastmath ⇒ hardware fmin/fmax, branchless.
47 @fastmath @inline function moncen(dl::T, dr::T) where {T}
48 dc = T(0.5)*(dl+dr); sgn = ifelse(dc>=zero(T), one(T), -one(T))
49 val = sgn*min(T(2)*min(abs(dl),abs(dr)), abs(dc))
50 ifelse(dl*dr <= zero(T), zero(T), val)
51 end
52 @fastmath @inline prim_slope(L::NTuple{9,T},M::NTuple{9,T},R::NTuple{9,T}) where {T} =
53 ntuple(i->T(0.5)*moncen(M[i]-L[i], R[i]-M[i]), 9)
54 @fastmath @inline function minmod(dl::T, dr::T) where {T}
55 sgn=ifelse(dl>=zero(T),one(T),-one(T))
56 ifelse(dl*dr<=zero(T),zero(T),sgn*min(abs(dl),abs(dr)))
57 end
58 @fastmath @inline prim_slope_minmod(L::NTuple{9,T},M::NTuple{9,T},R::NTuple{9,T}) where {T} =
59 ntuple(i->T(0.5)*minmod(M[i]-L[i],R[i]-M[i]),9)
60 @inline padd(q::NTuple{9,T}, s::NTuple{9,T}, a::T) where {T} = ntuple(i->q[i]+a*s[i], 9)
61 @inline dir_flux(q::NTuple{9,T}, dir::Int, γ::T) where {T} = rot_flux_from(phys_flux_x(rot_to(q,dir),γ), dir)
62
63 @inline function _drag_phi1(q::T) where {T}
64 abs(q) < T(1e-4) ? one(T) - q/T(2) + q*q/T(6) : -expm1(-q)/q
65 end
66
67 # Exact constant-force drag propagation of a no-drag Hancock face prediction.
68 # `q0` is the reconstructed face at t^n and `qh` is its no-drag t^(n+1/2)
69 # prediction. Only velocity has a direct Compton source; pressure and magnetic
70 # components retain their hyperbolic midpoint predictions.
71 @inline function drag_predict_face(q0::NTuple{9,T}, qh::NTuple{9,T}, impulse::T,
72 vx0::T, vy0::T, vz0::T) where {T}
73 decay = exp(-impulse)
74 phi1 = _drag_phi1(impulse)
75 drag_predict_face(q0, qh, decay, phi1, vx0, vy0, vz0)
76 end
77
78 @inline function drag_predict_face(q0::NTuple{9,T}, qh::NTuple{9,T}, decay::T,
79 phi1::T, vx0::T, vy0::T, vz0::T) where {T}
80 vx = vx0 + decay*(q0[2]-vx0) + phi1*(qh[2]-q0[2])
81 vy = vy0 + decay*(q0[3]-vy0) + phi1*(qh[3]-q0[3])
82 vz = vz0 + decay*(q0[4]-vz0) + phi1*(qh[4]-q0[4])
83 (qh[1],vx,vy,vz,qh[5],qh[6],qh[7],qh[8],qh[9])
84 end
85
86 # Coefficients of Gamma*integral_0^dt |w(t)|^2 dt for
87 # w(t)=exp(-Gamma*t)w0 + phi1(Gamma*t)*t*a. The polynomial branch avoids
88 # cancellation in f32 when the drag impulse is small.
89 @inline function _drag_work_coefficients(q::T) where {T}
90 if q < T(0.05)
91 q2=q*q; q3=q2*q; q4=q3*q
92 A=q*(one(T)-q+T(2/3)*q2-T(1/3)*q3+T(2/15)*q4)
93 C=q*(one(T)-q+T(7/12)*q2-T(1/4)*q3+T(31/360)*q4)
94 D=q*(T(1/3)-q/T(4)+T(7/60)*q2-q3/T(24)+T(31/2520)*q4)
95 return A,C,D
96 end
97 em1 = -expm1(-q)
98 em2 = -expm1(-T(2)*q)
99 A = T(0.5)*em2
100 C = em1*em1/q
101 D = (q-T(2)*em1+T(0.5)*em2)/(q*q)
102 A,C,D
103 end
104
105 struct DragCoefficients{T}
106 face_decay::T
107 face_phi1::T
108 decay::T
109 phi1::T
110 A::T
111 C::T
112 D::T
113 end
114
115 @inline function _drag_coefficients(q::T) where {T}
116 qh=T(0.5)*q
117 A,C,D=_drag_work_coefficients(q)
118 DragCoefficients(exp(-qh),_drag_phi1(qh),exp(-q),_drag_phi1(q),A,C,D)
119 end
120
121 @inline _no_drag_coefficients(::Type{T}) where {T} =
122 DragCoefficients(one(T),one(T),one(T),one(T),zero(T),zero(T),zero(T))
123
124 # Correct the conservative no-drag trial with the same exact constant-force
125 # propagator used at faces. Total energy loses the analytically integrated work
126 # done by drag; flux exchange remains conservative because the shared HLLD flux
127 # already used drag-aware midpoint states.
128 @inline function drag_correct_conserved(u0::NTuple{9,T}, trial::NTuple{9,T},
129 impulse::T, vx0::T, vy0::T, vz0::T,
130 smallr::T, smallE::T) where {T}
131 decay=exp(-impulse); phi1=_drag_phi1(impulse)
132 A,C,D=_drag_work_coefficients(impulse)
133 drag_correct_conserved(u0,trial,decay,phi1,A,C,D,vx0,vy0,vz0,smallr,smallE)
134 end
135
136 @inline function drag_correct_conserved(u0::NTuple{9,T}, trial::NTuple{9,T},
137 decay::T, phi1::T, A::T, C::T, D::T,
138 vx0::T, vy0::T, vz0::T,
139 smallr::T, smallE::T) where {T}
140 state,qwork,emin=_drag_correct_conserved_parts(u0,trial,decay,phi1,A,C,D,
141 vx0,vy0,vz0,smallr,smallE)
142 enew=max(trial[5]-qwork,emin)
143 return (state[1],state[2],state[3],state[4],enew,
144 state[6],state[7],state[8],state[9])
145 end
146
147 @inline function _drag_correct_conserved_parts(u0::NTuple{9,T}, trial::NTuple{9,T},
148 decay::T, phi1::T, A::T, C::T, D::T,
149 vx0::T, vy0::T, vz0::T,
150 smallr::T, smallE::T) where {T}
151 rn=max(trial[1],smallr); ro=max(u0[1],smallr)
152 vox=u0[2]/ro; voy=u0[3]/ro; voz=u0[4]/ro
153 vtx=trial[2]/rn; vty=trial[3]/rn; vtz=trial[4]/rn
154 dvx=vtx-vox; dvy=vty-voy; dvz=vtz-voz
155 wx=vox-vx0; wy=voy-vy0; wz=voz-vz0
156 vnx=vx0+decay*wx+phi1*dvx
157 vny=vy0+decay*wy+phi1*dvy
158 vnz=vz0+decay*wz+phi1*dvz
159 w2=wx*wx+wy*wy+wz*wz
160 wdv=wx*dvx+wy*dvy+wz*dvz
161 dv2=dvx*dvx+dvy*dvy+dvz*dvz
162 # Work on the gas is Gamma*integral v dot (v-v_target) dt. The second
163 # term vanishes in the CMB frame but keeps moving-target semantics exact.
164 iwx=(one(T)-decay)*wx+(one(T)-phi1)*dvx
165 iwy=(one(T)-decay)*wy+(one(T)-phi1)*dvy
166 iwz=(one(T)-decay)*wz+(one(T)-phi1)*dvz
167 qwork=rn*(A*w2+C*wdv+D*dv2+vx0*iwx+vy0*iwy+vz0*iwz)
168 nmx=rn*vnx; nmy=rn*vny; nmz=rn*vnz
169 ek=T(0.5)*(nmx*nmx+nmy*nmy+nmz*nmz)/rn
170 eb=T(0.5)*(trial[6]*trial[6]+trial[7]*trial[7]+trial[8]*trial[8])
171 state=(rn,nmx,nmy,nmz,trial[5],trial[6],trial[7],trial[8],trial[9])
172 return state,qwork,ek+eb+smallE
173 end
174
175 # MUSCL-Hancock half-step predictor: U^{n+1/2} = U^n - ½dt/dx · Σ_d (F(W+s_d)-F(W-s_d)).
176 @fastmath @inline function hancock(m0::NTuple{9,T}, sx,sy,sz, dtdx::T, γ::T) where {T}
177 u0 = prim2cons(m0,γ); h = T(0.5)*dtdx
178 fxp=dir_flux(padd(m0,sx,one(T)),1,γ); fxm=dir_flux(padd(m0,sx,-one(T)),1,γ)
179 fyp=dir_flux(padd(m0,sy,one(T)),2,γ); fym=dir_flux(padd(m0,sy,-one(T)),2,γ)
180 fzp=dir_flux(padd(m0,sz,one(T)),3,γ); fzm=dir_flux(padd(m0,sz,-one(T)),3,γ)
181 ntuple(i->u0[i]-h*((fxp[i]-fxm[i])+(fyp[i]-fym[i])+(fzp[i]-fzm[i])), 9)
182 end
183 # ── Reconstruction options: PLM (MonCen/minmod) and local PPM (CW84), unified as ──
184 # per-variable face OFFSETS (δL,δR) from the cell mean — face⁻ = m0+δL, face⁺ = m0+δR.
185 # Both use only a ±1 stencil, so each fits the cube's 2-cell halo with no resize.
186 # PLM: δR=+½·moncen, δL=-½·moncen (recovers prim_slope/padd exactly). PPM: the
187 # monotonized local parabola edges (Ustyugov-style local stencil + CW84 monotonize).
188 @inline _ppm_edges_unlim(qm::T,q0::T,qp::T) where {T} =
189 (slope=(qp-qm)*T(0.25); curve=(qm-T(2)*q0+qp)/T(12); (q0-slope+curve, q0+slope+curve))
190 @inline function _ppm_monotonize(qL::T,qa::T,qR::T) where {T}
191 dq=qR-qL; qmid=T(0.5)*(qL+qR); dq2_6=dq*dq/T(6); diff=(qa-qmid)*dq
192 if (qR-qa)*(qa-qL) <= zero(T); return (qa,qa)
193 elseif diff > dq2_6; return (T(3)*qa-T(2)*qR, qR)
194 elseif diff < -dq2_6; return (qL, T(3)*qa-T(2)*qL)
195 end
196 (qL,qR)
197 end
198 @inline function _ppm_edges(qm::T,q0::T,qp::T) where {T}
199 qL,qR=_ppm_edges_unlim(qm,q0,qp); _ppm_monotonize(qL,q0,qR)
200 end
201
202 const RECON_PLM = 0
203 const RECON_PPM = 1
204 const RECON_PLM_MINMOD = 2
205 recon_code(s::Symbol) = s === :plm ? RECON_PLM : s === :ppm ? RECON_PPM :
206 s === :plm_minmod ? RECON_PLM_MINMOD :
207 error("unknown recon :$s (have :plm, :plm_minmod, :ppm)")
208
209 # Per-axis face offsets (δL,δR) for the 9 primitives, given the ∓ neighbours.
210 # `rec` is a `Val` so the kernel specialises to ONE reconstruction — a runtime Int
211 # branch would keep both code paths live and spike GPU register pressure (~1.8× slower).
212 @inline function recon_offsets(mm::NTuple{9,T}, m0::NTuple{9,T}, mp::NTuple{9,T}, ::Val{RECON_PLM}) where {T}
213 s = prim_slope(mm, m0, mp) # ½·moncen per var
214 (ntuple(i->-s[i],9), s)
215 end
216 @inline function recon_offsets(mm::NTuple{9,T},m0::NTuple{9,T},mp::NTuple{9,T},::Val{RECON_PLM_MINMOD}) where {T}
217 s=prim_slope_minmod(mm,m0,mp) # ½·minmod per var
218 (ntuple(i->-s[i],9),s)
219 end
220 @inline function recon_offsets(mm::NTuple{9,T}, m0::NTuple{9,T}, mp::NTuple{9,T}, ::Val{RECON_PPM}) where {T}
221 e = ntuple(i->_ppm_edges(mm[i],m0[i],mp[i]), 9) # (qL,qR) per var
222 (ntuple(i->e[i][1]-m0[i],9), ntuple(i->e[i][2]-m0[i],9))
223 end
224 # MUSCL-Hancock ½-step predictor from per-axis edge offsets (PLM δR=+s,δL=-s ⇒ ≡ hancock).
225 @fastmath @inline function hancock_edges(m0::NTuple{9,T}, δLx,δRx,δLy,δRy,δLz,δRz, dtdx::T, γ::T) where {T}
226 u0=prim2cons(m0,γ); h=T(0.5)*dtdx
227 fxp=dir_flux(ntuple(i->m0[i]+δRx[i],9),1,γ); fxm=dir_flux(ntuple(i->m0[i]+δLx[i],9),1,γ)
228 fyp=dir_flux(ntuple(i->m0[i]+δRy[i],9),2,γ); fym=dir_flux(ntuple(i->m0[i]+δLy[i],9),2,γ)
229 fzp=dir_flux(ntuple(i->m0[i]+δRz[i],9),3,γ); fzm=dir_flux(ntuple(i->m0[i]+δLz[i],9),3,γ)
230 ntuple(i->u0[i]-h*((fxp[i]-fxm[i])+(fyp[i]-fym[i])+(fzp[i]-fzm[i])), 9)
231 end
232
233 # Dedner GLM pair: clean the normal field. Returns (bn*, ψ*).
234 @inline function glm_pair(bnL::T,bnR::T,psiL::T,psiR::T,ch::T) where {T}
235 bns = T(0.5)*(bnL+bnR)-T(0.5)*(psiR-psiL)/ch
236 psis = T(0.5)*(psiL+psiR)-T(0.5)*ch*(bnR-bnL)
237 bns, psis
238 end
239 @fastmath @inline function llf_x(L::NTuple{9,T},R::NTuple{9,T},γ::T,ch::T,fbn::T,fpsi::T) where {T}
240 fL=phys_flux_x(L,γ); fR=phys_flux_x(R,γ); uL=prim2cons(L,γ); uR=prim2cons(R,γ)
241 smax=max(abs(L[2])+fast_speed(L,γ,L[6]), abs(R[2])+fast_speed(R,γ,R[6]), ch)
242 f=ntuple(i->T(0.5)*(fL[i]+fR[i])-T(0.5)*smax*(uR[i]-uL[i]), 9)
243 (f[1],f[2],f[3],f[4],f[5],fbn,f[7],f[8],fpsi)
244 end
245 @fastmath @inline function hll_x(L::NTuple{9,T},R::NTuple{9,T},γ::T,ch::T,bn::T,fbn::T,fpsi::T) where {T}
246 cfL=fast_speed(L,γ,bn); cfR=fast_speed(R,γ,bn)
247 SL=min(min(L[2]-cfL, R[2]-cfR), zero(T)); SR=max(max(L[2]+cfL, R[2]+cfR), zero(T))
248 fL=phys_flux_x(L,γ); fR=phys_flux_x(R,γ); uL=prim2cons(L,γ); uR=prim2cons(R,γ)
249 invd=one(T)/(SR-SL)
250 f=ntuple(i->(SR*fL[i] - SL*fR[i] + SL*SR*(uR[i]-uL[i]))*invd, 9)
251 (f[1],f[2],f[3],f[4],f[5],fbn,f[7],f[8],fpsi)
252 end
253 # Miyoshi-Kusano single-star transverse state (lifted out of hlld_x so it's not a
254 # captured closure — closures can inflate GPU register pressure / hurt occupancy).
255 @fastmath @inline function _hlld_star(d::T,u::T,v::T,w::T,by::T,bz::T,S::T,bn::T,SM::T) where {T}
256 denK = d*(S-u)*(S-SM) - bn*bn
257 if abs(denK) < T(1e-12)
258 return (v, w, by, bz)
259 else
260 fac = bn*(SM-u)/denK
261 return (v-by*fac, w-bz*fac, by*(d*(S-u)*(S-u)-bn*bn)/denK, bz*(d*(S-u)*(S-u)-bn*bn)/denK)
262 end
263 end
264 @fastmath @inline function hlld_x(L::NTuple{9,T},R::NTuple{9,T},γ::T,ch::T,bn::T,fbn::T,fpsi::T) where {T}
265 dL,uL,vL,wL,pL=L[1],L[2],L[3],L[4],L[5]; byL,bzL=L[7],L[8]
266 dR,uR,vR,wR,pR=R[1],R[2],R[3],R[4],R[5]; byR,bzR=R[7],R[8]
267 b2L=bn*bn+byL*byL+bzL*bzL; b2R=bn*bn+byR*byR+bzR*bzR
268 ptL=pL+T(0.5)*b2L; ptR=pR+T(0.5)*b2R
269 EL=pL/(γ-one(T))+T(0.5)*dL*(uL*uL+vL*vL+wL*wL)+T(0.5)*b2L
270 ER=pR/(γ-one(T))+T(0.5)*dR*(uR*uR+vR*vR+wR*wR)+T(0.5)*b2R
271 cfL=fast_speed(L,γ,bn); cfR=fast_speed(R,γ,bn)
272 SL=min(min(uL,uR)-max(cfL,cfR), zero(T)); SR=max(max(uL,uR)+max(cfL,cfR), zero(T))
273 UL=(dL,dL*uL,dL*vL,dL*wL,EL,bn,byL,bzL,zero(T)); UR=(dR,dR*uR,dR*vR,dR*wR,ER,bn,byR,bzR,zero(T))
274 FL=phys_flux_x((dL,uL,vL,wL,pL,bn,byL,bzL,zero(T)),γ); FR=phys_flux_x((dR,uR,vR,wR,pR,bn,byR,bzR,zero(T)),γ)
275 den=(SR-uR)*dR-(SL-uL)*dL
276 SM=((SR-uR)*dR*uR-(SL-uL)*dL*uL-ptR+ptL)/den
277 pts=((SR-uR)*dR*ptL-(SL-uL)*dL*ptR+dL*dR*(SR-uR)*(SL-uL)*(uR-uL))/den
278 dLs=max(dL*(SL-uL)/(SL-SM),T(1e-12)); dRs=max(dR*(SR-uR)/(SR-SM),T(1e-12))
279 sqdLs=sqrt(dLs); sqdRs=sqrt(dRs); SLs=SM-abs(bn)/sqdLs; SRs=SM+abs(bn)/sqdRs
280 vyLs,vzLs,byLs,bzLs = _hlld_star(dL,uL,vL,wL,byL,bzL,SL,bn,SM)
281 vyRs,vzRs,byRs,bzRs = _hlld_star(dR,uR,vR,wR,byR,bzR,SR,bn,SM)
282 vdotbL=uL*bn+vL*byL+wL*bzL; vdotbLs=SM*bn+vyLs*byLs+vzLs*bzLs
283 ELs=((SL-uL)*EL-ptL*uL+pts*SM+bn*(vdotbL-vdotbLs))/(SL-SM)
284 vdotbR=uR*bn+vR*byR+wR*bzR; vdotbRs=SM*bn+vyRs*byRs+vzRs*bzRs
285 ERs=((SR-uR)*ER-ptR*uR+pts*SM+bn*(vdotbR-vdotbRs))/(SR-SM)
286 ULs=(dLs,dLs*SM,dLs*vyLs,dLs*vzLs,ELs,bn,byLs,bzLs,zero(T))
287 URs=(dRs,dRs*SM,dRs*vyRs,dRs*vzRs,ERs,bn,byRs,bzRs,zero(T))
288 sgn = bn>=zero(T) ? one(T) : -one(T); invsum=one(T)/(sqdLs+sqdRs)
289 vyss=(sqdLs*vyLs+sqdRs*vyRs+(byRs-byLs)*sgn)*invsum
290 vzss=(sqdLs*vzLs+sqdRs*vzRs+(bzRs-bzLs)*sgn)*invsum
291 byss=(sqdLs*byRs+sqdRs*byLs+sqdLs*sqdRs*(vyRs-vyLs)*sgn)*invsum
292 bzss=(sqdLs*bzRs+sqdRs*bzLs+sqdLs*sqdRs*(vzRs-vzLs)*sgn)*invsum
293 vdotbss=SM*bn+vyss*byss+vzss*bzss
294 ELss=ELs-sqdLs*(vdotbLs-vdotbss)*sgn; ERss=ERs+sqdRs*(vdotbRs-vdotbss)*sgn
295 ULss=(dLs,dLs*SM,dLs*vyss,dLs*vzss,ELss,bn,byss,bzss,zero(T))
296 URss=(dRs,dRs*SM,dRs*vyss,dRs*vzss,ERss,bn,byss,bzss,zero(T))
297 if SLs>=zero(T)
298 f=ntuple(i->FL[i]+SL*(ULs[i]-UL[i]), 9)
299 elseif SM>=zero(T)
300 f=ntuple(i->FL[i]+SLs*ULss[i]-(SLs-SL)*ULs[i]-SL*UL[i], 9)
301 elseif SRs>=zero(T)
302 f=ntuple(i->FR[i]+SRs*URss[i]-(SRs-SR)*URs[i]-SR*UR[i], 9)
303 else
304 f=ntuple(i->FR[i]+SR*(URs[i]-UR[i]), 9)
305 end
306 (f[1],f[2],f[3],f[4],f[5],fbn,f[7],f[8],fpsi)
307 end
308
309 const RSOLVE_HLLD = 0
310 const RSOLVE_HLL = 1
311 const RSOLVE_LLF = 2
312 riemann_code(s::Symbol) = s === :hlld ? RSOLVE_HLLD : s === :hll ? RSOLVE_HLL :
313 s === :llf ? RSOLVE_LLF :
314 error("unknown MHD riemann :$s (have :hlld, :hll, :llf)")
315
316 # Riemann dispatch in direction `dir`: rotate, clean GLM, solve selected MHD Riemann, rotate back.
317 @fastmath @inline function riemann(Lq::NTuple{9,T},Rq::NTuple{9,T},dir::Int,
318 γ::T,ch::T,smallr::T,pfl::T,llf_dmin::T,llf_pmin::T,rsol::Val{RS}) where {T,RS}
319 L0=rot_to(Lq,dir); R0=rot_to(Rq,dir)
320 L=(max(L0[1],smallr),L0[2],L0[3],L0[4],max(L0[5],pfl),L0[6],L0[7],L0[8],L0[9])
321 R=(max(R0[1],smallr),R0[2],R0[3],R0[4],max(R0[5],pfl),R0[6],R0[7],R0[8],R0[9])
322 bns,psis=glm_pair(L[6],R[6],L[9],R[9],ch); fbn=psis; fpsi=ch*ch*bns
323 Lc=(L[1],L[2],L[3],L[4],L[5],bns,L[7],L[8],L[9]); Rc=(R[1],R[2],R[3],R[4],R[5],bns,R[7],R[8],R[9])
324 uself=(llf_dmin>zero(T) && min(L[1],R[1])<llf_dmin)||(llf_pmin>zero(T) && min(L[5],R[5])<llf_pmin)
325 f = RS == RSOLVE_LLF ? llf_x(Lc,Rc,γ,ch,fbn,fpsi) :
326 RS == RSOLVE_HLL ? hll_x(Lc,Rc,γ,ch,bns,fbn,fpsi) :
327 uself ? llf_x(Lc,Rc,γ,ch,fbn,fpsi) : hlld_x(Lc,Rc,γ,ch,bns,fbn,fpsi)
328 rot_flux_from(f,dir)
329 end
330 # Constrained-transport Riemann dispatch (no GLM): the normal field is the single
331 # face-centred value shared by both states (the caller overwrites L[bn] and R[bn] with
332 # it), so b_n* = that value, the normal-field and ψ fluxes are identically zero, and the
333 # LLF fallback speed carries no cleaning speed ch (passed as zero). The transverse
334 # magnetic fluxes (f[6..8] excluding the normal slot) are the face EMFs used by the CT
335 # edge-EMF average. Same HLLD/HLL/LLF + llf_dmin/llf_pmin switching as `riemann`.
336 @fastmath @inline function riemann_ct(Lq::NTuple{9,T},Rq::NTuple{9,T},dir::Int,
337 γ::T,smallr::T,pfl::T,llf_dmin::T,llf_pmin::T,rsol::Val{RS}) where {T,RS}
338 L0=rot_to(Lq,dir); R0=rot_to(Rq,dir)
339 L=(max(L0[1],smallr),L0[2],L0[3],L0[4],max(L0[5],pfl),L0[6],L0[7],L0[8],zero(T))
340 R=(max(R0[1],smallr),R0[2],R0[3],R0[4],max(R0[5],pfl),L0[6],R0[7],R0[8],zero(T))
341 bns=L[6]
342 z=zero(T)
343 uself=(llf_dmin>z && min(L[1],R[1])<llf_dmin)||(llf_pmin>z && min(L[5],R[5])<llf_pmin)
344 f = RS == RSOLVE_LLF ? llf_x(L,R,γ,z,z,z) :
345 RS == RSOLVE_HLL ? hll_x(L,R,γ,z,bns,z,z) :
346 uself ? llf_x(L,R,γ,z,z,z) : hlld_x(L,R,γ,z,bns,z,z)
347 rot_flux_from(f,dir)
348 end
349 @inline _compensated_update!(residual,idx,x::T,dx::T,
350 ::Val{false}) where {T} = x+dx
351
352 @inline function _compensated_update!(residual,idx,x::T,dx::T,
353 ::Val{true}) where {T}
354 carry=@inbounds residual[idx]
355 increment=dx+carry
356 updated=x+increment
357 @inbounds residual[idx]=increment-(updated-x)
358 updated
359 end
360
361 @inline _compensated_floor_update!(residual,idx,x::T,dx::T,floor::T,
362 ::Val{false}) where {T} = max(x+dx,floor)
363
364 @inline function _compensated_floor_update!(residual,idx,x::T,dx::T,floor::T,
365 tracked::Val{true}) where {T}
366 updated=_compensated_update!(residual,idx,x,dx,tracked)
367 if updated < floor
368 @inbounds residual[idx]=zero(T)
369 return floor
370 end
371 return updated
372 end
373
374 @inline function drag_correct_conserved_compensated!(residual,slot,
375 u0::NTuple{9,T},trial::NTuple{9,T},impulse::T,
376 vx0::T,vy0::T,vz0::T,smallr::T,smallE::T,tracked::Val) where {T}
377 decay=exp(-impulse); phi1=_drag_phi1(impulse)
378 A,C,D=_drag_work_coefficients(impulse)
379 drag_correct_conserved_compensated!(residual,slot,u0,trial,decay,phi1,A,C,D,
380 vx0,vy0,vz0,smallr,smallE,tracked)
381 end
382
383 @inline function drag_correct_conserved_compensated!(residual,slot,
384 u0::NTuple{9,T},trial::NTuple{9,T},decay::T,phi1::T,A::T,C::T,D::T,
385 vx0::T,vy0::T,vz0::T,smallr::T,smallE::T,tracked::Val) where {T}
386 state,qwork,emin=_drag_correct_conserved_parts(u0,trial,decay,phi1,A,C,D,
387 vx0,vy0,vz0,smallr,smallE)
388 enew=_compensated_floor_update!(residual,slot,trial[5],-qwork,emin,tracked)
389 return (state[1],state[2],state[3],state[4],enew,
390 state[6],state[7],state[8],state[9])
391 end