-- Gkyl ------------------------------------------------------------------------
-- Plasma beach problem: http://ammar-hakim.org/sj/je/je8/je8-plasmabeach.html
--
-- 5 moment simulation of a wave propagating in a plasma that attenuates
-- (has a decreasing density) with x.
--------------------------------------------------------------------------------
local Moments = require("App.PlasmaOnCartGrid").Moments()
local Euler = require "Eq.Euler"
local Constants = require "Lib.Constants"
-- Global parameters.
local cfl = 0.4
local gasGamma = 5.0/3.0
local xlower = 0.0
local xupper = 1.0
local nx = 400
-- Compute coordinate of interior last edge.
local dx = (xupper-xlower)/nx
local xLastEdge = xupper-dx
-- Compute drive frequency.
local dx100 = (xupper-xlower)/100
local deltaT = dx100/Constants.SPEED_OF_LIGHT
local driveOmega = Constants.PI/10/deltaT
local tEnd = 5.0e-9
local nFrames = 100
epsilon0 = Constants.EPSILON0
mu0 = Constants.MU0
local function init(x)
local wpdt = 25*(1-x)^5 -- Plasma frequency.
local factor = deltaT^2*Constants.ELEMENTARY_CHARGE^2/(Constants.ELECTRON_MASS*Constants.EPSILON0)
local ne = wpdt^2/factor
local te = 1.0*Constants.EV2KELVIN -- Electron temperature [K].
local pre = ne*Constants.BOLTZMANN_CONSTANT*te
return Constants.ELECTRON_MASS*ne, 0, 0, 0, pre/(gasGamma-1),
Constants.PROTON_MASS*ne, 0, 0, 0, pre/(gasGamma-1),
0, 0, 0, 0, 0, 0, 0, 0
end
momentApp = Moments.App {
logToFile = true,
tEnd = tEnd,
nFrame = nFrames,
lower = {xlower},
upper = {xupper},
cells = {nx},
cfl = cfl,
timeStepper = "fvDimSplit",
maximumDt = 3e-12,
-- Decomposition for configuration space.
decompCuts = {1}, -- Cuts in each configuration direction.
useShared = false, -- If to use shared memory.
periodicDirs = {}, -- Periodic directions.
-- Electrons.
elc = Moments.Species {
charge = -Constants.ELEMENTARY_CHARGE, mass = Constants.ELECTRON_MASS,
equation = Euler { gasGamma = gasGamma },
equationInv = Euler { gasGamma = gasGamma, numericalFlux = "lax" },
forceInv = false,
-- Initial conditions.
init = function (t, xn)
local x = xn[1]
local rho_e, rhovx_e, rhovy_e, rhovz_e, u_e,
rho_i, rhovx_i, rhovy_i, rhovz_i, u_i,
Ex, Ey, Ez, Bx, By, Bz = init(x)
return rho_e, rhovx_e, rhovy_e, rhovz_e, u_e
end,
evolve = true, -- Evolve species?
bcx = { Moments.Species.bcCopy, Moments.Species.bcCopy },
},
field = Moments.Field {
epsilon0 = epsilon0, mu0 = mu0,
init = function (t, xn)
local x = xn[1]
local rho_e, rhovx_e, rhovy_e, rhovz_e, u_e,
rho_i, rhovx_i, rhovy_i, rhovz_i, u_i,
Ex, Ey, Ez, Bx, By, Bz = init(x)
return Ex, Ey, Ez, Bx, By, Bz
end,
evolve = true, -- Evolve field?
bcx = { Moments.Field.bcCopy, Moments.Field.bcCopy },
},
emSource = Moments.CollisionlessEmSource {
species = {"elc"},
timeStepper = "time-centered",
-- Additional source terms.
-- Not enabled here; for demo purposes.
-- Note: Dp are c arrays and use 0-based indices; xc and qbym are
-- Lua arrays and use 1-based indices.
hasAuxSourceFunction = true,
auxSourceFunction = function (self, xc, t, epsilon0, qbym, fDp, emDp, auxSrcDp)
local x = xc[1]
local nFluids = #qbym
local J0 = 1.0e-10 -- Amps/m^3.
-- Auxiliary source for currents.
for s=0, 0 do
auxSrcDp[s*3+0] = 0
auxSrcDp[s*3+1] = 0
auxSrcDp[s*3+2] = 0
end
-- Auxiliary source for E field.
auxSrcDp[nFluids*3+0] = 0
auxSrcDp[nFluids*3+1] = 0
auxSrcDp[nFluids*3+2] = 0
if (x>xLastEdge) then
auxSrcDp[nFluids*3+1] = -J0*math.sin(driveOmega*t)/Constants.EPSILON0
end
end,
},
}
-- Run application.
momentApp:run()