-- Gkyl ------------------------------------------------------------------------
-- 1x1v simulation of two stream instability
--------------------------------------------------------------------------------
-- App dependencies
--------------------------------------------------------------------------------
local Plasma = require("App.PlasmaOnCartGrid").VlasovMaxwell()
--------------------------------------------------------------------------------
-- Preamble
--------------------------------------------------------------------------------
local vth_e = 0.2             -- electron thermal velocity
local drift = 1.0             -- drift velocity
local perturbK = 0.5          -- perturbation wave-number
local perturbMag = 1.0e-6     -- distribution function perturbation

local function maxwellian(v, n, u, vth)   -- Maxwellian in 1x1v
   return 1/math.sqrt(2*math.pi*vth^2)*math.exp(-(v-u)^2/(2*vth^2))
end

plasmaApp = Plasma.App {
   -----------------------------------------------------------------------------
   -- Common
   -----------------------------------------------------------------------------
   logToFile = true,

   tEnd = 50.0,                  -- end time
   nFrame = 100,                 -- number of output frames
   lower = {-math.pi/perturbK},  -- configuration space lower left
   upper = {math.pi/perturbK},   -- configuration space upper right
   cells = {64},                 -- number of configuration space cells
   basis = "serendipity",        -- basis function
   polyOrder = 2,                -- polynomial order
   timeStepper = "rk3",          -- 3rd order Runge-Kutta

   -- MPI decomposition for configuration space
   decompCuts = {1},   -- number of cuts in each configuration direction
   useShared = false,  -- flag for MPI shared memory

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

   -----------------------------------------------------------------------------
   -- Electrons
   -----------------------------------------------------------------------------
   elc = Plasma.Species {
      charge = -1.0, mass = 1.0,
      -- velocity space grid
      lower = {-6.0},
      upper = {6.0},
      cells = {64},
      decompCuts = {1},
      -- initial conditions
      init = function (t, z)
	 local x, vx = z[1], z[2]
	 local f = maxwellian(vx, 0.5, drift, vth_e) +
            maxwellian(vx, 0.5, -drift, vth_e)
	 return (1+perturbMag*math.cos(perturbK*x))*f
      end,
      evolve = true,
      diagnostics = {"M0","M1i","M2ij","M3i","intM0","intM1i","intM2Flow","intM2Thermal"},
   },

   
   -----------------------------------------------------------------------------
   -- Field Solver
   -----------------------------------------------------------------------------
   field = Plasma.Field {
      epsilon0 = 1.0, mu0 = 1.0,
      init = function (t, z)
         local x = z[1]
	 local Ex = -perturbMag*math.sin(perturbK*x)/perturbK
	 return Ex, 0.0, 0.0, 0.0, 0.0, 0.0
      end,
      evolve = true,
   },
}
--------------------------------------------------------------------------------
-- Run the application
--------------------------------------------------------------------------------
plasmaApp:run()