-- Gkyl ------------------------------------------------------------------------
-- 2x2v simulation of instabilities driven by counter-streaming beams of plasma.
-- Ions are stationary, neutralizing background.
-- First 16 wave modes are perturbed and allowed to evolve;
-- many of the modes have comparable growth rates, e.g., two-stream and oblique.
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
-- App dependencies
--------------------------------------------------------------------------------
-- Load the Vlasov-Maxwell App
local Plasma = require("App.PlasmaOnCartGrid").VlasovMaxwell()
-- Pseudo-random number generator from SciLua package for initial conditions.
-- Specific rng is the mrg32k3a (multiple recursive generator) of Ecuyer99.
local prng = require "sci.prng"
local rng = prng.mrg32k3a()
--------------------------------------------------------------------------------
-- Preamble
--------------------------------------------------------------------------------
-- Constants
permitt = 1.0 -- Permittivity of free space
permeab = 1.0 -- Permeability of free space
lightSpeed = 1.0/math.sqrt(permitt*permeab) -- Speed of light
chargeElc = -1.0 -- Electron charge
massElc = 1.0 -- Electron mass
-- Initial conditions
-- 1 = Right-going beam; 2 = Left-going beam.
nElc1 = 0.5
nElc2 = 0.5
ud = 0.1 -- Drift velocity of beams
uxElc1 = 0.0
uyElc1 = ud
uxElc2 = 0.0
uyElc2 = -ud
R = 0.1 -- Ratio of thermal velocity to drift velocity
TElc1 = massElc*(R*ud)^2
TElc2 = massElc*(R*ud)^2
vthElc1 = math.sqrt(TElc1/massElc)
vthElc2 = math.sqrt(TElc2/massElc)
k0_TS = 6.135907273413176 -- Wavenumber of fastest growing two-stream mode
theta = 90.0/180.0*math.pi -- 0 deg is pure Weibel, 90 deg is pure two-stream
kx_TS = k0_TS*math.cos(theta)
ky_TS = k0_TS*math.sin(theta)
k0_Weibel = 2.31012970008316 -- Wavenumber of fastest growing Weibel mode
theta = 0.0/180.0*math.pi -- 0 deg is pure Weibel, 90 deg is pure two-stream
kx_Weibel = k0_Weibel*math.cos(theta)
ky_Weibel = k0_Weibel*math.sin(theta)
kx = k0_Weibel
ky = k0_TS/3.0
perturb_n = 1e-8
-- Perturbing the first 16 wave modes with random amplitudes and phases.
-- Note that loop goes from -N to N to sweep all possible phases.
N=16
P={}
for i=-N,N,1 do
P[i]={}
for j=-N,N,1 do
P[i][j]={}
for k=1,6,1 do
P[i][j][k]=rng:sample()
end
end
end
-- Domain size and number of cells
Lx = 2*math.pi/kx
Ly = 2*math.pi/ky
Nx = 16
Ny = 16
vLimElc = 3*ud -- Maximum velocity in velocity space
NvElc = 16
-- Maxwellian in 2x2v
local function maxwellian2D(n, vx, vy, ux, uy, vth)
local v2 = (vx - ux)^2 + (vy - uy)^2
return n/(2*math.pi*vth^2)*math.exp(-v2/(2*vth^2))
end
plasmaApp = Plasma.App {
--------------------------------------------------------------------------------
-- Common
--------------------------------------------------------------------------------
logToFile = true,
tEnd = 50.0, -- End time
nFrame = 1, -- Number of output frames
lower = {0.0,0.0}, -- Lower boundary of configuration space
upper = {Lx,Ly}, -- Upper boundary of configuration space
cells = {Nx,Ny}, -- Configuration space cells
basis = "serendipity", -- One of "serendipity", "maximal-order", or "tensor"
polyOrder = 2, -- Polynomial order
timeStepper = "rk3s4", -- One of "rk2", "rk3", or "rk3s4"
-- MPI decomposition for configuration space
decompCuts = {1,1}, -- Cuts in each configuration direction
useShared = false, -- shared memory is no longer supported
-- Boundary conditions for configuration space
periodicDirs = {1,2}, -- periodic directions (both x and y)
-- Integrated moment flag, compute integrated quantities 1000 times in simulation
calcIntQuantEvery = 0.001,
--------------------------------------------------------------------------------
-- Electrons
--------------------------------------------------------------------------------
elc = Plasma.Species {
charge = chargeElc, mass = massElc,
-- Velocity space grid
lower = {-vLimElc, -vLimElc},
upper = {vLimElc, vLimElc},
cells = {NvElc, NvElc},
-- Initial conditions
init = function (t, xn)
local x, y, vx, vy = xn[1], xn[2], xn[3], xn[4]
local fv = maxwellian2D(nElc1, vx, vy, uxElc1, uyElc1, vthElc1) +
maxwellian2D(nElc2, vx, vy, uxElc2, uyElc2, vthElc2)
return fv
end,
evolve = true,
diagnostics = {"M0","M1i","M2ij","M3i","intM0","intM1i","intM2Flow","intM2Thermal"},
},
--------------------------------------------------------------------------------
-- Field solver
--------------------------------------------------------------------------------
field = Plasma.Field {
epsilon0 = permitt, mu0 = permeab,
init = function (t, xn)
local x, y = xn[1], xn[2]
local E_x, E_y, B_z = 0.0, 0.0, 0.0
for i=-N,N,1 do
for j=-N,N,1 do
if i~=0 or j~=0 then
E_x = E_x + perturb_n*P[i][j][1]*math.sin(i*kx*x+j*ky*y+2*math.pi*P[i][j][2])
E_y = E_y + perturb_n*P[i][j][3]*math.sin(i*kx*x+j*ky*y+2*math.pi*P[i][j][4])
B_z = B_z + perturb_n*P[i][j][5]*math.sin(i*kx*x+j*ky*y+2*math.pi*P[i][j][6])
end
end
end
return E_x, E_y, 0.0, 0.0, 0.0, B_z
end,
evolve = true,
},
}
--------------------------------------------------------------------------------
-- Run application
--------------------------------------------------------------------------------
plasmaApp:run()