Fluid example¶

Physics model and initial conditions¶

This simulation models 2d fast magnetic reconnection in a two-fluid (electron-ion) plasma (see [Birn2001] for the classical paper on this topic and [Wang2015] for our earlier application). The electron and ion fluid spcies are evolved using the five-moment model

\begin{split}\begin{aligned} & \frac{\partial\rho_{s}}{\partial t}+\nabla\cdot\left(\rho_{s}\mathbf{u}_{s}\right)=0,\\ & \frac{\partial\left(\rho_{s}\mathbf{u}_{s}\right)}{\partial t}+\nabla p+\nabla\cdot\left(\rho_{s}\mathbf{u}_{s}\mathbf{u}_{s}\right)=n_{s}q_{s}\left(\mathbf{E}+\mathbf{u}_{s}\times\mathbf{B}\right),\\ & \frac{\partial\mathcal{E}_{s}}{\partial t}+\nabla\cdot\left[\mathbf{u}_{s}\left(p_{s}+\mathcal{E}_{s}\right)\right]=n_{s}q_{s}\mathbf{u}_{s}\cdot\mathbf{E},\end{aligned}\end{split}

where $$\rho_{s}$$, $$\mathbf{u}_{s}$$, $$p_{s}$$, and $$\mathcal{E}=\frac{3}{2}p_{s}+\frac{1}{2}\rho_{s}\mathbf{u}_{s}^{2}$$ are the mass density, velocity, thermal pressure, and total energy of the species $$s$$. The electromagnetic fields are evolved and coupled to the plasma species by the usual Maxwell’s equations

\begin{split}\begin{aligned} & \frac{1}{c^{2}}\frac{\partial\mathbf{E}}{\partial t}=\nabla\times\mathbf{B}-\mu_{0}\sum_{s}\mathbf{J}_{s},\\ & \frac{\partial\mathbf{B}}{\partial t}=-\nabla\times\mathbf{E}.\end{aligned}\end{split}

The algorithms are summarized in more detail in Moments App: Multifluid-moment-Maxwell model, as well as in [Hakim2006], [Hakim2008], and [Wang2020].

The simulation input file is adapted from those used in our earlier work presented in [Wang2015], though only a five-moment version is presented here. We employ a 2D simulation domain that is periodic in $$x$$ ($$-L_{x}/2<x<L_{x}/2$$) and is bounded by conducting walls in $$y$$ at $$y=\pm L_{y}/2$$. Here $$L_{x}=100d_{i0}$$, $$L_{y}=50d_{i0}$$, where $$d_{i0}=\sqrt{m_{i}/\mu_{0}n_{0}\left|e\right|^{2}}$$ is the ion inertia length due to $$n_{0}$$. The initial equilibrium is a single Harris sheet where magnetic field and number densities are specified by

$\mathbf{B}_{0}=B_{0}\tanh\left(y/\lambda_{B}\right)\hat{\mathbf{x}},$

and

$n_{e}=n_{i}=n_{0}\text{sech}^{2}\left(y/\lambda_{B}\right)+n_{b},$

respectively. The total current, $$\mathbf{J}_{0}=\nabla\times\mathbf{B}_{0}/\mu_{0}$$, is decomposed according to $$J_{ze0}/J_{zi0}=T_{i0}/T_{e0}$$, where the initial temperatures $$T_{e0}$$ and $$T_{i0}$$ are constant. A sinusoid perturbation is then applied on the magnetic field according to $$\delta\mathbf{B}=\hat{\mathbf{z}}\times\nabla\psi$$, where

$\psi=\delta\psi\cos\left(2\pi x/L_{x}\right)\cos\left(\pi y/L_{y}\right).$

