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 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
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
---------------------------------------------------------------------------------- App dependencies---------------------------------------------------------------------------------- Load the App and equations to be usedlocalMoments=require("App.PlasmaOnCartGrid").Moments()localEuler=require"Eq.Euler"---------------------------------------------------------------------------------- Preamble---------------------------------------------------------------------------------- Basic constantslocalgasGamma=5./3.-- gas gammalocalelcCharge=-1.0-- electron chargelocalionCharge=1.0-- ion chargelocalelcMass=1.0-- electron mass locallightSpeed=1.0-- light speedlocalmu0=1.0-- mu_0-- Problem-specific characteristic parameterslocaln0=1.0-- characteristic density variationlocalnb_n0=0.2-- background density vs. varation densitylocalplasmaBeta0=1.0-- thermal pressure vs. magnetic pressurelocallambda_di0=0.5-- transition layer thickness over ion inertial lengthlocalTi0_Te0=5.0-- ion temperature vs. electron temperaturelocalmassRatio=25.0-- ion mass vs. electron masslocalvAe0_lightSpeed=0.5-- electron Alfven speed vs. clocalperturbationLevel=0.1-- relative perturbation magnitude-- Derived parameters for problem setup or for diagnostic informationlocalepsilon0=1.0/lightSpeed^2/mu0-- epsilon_0localionMass=elcMass*massRatio-- ion masslocalvAe0=lightSpeed*vAe0_lightSpeed-- electron Alfven speedlocalvAlf0=vAe0*math.sqrt(elcMass/ionMass)-- plasma Alfven speedlocalB0=vAlf0*math.sqrt(n0*ionMass)-- background B fieldlocalOmegaCi0=ionCharge*B0/ionMass-- ion cyclotron frequencylocalOmegaPe0=math.sqrt(n0*elcCharge^2/(epsilon0*elcMass))-- electron plasma frequencylocalde0=lightSpeed/OmegaPe0-- electron inertial lengthlocalOmegaPi0=math.sqrt(n0*ionCharge^2/(epsilon0*ionMass))-- ion plasma frequencylocaldi0=lightSpeed/OmegaPi0-- ion inertial lengthlocallambda=lambda_di0*di0-- transition layer thicknesslocalpsi0=perturbationLevel*B0-- potential of perturbation B field-- Domain and time controllocalLx,Ly=25.6*di0,12.8*di0-- domain lengthslocalNx,Ny=128,64-- grid sizelocaltEnd=25.0/OmegaCi0-- end of simulationlocalnFrame=5-- number of output frames at t=tEnd---------------------------------------------------------------------------------- App construction--------------------------------------------------------------------------------localmomentApp=Moments.App{------------------------------------------------------------------------------- Common-----------------------------------------------------------------------------logToFile=true,-- will log to rt-5m-gem_0.logtEnd=tEnd,-- stop the simulation at t=tEndnFrame=nFrame,-- number of output frames at t=tEndlower={-Lx/2,-Ly/2},-- lower limit of domain in each directionupper={Lx/2,Ly/2},-- upper limit of domain in each directioncells={Nx,Ny},-- number of cells in each directiontimeStepper="fvDimSplit",-- always use fvDimSplit for fluid simulations,-- i.e., finite-volume dimensional-splittingperiodicDirs={1},-- periodic boundary condition directions------------------------------------------------------------------------------- Species------------------------------------------------------------------------------- electronselc=Moments.Species{charge=elcCharge,mass=elcMass,-- charge and massequation=Euler{gasGamma=gasGamma},-- default equation to be evolved-- using 2nd-order schemeequationInv=Euler{gasGamma=gasGamma,numericalFlux="lax"},-- backup equation to be solved to-- retake the time step in case of -- negative density/pressureinit=function(t,xn)-- initial condition to set the-- moments of this specieslocalx,y=xn[1],xn[2]localTeFrac=1.0/(1.0+Ti0_Te0)localsech2=(1.0/math.cosh(y/lambda))^2localn=n0*(sech2+nb_n0)localJz=-(B0/lambda)*sech2localTtotal=plasmaBeta0*(B0*B0)/2.0/n0localrho_e=n*elcMasslocalrhovx_e=0.0localrhovy_e=0.0localrhovz_e=(elcMass/elcCharge)*Jz*TeFraclocalthermalEnergy_e=n*Ttotal*TeFrac/(gasGamma-1.0)localkineticEnergy_e=0.5*rhovz_e*rhovz_e/rho_elocale_e=thermalEnergy_e+kineticEnergy_ereturnrho_e,rhovx_e,rhovy_e,rhovz_e,e_eend,evolve=true,-- whether to evolve this species-- in the hyperbolic step; it could-- still be evolved in the source-- update step, thoughbcy={Euler.bcWall,Euler.bcWall},-- boundary conditions along-- non-periodic directions; in this-- case, the y direction},-- ionsion=Moments.Species{charge=ionCharge,mass=ionMass,equation=Euler{gasGamma=gasGamma},equationInv=Euler{gasGamma=gasGamma,numericalFlux="lax"},init=function(t,xn)localx,y=xn[1],xn[2]localTiFrac=Ti0_Te0/(1.0+Ti0_Te0)localsech2=(1.0/math.cosh(y/lambda))^2localn=n0*(sech2+nb_n0)localJz=-(B0/lambda)*sech2localTtotal=plasmaBeta0*(B0*B0)/2.0/n0localrho_i=n*ionMasslocalrhovx_i=0.0localrhovy_i=0.0localrhovz_i=(ionMass/ionCharge)*Jz*TiFraclocalthermalEnergy_i=n*Ttotal*TiFrac/(gasGamma-1)localkineticEnergy_i=0.5*rhovz_i*rhovz_i/rho_ilocale_i=thermalEnergy_i+kineticEnergy_ireturnrho_i,rhovx_i,rhovy_i,rhovz_i,e_iend,evolve=true,bcy={Euler.bcWall,Euler.bcWall},},------------------------------------------------------------------------------- Electromagnetic field-----------------------------------------------------------------------------field=Moments.Field{epsilon0=epsilon0,mu0=mu0,init=function(t,xn)localx,y=xn[1],xn[2]localpi=math.pilocalsin=math.sinlocalcos=math.coslocalEx,Ey,Ez=0.0,0.0,0.0localBx=B0*math.tanh(y/lambda)localBy=0.0localBz=0.0localdBx=-psi0*(pi/Ly)*cos(2*pi*x/Lx)*sin(pi*y/Ly)localdBy=psi0*(2*pi/Lx)*sin(2*pi*x/Lx)*cos(pi*y/Ly)localdBz=0.0Bx=Bx+dBxBy=By+dByBz=Bz+dBzreturnEx,Ey,Ez,Bx,By,Bzend,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 stepbcy={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 involvedtimeStepper="direct",-- time-stepper; one of "time-centered",-- "direct", and "exact"},}---------------------------------------------------------------------------------- Run the App--------------------------------------------------------------------------------momentApp:run()
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.
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}\).
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\) .
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.
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
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
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