Swift generator for three-dimensional magnetohydrodynamic turbulence

Magnetohydrodynamic turbulence is central to laboratory and astrophysical plasmas, and is invoked for interpreting many observed scalings. Verifying predicted scaling law behavior requires extreme-resolution direct numerical simulations (DNS), with needed computing resources excluding systematic parameter surveys. We here present an analytic generator of realistically looking turbulent magnetic fields, that computes three-dimensional (3D) $O\left({1000}^3\right)$ solenoidal vector fields in minutes to hours on desktop computers. Our model is inspired by recent developments in 3D incompressible fluid turbulence theory, where a Gaussian white noise vector subjected to a nonlinear transformation results in an intermittent, multifractal random field. Our $B\times C$ model has only few parameters that have clear geometric interpretations. We directly compare a (costly) DNS with a swiftly $B\times C$-generated realization, in terms of its (1) characteristic sheetlike structures of current density, (2) volume-filling aspects across current intensity, (3) power-spectral behavior, (4) probability distribution functions of increments for magnetic field and current density, structure functions, and spectra of exponents, and (5) partial variance of increments. The model even allows to mimic time-evolving magnetic and current density distributions and can be used for synthetic observations on 3D turbulent data cubes.

Figure 1

From deterministic spirals in two dimensions to random sheets in three dimensions. (Top row) 2D setup: Using the deterministic $r_p$ field in the left panel and the $\theta$ field in the middle panel, we construct with (7) the spiral-shaped field $S$ on the right. This $S$ could mimic an isolated eddy. (Bottom row) 3D setup, generalizing the top row: Using the random $R$ field [norm of (3)] in the left panel and the $\theta$ field (8) in the middle panel, we construct similarly the field $S$ with swirling sheets on the right. This $S$ is used in (2) to mimic a distribution of eddies.

Figure 2

Resources required to generate a magnetic field realization with our $B\times C$ code (written in Python). Computing time as a function of resolution is plotted in blue (left $y$ axis), and the required RAM memory in red (right $y$ axis), performed with a 40-logical-cores desktop. Continuous lines are for double precision (float64) calculations, and dashed lines for single precision (float32). Hence, ${1024}^3$ data are generated in about 10 to 25 min, depending on the precision needed.

Figure 3

Figures 3, 4, 5 are visual comparisons of fields generated with a DNS (left column) to fields generated with our $B\times C$ model (right column). In the present figure, (a) in the top row are slices of $B$, the norm of the magnetic field, (b) in the middle row are slices of $j$ in logarithmic scale, the norm of the current density, and (c) in the bottom row are slices of $j_x$, the $x$ component of the current density, which shows some vector information (orientation of $\vec{j}$). From Figs. 3, 4, 5 we conclude that, while $B\times C$ fields are generated using several orders of magnitude fewer resources, they have a similar visual aspect to the DNS.

Figure 4

Continuation of Fig. 3. Isocontours of $j$, the norm of the current density, are shown for values of $60\%$ (top row), $30\%$ (middle row), and $10\%$ (bottom row) of the maximal value. The volume filling and the shape of the structures of the $B\times C$ field at different amplitudes of $j$ matches qualitatively that of the DNS.

Figure 5

Continuation of Figs. 3 and 4. (a) In the top row are 3D vector visualizations, and (b) in the bottom row are 2D vector cuts. While Figs. 3 and 4 show the distribution of the sheets throughout the whole space, here we have zoomed on specific regions to reveal finer details of some clusters of sheetlike structures.

Figure 6

Statistical comparison of the DNS and $B\times C$ fields. In the top panel we plot the power spectrum $P\left(k\right)$ of the norm of the magnetic field, from three orthogonal 2D slices passing through the center of the data cube. Continuous lines correspond to $B\times C$ and dotted lines to the DNS, where in red the data used are from the slice with fixed $x=512$ (the resolution being $N=1024$), in green with fixed $y=512$ and in blue with fixed $z=512$. Red, green, and blue curves of a given data set match because the fields are statistically isotropic. The bottom panel is the same with the current density field. The important point is that the power spectra of the $B\times C$ fields have the characteristic shape of turbulent fields, with a clear power-law inertial range delimited by a large-scale cutoff at small $k$ and a small-scale cutoff at large $k$.

Figure 7

Same as Fig. 6, but where the spectra have been compensated.

Figure 8

Three illustrations of how the power spectrum of the norm of a $B\times C$-generated magnetic field changes when varying some parameters of the model (namely, $h, \eta$, and $L_R$ from top to bottom panels), while keeping the other parameters to their reference values. These examples were made at a resolution $N=512$, and, as in the rest of the paper, the parameter $L_R$ is measured in box-size units and $dx=1/N$. The black arrows suggest how tweaking these parameters may help fitting a given power spectrum.

Figure 9

(a) Top row: PDFs of increments of the norm of the magnetic field generated with our DNS (left column) and that generated with our $B\times C$ model (right column), at lags $\ell =4,7,10,13,17,22,30$. The dotted black curves correspond to unit-variance Gaussian PDFs. As the lag decreases, the curves deviate from Gaussianity, which is characteristic of intermittency. (b) Bottom row: Same plots using the norm of the current density field instead of $B$. The fact that the plots on the left and right columns look like each other indicates that $B\times C$ generated fields have rather realistic statistical properties.

Figure 10

(a) Top row: The first seven structure functions $S_n\left(\ell \right)$ of the norm of the magnetic field generated with our DNS (left column) and that generated with our $B\times C$ model (right column), obtained using the PDFs from Fig. 9. (b) Bottom row: Spectra of exponents resulting from fitting the power-law behaviors of the structure functions shown in the top row. The dashed blue lines correspond to the spectrum of exponents of a nonintermittent field (Kolmogorov scaling), while the dashed black curves correspond to the DNS and $B\times C$ data in the left and right panels, respectively. The pronounced departures of the black curves from the linear laws manifest intermittency, and the similarity between the plots of both columns highlights $B\times C$'s ability to mimic DNS data.

Figure 11

(a) An example of our 1D data path in a $X-Y$ cut used in a PVI quantification, and (b) the locally sampled data series below as function of path parameter $s$.

Figure 12

PVIs computed for $\mathrm{\Delta }s=1,10,50$ (top, middle, bottom). The dashed horizontal lines indicated for each cut (marked in the same color) represent the threshold of PVI $\theta =3\sigma$, calculated separately for each PVI series.

Figure 13

${\mathrm{PVI}}^2$ and the square of the J normalized to its mean value is plotted as a spatial signal with lag $\mathrm{\Delta }s=1$ for the $Z$ cut in (a) the DNS data and (b) the $B\times C$ data.

Figure 14

Kernel density estimate of the joint PDF of the magnitude of the current $\mathbf{J}/{\mathbf{J}}_{rms}$ and the PVI values for the detected PVI events. The PVI signal has been computed on a spatial separation $\mathrm{\Delta }s=1$ and the Pearson correlation coefficient is 0.62 for DNS (a) and 0.71 for $B\times C$ (b).