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
../_images/rt-5m-gem_elc_frame5.png

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
../_images/rt-5m-gem_elc_frame5_rhovze.png

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
../_images/rt-5m-gem_fieldEnergy.png

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