Input file¶

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 -------------------------------------------------------------------------------- -- App dependencies -------------------------------------------------------------------------------- -- Load the App and equations to be used local Moments = require("App.PlasmaOnCartGrid").Moments() local Euler = require "Eq.Euler" -------------------------------------------------------------------------------- -- Preamble -------------------------------------------------------------------------------- -- Basic constants local gasGamma = 5./3. -- gas gamma local elcCharge = -1.0 -- electron charge local ionCharge = 1.0 -- ion charge local elcMass = 1.0 -- electron mass local lightSpeed = 1.0 -- light speed local mu0 = 1.0 -- mu_0 -- Problem-specific characteristic parameters local n0 = 1.0 -- characteristic density variation local nb_n0 = 0.2 -- background density vs. varation density local plasmaBeta0 = 1.0 -- thermal pressure vs. magnetic pressure local lambda_di0 = 0.5 -- transition layer thickness over ion inertial length local Ti0_Te0 = 5.0 -- ion temperature vs. electron temperature local massRatio = 25.0 -- ion mass vs. electron mass local vAe0_lightSpeed = 0.5 -- electron Alfven speed vs. c local perturbationLevel = 0.1 -- relative perturbation magnitude -- Derived parameters for problem setup or for diagnostic information local epsilon0 = 1.0 / lightSpeed^2 / mu0 -- epsilon_0 local ionMass = elcMass * massRatio -- ion mass local vAe0 = lightSpeed * vAe0_lightSpeed -- electron Alfven speed local vAlf0 = vAe0 * math.sqrt(elcMass / ionMass) -- plasma Alfven speed local B0 = vAlf0 * math.sqrt(n0 * ionMass) -- background B field local OmegaCi0 = ionCharge * B0 / ionMass -- ion cyclotron frequency local OmegaPe0 = math.sqrt(n0 * elcCharge^2 / (epsilon0 * elcMass)) -- electron plasma frequency local de0 = lightSpeed / OmegaPe0 -- electron inertial length local OmegaPi0 = math.sqrt(n0 * ionCharge^2 / (epsilon0 * ionMass)) -- ion plasma frequency local di0 = lightSpeed / OmegaPi0 -- ion inertial length local lambda = lambda_di0 * di0 -- transition layer thickness local psi0 = perturbationLevel * B0 -- potential of perturbation B field -- Domain and time control local Lx, Ly = 25.6 * di0, 12.8 * di0 -- domain lengths local Nx, Ny = 128, 64 -- grid size local tEnd = 25.0 / OmegaCi0 -- end of simulation local nFrame = 5 -- number of output frames at t=tEnd -------------------------------------------------------------------------------- -- App construction -------------------------------------------------------------------------------- local momentApp = Moments.App { ----------------------------------------------------------------------------- -- Common ----------------------------------------------------------------------------- logToFile = true, -- will log to rt-5m-gem_0.log tEnd = tEnd, -- stop the simulation at t=tEnd nFrame = nFrame, -- number of output frames at t=tEnd lower = {-Lx/2, -Ly/2}, -- lower limit of domain in each direction upper = {Lx/2, Ly/2}, -- upper limit of domain in each direction cells = {Nx, Ny}, -- number of cells in each direction timeStepper = "fvDimSplit", -- always use fvDimSplit for fluid simulations, -- i.e., finite-volume dimensional-splitting periodicDirs = {1}, -- periodic boundary condition directions ----------------------------------------------------------------------------- -- Species ----------------------------------------------------------------------------- -- electrons elc = Moments.Species { charge = elcCharge, mass = elcMass, -- charge and mass equation = Euler { gasGamma = gasGamma }, -- default equation to be evolved -- using 2nd-order scheme equationInv = Euler { gasGamma = gasGamma, numericalFlux = "lax" }, -- backup equation to be solved to -- retake the time step in case of -- negative density/pressure init = function (t, xn) -- initial condition to set the -- moments of this species local x, y = xn[1], xn[2] local TeFrac = 1.0 / (1.0 + Ti0_Te0) local sech2 = (1.0 / math.cosh(y / lambda))^2 local n = n0 * (sech2 + nb_n0) local Jz = -(B0 / lambda) * sech2 local Ttotal = plasmaBeta0 * (B0 * B0) / 2.0 /n0 local rho_e = n * elcMass local rhovx_e = 0.0 local rhovy_e = 0.0 local rhovz_e = (elcMass / elcCharge) * Jz * TeFrac local thermalEnergy_e = n * Ttotal * TeFrac / (gasGamma-1.0) local kineticEnergy_e = 0.5 * rhovz_e * rhovz_e / rho_e local e_e = thermalEnergy_e + kineticEnergy_e return rho_e, rhovx_e, rhovy_e, rhovz_e, e_e end, evolve = true, -- whether to evolve this species -- in the hyperbolic step; it could -- still be evolved in the source -- update step, though bcy = { Euler.bcWall, Euler.bcWall }, -- boundary conditions along -- non-periodic directions; in this -- case, the y direction }, -- ions ion = Moments.Species { charge = ionCharge, mass = ionMass, equation = Euler { gasGamma = gasGamma }, equationInv = Euler { gasGamma = gasGamma, numericalFlux = "lax" }, init = function (t, xn) local x, y = xn[1], xn[2] local TiFrac = Ti0_Te0 / (1.0 + Ti0_Te0) local sech2 = (1.0 / math.cosh(y / lambda))^2 local n = n0 * (sech2 + nb_n0) local Jz = -(B0 / lambda) * sech2 local Ttotal = plasmaBeta0 * (B0 * B0)/ 2.0 / n0 local rho_i = n*ionMass local rhovx_i = 0.0 local rhovy_i = 0.0 local rhovz_i = (ionMass / ionCharge) * Jz * TiFrac local thermalEnergy_i = n * Ttotal * TiFrac / (gasGamma - 1) local kineticEnergy_i = 0.5 * rhovz_i * rhovz_i / rho_i local e_i = thermalEnergy_i + kineticEnergy_i return rho_i, rhovx_i, rhovy_i, rhovz_i, e_i end, evolve = true, bcy = { Euler.bcWall, Euler.bcWall }, }, ----------------------------------------------------------------------------- -- Electromagnetic field ----------------------------------------------------------------------------- field = Moments.Field { epsilon0 = epsilon0, mu0 = mu0, init = function (t, xn) local x, y = xn[1], xn[2] local pi = math.pi local sin = math.sin local cos = math.cos local Ex, Ey, Ez = 0.0, 0.0, 0.0 local Bx = B0 * math.tanh(y / lambda) local By = 0.0 local Bz = 0.0 local dBx = -psi0 *(pi / Ly) * cos(2 * pi * x / Lx) * sin(pi * y / Ly) local dBy = psi0 * (2 * pi / Lx) * sin(2 * pi * x / Lx) * cos(pi * y / Ly) local dBz = 0.0 Bx = Bx + dBx By = By + dBy Bz = Bz + dBz return Ex, Ey, Ez, Bx, By, Bz end, evolve = true, -- whether to evolve the field in the hyperbolic step; even -- if this is set to false, the E field could still be -- evolved in the source update step bcy = { Moments.Field.bcReflect, Moments.Field.bcReflect }, -- boundary conditions along the nonperiodic directions, -- in this case, the y direction. }, ----------------------------------------------------------------------------- -- Source equations that couple the plasma fluids and EM fields ----------------------------------------------------------------------------- emSource = Moments.CollisionlessEmSource { species = {"elc", "ion"}, -- plasma species involved timeStepper = "direct", -- time-stepper; one of "time-centered", -- "direct", and "exact" }, } -------------------------------------------------------------------------------- -- Run the App -------------------------------------------------------------------------------- momentApp:run() 

