-- Gkyl ------------------------------------------------------------------------
local Vlasov = (require "App.PlasmaOnCartGrid").VlasovMaxwell()
epsilon0 = 1.0 -- permittivity of free space
mu0 = 1.0 -- pemiability of free space
lightSpeed = 1/math.sqrt(mu0*epsilon0) -- speed of light
Te_Ti = 45.0 -- ratio of electron to ion temperature
n0 = 1.0 -- initial number density
n1 = 2.0 -- initial number density on high side
elcTemp = 1.0e-2 -- electron temperature
elcMass = 1.0 -- electron mass
elcCharge = -1.0 -- electron charge
ionTemp = elcTemp/Te_Ti -- ion temperature
ionMass = 1836.2 -- ion mass
ionCharge = 1.0 -- ion charge
ionFrac = 0.01
AlTemp = elcTemp/Te_Ti -- Al temperature
AlMass = 49577. -- Al mass
AlCharge = 13.0 -- Al charge
AlFrac = 1. - ionFrac
nfrac = math.abs(elcCharge/(ionFrac*ionCharge + AlFrac*AlCharge)) -- Maintain charge neutrality
-- thermal speeds
cs = math.sqrt(elcTemp/ionMass)
vtElc = math.sqrt(elcTemp/elcMass)
vtIon = math.sqrt(ionTemp/ionMass)
vtAl = math.sqrt(AlTemp/ionMass)
-- plasma frequency and Debye length
wpe = math.sqrt(elcCharge^2*n0/(epsilon0*elcMass))
wpi = math.sqrt(ionCharge^2*n0*ionFrac*nfrac/(epsilon0*ionMass))
wpAl = math.sqrt(AlCharge^2*n0*AlFrac*nfrac/(epsilon0*AlMass))
lambdaD = vtElc/wpe
-- domain size and simulation time
LX = 200*lambdaD
vlasovApp = Vlasov.App {
logToFile = true,
tEnd = 1720.0/wpe, -- end time
nFrame = 40, -- number of output frames
lower = {0.0}, -- configuration space lower left
upper = {LX}, -- configuration space upper right
cells = {512}, -- configuration space cells
basis = "serendipity", -- one of "serendipity" or "maximal-order"
polyOrder = 2, -- polynomial order
timeStepper = "rk3", -- one of "rk2" or "rk3"
-- decomposition for configuration space
decompCuts = {256}, -- cuts in each configuration direction
useShared = false, -- if to use shared memory
-- boundary conditions for configuration space
periodicDirs = {1}, -- periodic directions
-- electrons
elc = Vlasov.Species {
charge = elcCharge, mass = elcMass,
-- velocity space grid
lower = {-8.0*vtElc},
upper = {8.0*vtElc},
cells = {64},
decompCuts = {1},
-- initial conditions
init = function (t, xn)
local x, v = xn[1], xn[2]
local sloc1 = 0.25*LX
local sloc2 = 0.75*LX
local fv = n1/math.sqrt(2*math.pi*elcTemp/elcMass)*math.exp(-elcMass*v^2/(2.0*elcTemp))
if x>sloc1 and x<sloc2 then
fv = n0/math.sqrt(2*math.pi*elcTemp/elcMass)*math.exp(-elcMass*v^2/(2.0*elcTemp))
end
return fv
end,
evolve = true, -- evolve species?
diagnostics = { "M0", "M1i", "M2" }
},
-- protons
ion = Vlasov.Species {
charge = ionCharge, mass = ionMass,
-- velocity space grid
lower = {-32.0*vtIon},
upper = {32.0*vtIon},
cells = {256},
decompCuts = {1},
-- initial conditions
init = function (t, xn)
local x, v = xn[1], xn[2]
local sloc1 = 0.25*LX
local sloc2 = 0.75*LX
local fv = n1*nfrac*ionFrac/math.sqrt(2*math.pi*ionTemp/ionMass)*math.exp(-ionMass*v^2/(2.0*ionTemp))
if x>sloc1 and x<sloc2 then
fv = n0*nfrac*ionFrac/math.sqrt(2*math.pi*ionTemp/ionMass)*math.exp(-ionMass*v^2/(2.0*ionTemp))
end
return fv
end,
evolve = true, -- evolve species?
diagnostics = { "M0", "M1i", "M2" }
},
-- Aluminum
Al = Vlasov.Species {
charge = AlCharge, mass = AlMass,
-- velocity space grid
lower = {-64.0*vtAl},
upper = {64.0*vtAl},
cells = {512},
decompCuts = {1},
-- initial conditions
init = function (t, xn)
local x, v = xn[1], xn[2]
local sloc1 = 0.25*LX
local sloc2 = 0.75*LX
local fv = n1*nfrac*AlFrac/math.sqrt(2*math.pi*AlTemp/AlMass)*math.exp(-AlMass*v^2/(2.0*AlTemp))
if x>sloc1 and x<sloc2 then
fv = n0*nfrac*AlFrac/math.sqrt(2*math.pi*AlTemp/AlMass)*math.exp(-AlMass*v^2/(2.0*AlTemp))
end
return fv
end,
evolve = true, -- evolve species?
diagnostics = { "M0", "M1i", "M2" }
},
-- field solver
field = Vlasov.Field {
epsilon0 = 1.0, mu0 = 1.0,
init = function (t, xn)
return 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
end,
evolve = true, -- evolve field?
},
}
-- run application
vlasovApp:run()