Appearance
adam_flume_mhd_library
ADAM, FLUME pointwise ideal MHD physics shared by the CPU and FNL backends.
Same contract as adam_flume_euler_library: every routine is pure, takes explicit-size dummies and is tagged !$acc routine seq + !$omp declare target, and never branches on the model (issue #41, D-5). Code units (issue #41, D-2): B is rationalised (B_SI / sqrt(mu0)), so the magnetic pressure is |B|^2 / 2 and
E = p/(gamma-1) + rho |u|^2 / 2 + |B|^2 / 2, p = (gamma-1) (E - rho |u|^2 / 2 - |B|^2 / 2), T = p / (rho R)The routines work on the first NV_MHD conservative variables: the GLM variant (NV_MHD_GLM) passes its whole state, whose leading NV_MHD entries are the same (sequence association); psi is not part of the energy (D-3).
Two model variants (issue #41, D-1/D-14): without divergence control (mhd_*, NV_MHD) and with mixed GLM cleaning (mhd_glm_*, NV_MHD_GLM, the cleaning speed ch an explicit argument). Their eigensystems share the 7x7 Roe-Balsara core of Stone et al. (2008, ApJS 178, appendix B) evaluated at the arithmetic face average, a physical state (so the Roe-matrix factors are X = 0, Y = 1): the B_n row is decoupled (speed 0) without cleaning, the (B_n, psi) pair is the 2x2 block of speeds -+c_h with GLM (block-diagonal eigenvectors). The face split and the back-projection follow the Euler library (compute_face_split_fluxes), each variant at its own compile-time size.
Source: src/app/flume/common/adam_flume_mhd_library.F90
Dependencies
Contents
- mhd_conservative_to_auxiliary
- mhd_eigenvalues
- mhd_eigenvectors
- mhd_face_average
- mhd_face_flux_back_projection
- mhd_face_split_fluxes
- mhd_flux
- mhd_glm_eigenvalues
- mhd_glm_eigenvectors
- mhd_glm_face_flux_back_projection
- mhd_glm_face_split_fluxes
- mhd_glm_flux
- mhd_primitive_to_conservative
- mhd_eigenvectors_core
- mhd_frame_indexes
- mhd_wave_eigenvalues
- mhd_fast_speed
- mhd_sum3
Variables
| Name | Type | Attributes | Description |
|---|---|---|---|
EPS_BT | real(kind=R8P) | parameter | Degenerate transverse field: |B_t| <= EPS_BT max(|B|, sqrt(rho) a). |
EPS_FS | real(kind=R8P) | parameter | Degenerate fast-slow pair (triple umbilic): c_f^2 - c_s^2 <= EPS_FS c_f^2. |
Subroutines
mhd_conservative_to_auxiliary
Compute the auxiliary (primitive and derived) variables of a cell from its conservative variables.
qa(IA_H) is the total specific enthalpy including the magnetic pressure, (E + p + |B|^2 / 2) / rho; the magnetic field is copied into qa(IA_BX:IA_BZ) so the eigen-routines read one array.
Attributes: pure
fortran
subroutine mhd_conservative_to_auxiliary(gamma, R, q, qa)Arguments
| Name | Type | Intent | Attributes | Description |
|---|---|---|---|---|
gamma | real(kind=R8P) | in | Specific heats ratio. | |
R | real(kind=R8P) | in | Gas constant. | |
q | real(kind=R8P) | in | Conservative variables. | |
qa | real(kind=R8P) | out | Auxiliary variables. |
Call graph
mhd_eigenvalues
Compute the eigenvalues of the MHD system without divergence control in direction d: (u_n - c_f, u_n - c_a, u_n - c_s, u_n, u_n + c_s, u_n + c_a, u_n + c_f, 0), the last one of the decoupled B_n.
Attributes: pure
fortran
subroutine mhd_eigenvalues(d, qa, lambda)Arguments
| Name | Type | Intent | Attributes | Description |
|---|---|---|---|---|
d | integer(kind=I4P) | in | Direction, 1=x, 2=y, 3=z. | |
qa | real(kind=R8P) | in | Auxiliary variables. | |
lambda | real(kind=R8P) | out | Eigenvalues. |
Call graph
mhd_eigenvectors
Compute the left (rows) and right (columns) eigenvectors of the MHD system without divergence control in direction d at state qa: the 7x7 Roe-Balsara core on (rho, rho u_n, rho u_t1, rho u_t2, E, B_t1, B_t2) and the decoupled B_n (identity row and column).
Attributes: pure
fortran
subroutine mhd_eigenvectors(gamma, d, qa, el, er)Arguments
| Name | Type | Intent | Attributes | Description |
|---|---|---|---|---|
gamma | real(kind=R8P) | in | Specific heats ratio. | |
d | integer(kind=I4P) | in | Direction, 1=x, 2=y, 3=z. | |
qa | real(kind=R8P) | in | Auxiliary variables (e.g. a face average). | |
el | real(kind=R8P) | out | Left eigenvectors, el(k,:) = l_k. | |
er | real(kind=R8P) | out | Right eigenvectors, er(:,k) = r_k. |
Call graph
mhd_face_average
Compute the arithmetic average of two states for the face eigenvectors (issue #41, section 3.3): density, velocity, pressure and magnetic field averaged, speed of sound and total enthalpy recomputed from them, so the average is a physical state and its eigensystem is exact (no Roe-matrix X, Y factors).
Attributes: pure
fortran
subroutine mhd_face_average(gamma, qaL, qaR, avg)Arguments
| Name | Type | Intent | Attributes | Description |
|---|---|---|---|---|
gamma | real(kind=R8P) | in | Specific heats ratio. | |
qaL | real(kind=R8P) | in | Left state auxiliary variables. | |
qaR | real(kind=R8P) | in | Right state auxiliary variables. | |
avg | real(kind=R8P) | out | Average state auxiliary variables. |
Call graph
mhd_face_flux_back_projection
Return the face flux in conservative variables from the reconstructed split fields, F = R (f+ + f-).
Attributes: pure
fortran
subroutine mhd_face_flux_back_projection(is_characteristic, er, vr, flux)Arguments
| Name | Type | Intent | Attributes | Description |
|---|---|---|---|---|
is_characteristic | logical | in | Reconstruction variables: characteristic or conservative. | |
er | real(kind=R8P) | in | Right eigenvectors, er(:,k) = r_k (unused if conservative). | |
vr | real(kind=R8P) | in | Reconstructed split fields at the face. | |
flux | real(kind=R8P) | out | Face flux. |
Call graph
mhd_face_split_fluxes
Project and Lax-Friedrichs-split the stencil of face i+1/2 (MHD without divergence control), ready for the WENO upwind reconstruction: the algorithm and layout of the Euler compute_face_split_fluxes, with the eigenvectors of the arithmetic face average. The B_n row has speed 0 in both variants, so its face flux is exactly zero.
Attributes: pure
fortran
subroutine mhd_face_split_fluxes(gamma, d, S, is_characteristic, qs, qas, fsplit, er)Arguments
| Name | Type | Intent | Attributes | Description |
|---|---|---|---|---|
gamma | real(kind=R8P) | in | Specific heats ratio. | |
d | integer(kind=I4P) | in | Direction, 1=x, 2=y, 3=z. | |
S | integer(kind=I4P) | in | WENO stencil half-width, S <= S_MAX. | |
is_characteristic | logical | in | Characteristic (or conservative) variables. | |
qs | real(kind=R8P) | in | Stencil conservative variables. | |
qas | real(kind=R8P) | in | Stencil auxiliary variables. | |
fsplit | real(kind=R8P) | out | Split fields in the WENO upwind layout. | |
er | real(kind=R8P) | out | Right eigenvectors (identity if conservative). |
Call graph
mhd_flux
Compute the physical flux of the MHD system without divergence control in direction d (issue #41, section 3.1): (rho u_n, rho u u_n + p_t n - B B_n, (E + p_t) u_n - (u.B) B_n, B u_n - u B_n), the B_n component zero (induction dB/dt = curl(u x B) = -div(u B - B u): the flux of B_j along n is u_n B_j - B_n u_j).
Attributes: pure
fortran
subroutine mhd_flux(d, q, qa, f)Arguments
| Name | Type | Intent | Attributes | Description |
|---|---|---|---|---|
d | integer(kind=I4P) | in | Direction, 1=x, 2=y, 3=z. | |
q | real(kind=R8P) | in | Conservative variables. | |
qa | real(kind=R8P) | in | Auxiliary variables. | |
f | real(kind=R8P) | out | Physical flux. |
Call graph
mhd_glm_eigenvalues
Compute the eigenvalues of the MHD system with GLM cleaning in direction d: the 7 MHD waves and the (B_n, psi) pair -c_h, +c_h.
Attributes: pure
fortran
subroutine mhd_glm_eigenvalues(ch, d, qa, lambda)Arguments
| Name | Type | Intent | Attributes | Description |
|---|---|---|---|---|
ch | real(kind=R8P) | in | GLM cleaning speed. | |
d | integer(kind=I4P) | in | Direction, 1=x, 2=y, 3=z. | |
qa | real(kind=R8P) | in | Auxiliary variables. | |
lambda | real(kind=R8P) | out | Eigenvalues. |
Call graph
mhd_glm_eigenvectors
Compute the left (rows) and right (columns) eigenvectors of the MHD system with GLM cleaning in direction d: block-diagonal, the 7x7 Roe-Balsara core and the (B_n, psi) pair of the flux (psi, c_h^2 B_n), r = (1, -+c_h), l = (1, -+1/c_h) / 2 (issue #41, section 3.3).
Attributes: pure
fortran
subroutine mhd_glm_eigenvectors(ch, gamma, d, qa, el, er)Arguments
| Name | Type | Intent | Attributes | Description |
|---|---|---|---|---|
ch | real(kind=R8P) | in | GLM cleaning speed. | |
gamma | real(kind=R8P) | in | Specific heats ratio. | |
d | integer(kind=I4P) | in | Direction, 1=x, 2=y, 3=z. | |
qa | real(kind=R8P) | in | Auxiliary variables (e.g. a face average). | |
el | real(kind=R8P) | out | Left eigenvectors, el(k,:) = l_k. | |
er | real(kind=R8P) | out | Right eigenvectors, er(:,k) = r_k. |
Call graph
mhd_glm_face_flux_back_projection
Return the face flux in conservative variables from the reconstructed split fields, F = R (f+ + f-) (GLM).
Attributes: pure
fortran
subroutine mhd_glm_face_flux_back_projection(is_characteristic, er, vr, flux)Arguments
| Name | Type | Intent | Attributes | Description |
|---|---|---|---|---|
is_characteristic | logical | in | Reconstruction variables: characteristic or conservative. | |
er | real(kind=R8P) | in | Right eigenvectors, er(:,k) = r_k (unused if conservative). | |
vr | real(kind=R8P) | in | Reconstructed split fields at the face. | |
flux | real(kind=R8P) | out | Face flux. |
Call graph
mhd_glm_face_split_fluxes
Project and Lax-Friedrichs-split the stencil of face i+1/2 (MHD with GLM cleaning): as mhd_face_split_fluxes, with the (B_n, psi) pair at speeds -+c_h (conservative variant: one speed, max(|u_n| + c_f, c_h)).
Attributes: pure
fortran
subroutine mhd_glm_face_split_fluxes(ch, gamma, d, S, is_characteristic, qs, qas, fsplit, er)Arguments
| Name | Type | Intent | Attributes | Description |
|---|---|---|---|---|
ch | real(kind=R8P) | in | GLM cleaning speed. | |
gamma | real(kind=R8P) | in | Specific heats ratio. | |
d | integer(kind=I4P) | in | Direction, 1=x, 2=y, 3=z. | |
S | integer(kind=I4P) | in | WENO stencil half-width, S <= S_MAX. | |
is_characteristic | logical | in | Characteristic (or conservative) variables. | |
qs | real(kind=R8P) | in | Stencil conservative variables. | |
qas | real(kind=R8P) | in | Stencil auxiliary variables. | |
fsplit | real(kind=R8P) | out | Split fields in the WENO upwind layout. | |
er | real(kind=R8P) | out | Right eigenvectors (identity if conservative). |
Call graph
mhd_glm_flux
Compute the physical flux of the MHD system with GLM cleaning in direction d: the MHD flux with the B_n component psi and the psi component c_h^2 B_n (mixed GLM, issue #41, section 3.1).
Attributes: pure
fortran
subroutine mhd_glm_flux(ch, d, q, qa, f)Arguments
| Name | Type | Intent | Attributes | Description |
|---|---|---|---|---|
ch | real(kind=R8P) | in | GLM cleaning speed. | |
d | integer(kind=I4P) | in | Direction, 1=x, 2=y, 3=z. | |
q | real(kind=R8P) | in | Conservative variables. | |
qa | real(kind=R8P) | in | Auxiliary variables. | |
f | real(kind=R8P) | out | Physical flux. |
Call graph
mhd_primitive_to_conservative
Compute the conservative variables of a cell from its primitive variables (psi, if any, is set by the caller).
Attributes: pure
fortran
subroutine mhd_primitive_to_conservative(gamma, r, u, v, w, p, bx, by, bz, q)Arguments
| Name | Type | Intent | Attributes | Description |
|---|---|---|---|---|
gamma | real(kind=R8P) | in | Specific heats ratio. | |
r | real(kind=R8P) | in | Density. | |
u | real(kind=R8P) | in | Velocity components. | |
v | real(kind=R8P) | in | Velocity components. | |
w | real(kind=R8P) | in | Velocity components. | |
p | real(kind=R8P) | in | Pressure. | |
bx | real(kind=R8P) | in | Magnetic field components. | |
by | real(kind=R8P) | in | Magnetic field components. | |
bz | real(kind=R8P) | in | Magnetic field components. | |
q | real(kind=R8P) | out | Conservative variables. |
Call graph
mhd_eigenvectors_core
Compute the 7x7 Roe-Balsara eigenvectors of ideal MHD in direction d at the physical state qa (Stone et al. 2008, ApJS 178, appendix B, conservative variables, with the Roe-matrix factors X = 0, Y = 1 of a physical state). Core variables, in the full-state indexes returned in mp: (rho, rho u_n, rho u_t1, rho u_t2, E, B_t1, B_t2), tangents taken cyclically as in the Euler library; B_n enters as a parameter. Waves ordered (u_n - c_f, u_n - c_a, u_n - c_s, u_n, u_n + c_s, u_n + c_a, u_n + c_f); l7(k,:) is the row l_k, r7(:,k) the column r_k, l7 r7 = I.
Degeneracies (named thresholds, issue #41, section 3.3): a transverse field below EPS_BT of the field (or of sqrt(rho) a) takes the direction (beta_t1, beta_t2) = (1, 0); a fast-slow separation c_f^2 - c_s^2 below EPS_FS c_f^2 (triple umbilic) takes (alpha_f, alpha_s) = (1, 0); alpha_f, alpha_s are clamped to [0, 1].
Attributes: pure
fortran
subroutine mhd_eigenvectors_core(gamma, d, qa, l7, r7, mp)Arguments
| Name | Type | Intent | Attributes | Description |
|---|---|---|---|---|
gamma | real(kind=R8P) | in | Specific heats ratio. | |
d | integer(kind=I4P) | in | Direction, 1=x, 2=y, 3=z. | |
qa | real(kind=R8P) | in | Auxiliary variables. | |
l7 | real(kind=R8P) | out | Left eigenvectors, rows. | |
r7 | real(kind=R8P) | out | Right eigenvectors, columns. | |
mp | integer(kind=I4P) | out | Full-state index of each core variable. |
Call graph
mhd_frame_indexes
Return the conservative variables in the frame order of direction d, (rho, rho u_n, rho u_t1, rho u_t2, E, B_n, B_t1, B_t2), tangents cyclic as in mhd_eigenvectors_core: the projections sum in this order, so a problem rotated from x to y or z sums the same terms in the same order (bitwise invariant, issue #41, MV-4).
Attributes: pure
fortran
subroutine mhd_frame_indexes(d, pv)Arguments
| Name | Type | Intent | Attributes | Description |
|---|---|---|---|---|
d | integer(kind=I4P) | in | Direction, 1=x, 2=y, 3=z. | |
pv | integer(kind=I4P) | out | Full-state index of each frame variable. |
Call graph
mhd_wave_eigenvalues
Compute the 7 MHD wave speeds in direction d, (u_n - c_f, u_n - c_a, u_n - c_s, u_n, u_n + c_s, u_n + c_a, u_n + c_f), the fast speed as mhd_fast_speed, the slow one as a^2 b_n^2 / (rho c_f^2) (no cancellation).
Attributes: pure
fortran
subroutine mhd_wave_eigenvalues(d, qa, lambda)Arguments
| Name | Type | Intent | Attributes | Description |
|---|---|---|---|---|
d | integer(kind=I4P) | in | Direction, 1=x, 2=y, 3=z. | |
qa | real(kind=R8P) | in | Auxiliary variables. | |
lambda | real(kind=R8P) | out | Wave speeds. |
Call graph
Functions
mhd_fast_speed
Return the fast magnetosonic speed of a cell along direction d.
cf^2 = (a^2 + b^2 + sqrt((a^2 - b^2)^2 + 4 a^2 bt^2)) / 2 with b^2 = |B|^2 / rho and bt^2 its part transverse to d: the discriminant (a^2 + b^2)^2 - 4 a^2 bn^2 written as a sum of non-negative terms, so round-off cannot make it negative when B is aligned with d.
Attributes: pure
Returns: real(kind=R8P)
fortran
function mhd_fast_speed(d, qa) result(cf)Arguments
| Name | Type | Intent | Attributes | Description |
|---|---|---|---|---|
d | integer(kind=I4P) | in | Direction, 1=x, 2=y, 3=z. | |
qa | real(kind=R8P) | in | Auxiliary variables. |
Call graph
mhd_sum3
Return a + b + c summed in increasing order, so the result does not depend on the order of the arguments.
Why: |u|^2, |B|^2, u.B summed in the fixed x, y, z order round differently when the components are permuted, so a problem rotated from x to y or z (cyclic tangents) drifts at round-off (issue #41, MV-4). Sorting the three terms (branchless min/max) makes the sum a function of the set of terms: rotations are bitwise invariant.
Attributes: pure
Returns: real(kind=R8P)
fortran
function mhd_sum3(a, b, c) result(s)Arguments
| Name | Type | Intent | Attributes | Description |
|---|---|---|---|---|
a | real(kind=R8P) | in | Terms. | |
b | real(kind=R8P) | in | Terms. | |
c | real(kind=R8P) | in | Terms. |
Call graph