Running the simulation¶

The input file rt-5m-gem.lua can be run using the gkyl executable

mpirun -n 4 ~/gkylsoft/gkyl/bin/gkyl inputFiles/rt-5m-gem.lua


assuming gkyl has been installed in the user’s home directory. When running this simulation, a user should see the following output

Fri Sep 18 2020 15:53:11.000000000
Gkyl built with e0fd133198bd+
Gkyl built on Sep 11 2020 17:59:49
Initializing PlasmaOnCartGrid simulation ...
Using CFL number 1
Initializing completed in 0.0723601 sec

Starting main loop of PlasmaOnCartGrid simulation ...

** Time step too large! Aborting this step! ** Time step 1250 too large! Will retake with dt 1
Step 0 at time 0. Time step 1. Completed 0%
0123456789 Step   125 at time 125. Time step 1. Completed 10%
0123456789 Step   250 at time 250. Time step 1. Completed 20%
0123456789 Step   375 at time 375. Time step 1. Completed 30%
0123


Postprocessing¶

The output of this simulation is the following set of files:

• Electron fluid moments: rt-5m-gem_elc_#.bp.
• Ion fluid moments: rt-5m-gem_ion_#.bp.
• Electromagnetic fields: rt-5m-gem_field_#.bp.
• Field energy: rt-5m-gem_fieldEnergy.bp.
• Diagnostic integrated moments: rt-5m-gem_elc_intMom.bp and rt-5m-gem_in_intMom.bp.

