Fluid example¶
Contents
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
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
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
and
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
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_elc_#.bp
. - Electromagnetic fields:
rt-5m-gem_field_#.bp
. - Field energy:
rt-5m-gem_fieldEnergy.bp
. - Diagnostic integrated moments:
rt-5m-gem_elc_intMom.bp
andrt-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
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
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
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 |