-- Program to solve Maxwell equations

-- CFL number
cfl = 0.75

-- computational domain
grid = Grid.RectCart2D {
   lower = {-1.0, -1.0},
   upper = {1.0, 1.0},
   cells = {151, 151},
}

-- solution
q = DataStruct.Field2D {
   onGrid = grid,
   -- [Ex, Ey, Ez, Bx, By, Bz, phi_e, phi_m]
   numComponents = 8,
   ghost = {2, 2},
}
qX = DataStruct.Field2D { -- to store results from X sweep
   onGrid = grid,
   -- [Ex, Ey, Ez, Bx, By, Bz, phi_e, phi_m]
   numComponents = 8,
   ghost = {2, 2},
}
-- updated solution
qNew = DataStruct.Field2D {
   onGrid = grid,
   -- [Ex, Ey, Ez, Bx, By, Bz, phi_e, phi_m]
   numComponents = 8,
   ghost = {2, 2},
}

-- create duplicate copy in case we need to take step again
qNewDup = qNew:duplicate()

-- create an alias to in-plane electric field
eFldAlias = qNew:alias(0, 2) -- [Ex, Ey]

-- initial condition to apply
function init(x,y,z)
   local xc, yc = 0.0, 0.0
   local rad2 = (x-xc)^2 + (y-yc)^2
   return 0.0, 0.0, math.exp(-25*rad2), 0.0, 0.0, 0.0, 0.0, 0.0
end

-- apply initial conditions
q:set(init)
qNew:copy(q)

-- write initial conditions
q:write("q_0.h5", 0.0)

-- define equation to solve
maxwellEqn = HyperEquation.PhMaxwell {
   -- speed of light
   lightSpeed = 1.0,
   -- factor for electric field correction potential speed
   elcErrorSpeedFactor = 1.0,
   -- factor for magnetic field correction potential speed
   mgnErrorSpeedFactor = 1.0,
}


-- updater for Euler equations
maxSlvrDir0 = Updater.WavePropagation2D {
   onGrid = grid,
   equation = maxwellEqn,
   -- one of no-limiter, min-mod, superbee, van-leer, monotonized-centered, beam-warming
   limiter = "no-limiter", 
   cfl = cfl,
   cflm = cfl*1.01,
   updateDirections = {0} -- directions to update
}
-- set input/output arrays (these do not change so set it once)
maxSlvrDir0:setIn( {q} )
maxSlvrDir0:setOut( {qX} )

maxSlvrDir1 = Updater.WavePropagation2D {
   onGrid = grid,
   equation = maxwellEqn,
   -- one of no-limiter, min-mod, superbee, van-leer, monotonized-centered, beam-warming
   limiter = "no-limiter", 
   cfl = cfl,
   cflm = cfl*1.01,
   updateDirections = {1} -- directions to update
}
maxSlvrDir1:setIn( {qX} )
maxSlvrDir1:setOut( {qNew} )

-- boundary condition to apply
bcElc = BoundaryCondition.ZeroTangent { components = {0, 1, 2} }
bcMgn = BoundaryCondition.ZeroNormal { components = {3, 4, 5} }
potBc = BoundaryCondition.Copy { components = {6, 7}, fact = {-1, -1} }

-- create boundary condition object
function createBc(myDir, myEdge)
   local bc = Updater.Bc2D {
      onGrid = grid,
      -- boundary conditions to apply
      boundaryConditions = {bcElc, bcMgn, potBc},
      -- direction to apply
      dir = myDir,
      -- edge to apply on
      edge = myEdge,
   }
   return bc
end

-- create updaters to apply boundary conditions
bcLeft = createBc(0, "lower")
bcRight = createBc(0, "upper")
bcBottom = createBc(1, "lower")
bcTop = createBc(1, "upper")

-- apply boundary conditions
function applyBcDir(dir, fld, tCurr, dt)
   if dir == 0 then
      for i,bc in ipairs({bcLeft, bcRight}) do
	 bc:setOut( {fld} )
	 bc:setCurrTime(tCurr)
	 bc:advance(tCurr+dt)
      end
   else
      for i,bc in ipairs({bcTop, bcBottom}) do
	 bc:setOut( {fld} )
	 bc:setCurrTime(tCurr)
	 bc:advance(tCurr+dt)
      end
   end
end

function solveMaxwell(tCurr, dt)
   -- X-direction update
   maxSlvrDir0:setCurrTime(tCurr)
   local s0, dt0 = maxSlvrDir0:advance(tCurr+dt)
   applyBcDir(0, qX, tCurr, dt)

   -- Y-direction update
   maxSlvrDir1:setCurrTime(tCurr)
   local s1, dt1 = maxSlvrDir1:advance(tCurr+dt)
   applyBcDir(1, qNew, tCurr, dt)

   return s1 and s2, math.min(dt0, dt1)
end

-- function to advance solution from tStart to tEnd
function advanceFrame(tStart, tEnd, initDt)
   local step = 1
   local status, dtSuggested
   local tCurr = tStart
   local myDt = initDt
   while true do
      qNewDup:copy(qNew)
      if (tCurr+myDt > tEnd) then -- hit tEnd exactly
	 myDt = tEnd-tCurr
      end

      Lucee.logInfo (string.format("  Taking step %d at time %g with dt %g", step, tCurr, myDt))
      status, dtSuggested = solveMaxwell(tCurr, myDt)

      if (status == false) then
	 Lucee.logInfo (string.format("  ** Time step %g too large! Will retake with dt %g", myDt, dtSuggested))
	 myDt = dtSuggested
	 qNew:copy(qNewDup)
      else
	 q:copy(qNew)
	 tCurr = tCurr + myDt
	 step = step + 1
	 -- check if done
	 if (tCurr >= tEnd) then
	    break
	 end
      end
   end

   return dtSuggested
end

dtSuggested = 100.0 -- initial suggested time-step
-- parameters to control time-stepping
tStart = 0.0
tEnd = 3.0

nFrames = 60
tFrame = (tEnd-tStart)/nFrames -- time between frames

tCurr = tStart
for frame = 1, nFrames do
   Lucee.logInfo (string.format("-- Advancing solution from %g to %g", tCurr, tCurr+tFrame))
   dtSuggested = advanceFrame(tCurr, tCurr+tFrame, dtSuggested)
   q:write( string.format("q_%d.h5", frame), tCurr+tFrame )
   tCurr = tCurr+tFrame
   Lucee.logInfo ("")
end