pair plasma, low frequency "Alfven wave", w=k*vA, small k and w;w is well resolved, none of wce and wpe is resolved; k is well resolved, none of rce and de is resolved; small k requires longer wavelength, thus longer box; perpendicular propagation, Ey, vx, Bz nonzero; vy small;

rt-5m-uniform_r70.lua

-- Gkyl ------------------------------------------------------------------------
local Moments = require("App.PlasmaOnCartGrid").Moments
local Euler = require "Eq.Euler"
local Constants = require "Lib.Constants"
local Mpi = require "Comm.Mpi"

local rank = Mpi.Comm_rank(Mpi.COMM_WORLD)
local seed = rank^2 + 3
math.randomseed(seed)

----------------------
-- HELPER FUNCTIONS --
----------------------
local function log(...)
   if rank == 0 then
      print(string.format(...))
   end
end

----------------
-- PARAMETERS --
----------------
local lightSpeed = 1.
local mu0 = 1.
local epsilon0 = 1 / mu0 / (lightSpeed^2)
local gasGamma = 5 / 3

local ionMass = 1.
local ionCharge = 1.

local p0_pmag0 = 1
local wp0_e_wc0_e = 20
local ionMass_eleMass = 1
local p0_i_p0_e = 1

local Bpert_B0 = 1e-2

local eleMass = ionMass / (ionMass_eleMass)
local eleCharge = -ionCharge

local n0 = 1.
local wp0_e = math.sqrt(n0 * eleCharge^2 / eleMass / epsilon0)
local wc0_e = wp0_e / wp0_e_wc0_e
local B0 = wc0_e * eleMass / eleCharge

-- local B0 = 1.
-- local wc0_e = B0 * eleCharge / eleMass
-- local wp0_e = math.abs(wc0_e) * wp0_e_wc0_e
-- local n0 = wp0_e^2 * (eleMass * epsilon0 / eleCharge^2)

local pmag0 = B0^2/mu0/2
local p0 = pmag0 * p0_pmag0
local p0_e = p0 / (p0_i_p0_e + 1)
local p0_i = p0 - p0_e

local Bpert = Bpert_B0 * B0

local wp0_i = math.sqrt(n0 * ionCharge^2 / ionMass / epsilon0)
local d0_i = lightSpeed / wp0_i
local d0_e = lightSpeed / wp0_e

local Lx = d0_i * 1e6
local xlo, xup, Nx =  -Lx/2, Lx/2, 64
local lower = {xlo}
local upper = {xup}
local cells = {Nx}

local wc0_i = B0 * ionCharge / ionMass
local tEnd = 10000 / math.abs(wc0_i)
local nFrame = 1000
local cfl = 1

-- One Alfven mode (v1x, E1y, B1z).
-- In the complex plane,
--    v1y/E1y = 1/Bz0,
--    B1y/E1y = k/w,
local vA0 = math.abs(B0 / math.sqrt(mu0 * n0 * (ionMass + eleMass)))
local kx = 2 * math.pi / Lx
local w = math.sqrt(kx^2 * vA0^2 / (1 + (vA0/lightSpeed)^2))

-- diagnostic parameters
local rho0_e = n0 * eleMass
local rho0_i = n0 * ionMass

local dx = Lx / Nx
local vt0_e = math.sqrt(3 * p0_e / rho0_e)
local rt0_e = vt0_e / wc0_e
local lambda_De = math.sqrt(epsilon0*p0_e/rho0_e^2*(eleMass/eleCharge)^2)
local vA0_e = B0 / math.sqrt(mu0 * rho0_e)
local wp0 = math.sqrt(wp0_e^2 + wp0_i^2)

local maximumDt = 100
local dt = math.min(cfl * dx / lightSpeed, maximumDt)
local dt = cfl * dx / lightSpeed

log("%30s = %g", "lightSpeed", lightSpeed)
log("%30s = %g", "mu0", mu0)
log("%30s = %g", "epsilon0", epsilon0)
log("%30s = %g", "gasGamma", gasGamma)

log("%30s = %g", "eleMass", eleMass)
log("%30s = %g", "eleCharge", eleCharge)
log("%30s = %g", "ionMass", ionMass)
log("%30s = %g", "ionCharge", ionCharge)

log("%30s = %g", "n0", n0)
log("%30s = %g", "p0_e", p0_e)
log("%30s = %g", "p0_i", p0_i)
log("%30s = %g", "B0", B0)

log("%30s = %g", "dx", dx)
log("%30s = %g", "dx/d0_e", dx / d0_e)
log("%30s = %g", "dx/lambda_De", dx / lambda_De)
log("%30s = %g", "dx/rt0_e", dx / rt0_e)
log("%30s = %g", "dx/rt1_e", dx / rt0_e * Bpert_B0)

