-- Gkyl --------------------------------------------------------------
-- 1x1v simulation of ionization in simple test problem
-- with spatially constant fluid moments and periodic BCs.
----------------------------------------------------------------------
local Plasma = require("App.PlasmaOnCartGrid").VlasovMaxwell()
local Constants = require "Lib.Constants"

-- SI units.
local eV = Constants.ELEMENTARY_CHARGE
local eps0, mu0 = Constants.EPSILON0, 1.257e-6

local qe, qi = -eV, eV
local me, mi = 9.109383e-31, 1.6726218e-27
local n0 = 1.0e17
local ne, ni, nn = n0, n0, 0.01*n0
local vde, vdi = 0.0, 0.0
local Te, Ti = 20*eV, 1*eV

local vte, vti = math.sqrt(Te/me), math.sqrt(Ti/mi)
local uB = math.sqrt(Te/mi)
local omega_pe = math.sqrt((ne * qe^2)/(eps0*me))
local lambda_D = math.sqrt((eps0 * Te)/(ne * qe^2))

-- Initialization function.
local function maxwellian(n, vd, vth, v)
   return n / math.sqrt(2*math.pi*vth*vth) * 
      math.exp(-(v-vd)^2/(2*vth*vth))
end

plasmaApp = Plasma.App {
   logToFile = false,

   tEnd = 1/omega_pe,       -- End time.
   nFrame = 1,              -- Number of output frames.
   lower = {0.0},            -- Configuration space lower left.
   upper = {128.0*lambda_D}, -- Configuration space upper right.
   cells = {64},            -- Configuration space cells.
   basis = "serendipity",    -- One of "serendipity" or "maximal-order".
   polyOrder = 3,            -- Polynomial order.
   cflFrac = 1,              -- CFL "fraction". Usually 1.0.
   timeStepper = "rk3",      -- One of "rk2" or "rk3".

   -- 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.
      
   -- Electrons.
   elc = Plasma.Species {
      charge = qe, mass = me,
      -- Velocity space grid.
      lower = {-6.0*vte},
      upper = {6.0*vte},
      cells = {16},
      decompCuts = {1},
      -- Initial conditions.
      init = function (t, xn)
         local x, v = xn[1], xn[2]
         return maxwellian(ne, 0, vte, v)
      end,
      evolve = true,
      diagnostics = { "M0", "M1i", "M2", "vtSq", "intM0", "intM1i", "intM2Flow", "intM2Thermal" },
      ionization = Plasma.Ionization {
         collideWith = {"neut"},        -- species to collide with
      	 electrons = "elc",             -- define name for electron species
      	 neutrals = "neut",             -- define name for neutral species
      	 elemCharge = eV,               -- define elementary charge
      	 elcMass = me,                  -- electron mass
         plasma = "H",                  -- ion species element
      },
   },

      -- Ions
   ion = Plasma.Species {
      charge = qi, mass = mi,
      -- Velocity space grid.
      lower = {-6.0*uB},
      upper = {6.0*uB},
      cells = {16},
      decompCuts = {1},
      -- Initial conditions.
      init = function (t, xn)
         local x, v = xn[1], xn[2]
         return maxwellian(ni, 0, vti, v)
      end,
      evolve = true,
      diagnostics = { "M0", "M1i", "M2", "intM0", "intM1i", "intM2Flow", "intM2Thermal" },
      ionization = Plasma.Ionization {
         collideWith = {"neut"},        -- species to collide with
      	 electrons = "elc",             -- define name for electron species
      	 neutrals = "neut",             -- define name for neutral species
      	 elemCharge = eV,               -- define elementary charge
      	 elcMass = me,                  -- electron mass
         plasma = "H",                  -- ion species element

      }
   },

   neut = Plasma.Species {
      charge = 0.0, mass = mi,
      -- Velocity space grid
      lower = {-6.0*uB},
      upper = {6.0*uB},
      cells = {16},
      decompCuts = {1},
      -- Initial conditions.
      init = function (t, xn)
         local x, v = xn[1], xn[2]
         return maxwellian(nn, 0, vti, v)
      end,
      evolve = true,
      diagnostics = { "M0", "M1i", "M2", "intM0", "intM1i", "intM2Flow", "intM2Thermal" },
      ionization = Plasma.Ionization {
         collideWith = {"elc"},         -- species to collide with
      	 electrons = "elc",             -- define name for electron species
      	 neutrals = "neut",             -- define name for neutral species
      	 elemCharge = eV,               -- define elementary charge
      	 elcMass = me,                  -- electron mass
         plasma = "H",                  -- ion species element  
      }
   },
   
   -- Field solver.
   field = Plasma.Field {
      epsilon0 = eps0, mu0 = mu0,
      init = function (t, xn)
         return 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
      end,
      evolve = true,
   },
}
-- Run application.
plasmaApp:run()