Snapshots (frames) are labeled by the number # at the end of the file name, while integrated diagnostics that are computed as a time-series, such as the field energy, are written out as a single file. Since nFrame=1 in the input file, the only frames that are output are 0, corresponding to the initial condition, and 1, corresponding to the end of the simulation.

Plotting all physical quantities in a frame¶

The following command plot all five moment term in the 5th electron output frame. Note that the five moment terms are ordered as $$\rho,\,\rho v_x\, \rho v_y, \rho v_z, \mathcal{E}$$.

pgkyl rt-5m-gem_elc_5.bp plot --fix-aspect


All five moment terms of the electron fluid.

Plotting one physical quantity¶

The following command plot the 3rd moment term in the 5th electron output frame., The 3rd moment is the out-of-plane momentum, which differs from the electron out-of-plane current by an coefficient $$q_e/m_e$$ .

pgkyl rt-5m-gem_elc_5.bp -c 3 plot --fix-aspect


Out-of-plane electron current, $$\rho u_{e,z}$$, in the 5th (last) output frame.

Plotting integrated quantities over time¶

Finally, since we performed this simulation at two different resolutions, and interesting diagnostic to look at is a comparison of integrated quantities between the two simulations. For ease of plotting we have moved the data from the two simulations to two different folders, res1 ($$16^2 \times 16^2$$) and res2 ($$32^2 \times 32^2$$). Here, we are being agnostic on what a user might have named these two different simulations and labeling them ourselves with postgkyl.

pgkyl rt-5m-gem_fieldEnergy.bp select -c 5 plot --logy


Integrated magnetic field energy, $$|B_z|^2$$, plotted as a function of time.

References¶

 [Hakim2006] Hakim, A., Loverich, J., & Shumlak, U. (2006). A high resolution wave propagation scheme for ideal Two-Fluid plasma equations. Journal of Computational Physics, 219, 418–442. https://doi.org/10.1016/j.jcp.2006.03.036
 [Hakim2008] Hakim, A. H. (2008). Extended MHD modeling with the ten-moment equations. Journal of Fusion Energy, 27, 36–43. https://doi.org/10.1007/s10894-007-9116-z
 [Wang2020] Wang, L., Hakim, A. H., Ng, J., Dong, C., & Germaschewski, K. (2020). Exact and locally implicit source term solvers for multifluid-Maxwell systems. Journal of Computational Physics. https://doi.org/10.1016/j.jcp.2020.109510
 [Birn2001] Birn, L., et al. (2001). Geospace Environmental Modeling ({GEM}) magnetic reconnection challenge. Journal of Geophysical Research, 106, 3715.
 [Wang2015] (1, 2) Wang, L., Hakim, A. H. A. H., Bhattacharjee, A., & Germaschewski, K. (2015). Comparison of multi-fluid moment models with particle-in-cell simulations of collisionless magnetic reconnection. Physics of Plasmas, 22(1), 012108. https://doi.org/10.1063/1.4906063