log("%30s = %g", "dt", dt)
log("%30s = %g", "dt*wp0_e", dt * wp0_e)
log("%30s = %g", "dt*wc0_e", math.abs(dt * wc0_e))
log("%30s = %g", "dt*wc1_e", math.abs(dt * wc0_e * Bpert_B0))

log("%30s = %g", "Lx", Lx)
log("%30s = %g", "Nx", Nx)
log("%30s = %g", "tEnd", tEnd)
log("%30s = %g", "nFrame", nFrame)

log("%30s = %g", "w/wc0_e", w/wc0_e)
log("%30s = %g", "w/wp0_e", w/wp0_e)
log("%30s = %g", "vA0/lightSpeed", vA0/lightSpeed)
log("%30s = %g", "dt*w", math.abs(dt * w))

-----------------------
-- INITIAL CONDITION --
-----------------------

function calcV(x, y, q, m)
   local Epert = Bpert * w / kx
   local Vpert = Epert / B0
   local vx = 0
   local vy = Vpert * math.cos(kx * x)
   local vz = 0
   return vx, vy, vz
end
function calcE(x, y)
   local Epert = Bpert * w / kx
   local Ex = 0
   local Ey = Epert * math.cos(kx * x)
   local Ez = 0
   return Ex, Ey, Ez
end
function calcB(x, y)
   local Bx = 0
   local By = 0
   local Bz = B0 + Bpert * math.cos(kx * x)
   print(Bz)
   return Bx, By, Bz
end
 
---------
-- APP --
---------
momentApp = Moments.App {
   logToFile = true,

   lower = lower,
   upper = upper,
   cells = cells,
   tEnd = tEnd,
   nFrame = nFrame,
   timeStepper = "fvDimSplit",
--   cfl = cfl,
   maximumDt = maximumDt,

   -- decomposition for configuration space
   decompCuts = {1}, -- cuts in each configuration direction
   useShared = false, -- if to use shared memory

   -- boundary conditions for configuration space
   periodicDirs = {1}, -- periodic directions

   ele = Moments.Species {
      charge = eleCharge, mass = eleMass,
      equation = Euler { gasGamma = gasGamma },
      equationInv = Euler { gasGamma = gasGamma, numericalFlux = "lax" },
      forceInv = true,
      init = function (t, xn)
         local x, y = xn[1], xn[2]
         local rho_e = rho0_e
         local vx, vy, vz = calcV(x, y, eleCharge, eleMass)
         local p_e = p0_e

         local rhovx_e = rho_e * vx
         local rhovy_e = rho_e * vy
         local rhovz_e = rho_e * vz
         local u_e = p_e / (gasGamma - 1)
            + 0.5 * (rhovx_e^2 + rhovy_e^2 + rhovz_e^2) / rho_e
         return rho_e, rhovx_e, rhovy_e, rhovz_e, u_e
      end,
      evolve = false,
   },

   ion = Moments.Species {
      charge = ionCharge, mass = ionMass,
      equation = Euler { gasGamma = gasGamma },
      equationInv = Euler { gasGamma = gasGamma, numericalFlux = "lax" },
      forceInv = true,
      init = function (t, xn)
         local x, y = xn[1], xn[2]
         local rho_i = rho0_i
         local vx, vy, vz = calcV(x, y, ionCharge, ionMass)
         local p_i = p0_i

         local rhovx_i = rho_i * vx
         local rhovy_i = rho_i * vy
         local rhovz_i = rho_i * vz
         local u_i = p_i / (gasGamma - 1)
            + 0.5 * (rhovx_i^2 + rhovy_i^2 + rhovz_i^2) / rho_i
         return rho_i, rhovx_i, rhovy_i, rhovz_i, u_i
      end,
      evolve = false,
   },

   field = Moments.Field {
      epsilon0 = epsilon0, mu0 = mu0,
      limiter = "zero",
      -- limiter = "no-limiter",
      init = function (t, xn)
         local x, y = xn[1], xn[2]
         local Bx, By, Bz = calcB(x, y)

         local Ex, Ey, Ez = calcE(x, y)
         -- local vx, vy, vz = ???
         -- local Ex = -vy * Bz + vz * By
         -- local Ey = -vz * Bx + vx * Bz
         -- local Ez = -vx * By + vy * Bx
         return Ex, Ey, Ez, Bx, By, Bz
      end,
      evolve = true,
   },

   emSource = Moments.CollisionlessEmSource {
      species = {"ele", "ion"},
      evolve = {true, true},
      timeStepper = "direct",
      -- timeStepper = "exact",
   },
}

momentApp:run()