Skip to content

Instantly share code, notes, and snippets.

@liangwang0734

liangwang0734/diff.patch

Last active October 11, 2018 20:57
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save liangwang0734/c0e2542bade8a44946078025271b9717 to your computer and use it in GitHub Desktop.
Save liangwang0734/c0e2542bade8a44946078025271b9717 to your computer and use it in GitHub Desktop.
Download ZIP
gkeyll v1 fake vs real 2d
Raw
diff.patch
--- fake_2d/gkeyll_2d_fake.lua 2018-10-11 16:37:19.000000000 -0400
+++ real_2d/gkeyll_2d_real.lua 2018-10-11 16:39:50.000000000 -0400
@@ -54,10 +54,9 @@
-- domain and grid parameters
xlo, xup, nx = -8. * di0, 8. * di0, 64
ylo, yup, ny = -16. * di0, 16. * di0, 128
- zlo, zup, nz = -0.5 * di0, 0.5 * di0, 1
- lower = {xlo, ylo, zlo}
- upper = {xup, yup, zup}
- cells = {nx, ny, nz}
+ lower = {xlo, ylo}
+ upper = {xup, yup}
+ cells = {nx, ny}
-- computational options
cfl = 0.9
cflm = 1.1 * cfl
@@ -119,8 +118,8 @@
----------------------------
-- DECOMPOSITION AND GRID --
----------------------------
-decomposition = DecompRegionCalc3D.CartGeneral {}
-grid = Grid.RectCart3D {
+decomposition = DecompRegionCalc2D.CartGeneral {}
+grid = Grid.RectCart2D {
lower = lower,
upper = upper,
cells = cells,
@@ -135,22 +134,19 @@
if not numComponents then
numComponents = 18
end
- return DataStruct.Field3D {
+ return DataStruct.Field2D {
onGrid = grid,
numComponents = numComponents,
ghost = {2, 2},
- -- writeGhost = {2, 2},
+ --writeGhost = {2, 2},
}
end
q = createData()
-qNew = createData()
+qX = createData()
+qNew = q
+qY = qNew
qDup = createData()
--- reuse array for dimensional-splitting
-qX = qNew
-qY = q
-qZ = qNew
-
getFields = function(q)
return q:alias(0, 5), q:alias(5, 10), q:alias(10, 18)
end
@@ -158,16 +154,11 @@
elcNew,ionNew,emfNew = getFields(qNew)
elcX,ionX,emfX = getFields(qX)
elcY,ionY,emfY = getFields(qY)
-elcZ,ionZ,emfZ = getFields(qZ)
-------------------------
-- BOUNDARY CONDITIONS --
-------------------------
function applyBc(q, tCurr, tEnd, dir)
- if (dir == 2) then
- q:applyCopyBc(2, "lower")
- q:applyCopyBc(2, "upper")
- end
q:sync()
end
@@ -188,7 +179,7 @@
}
createSlvr = function(equation, input, output, limiter, updateDirections)
- local slvr = Updater.WavePropagation3D {
+ local slvr = Updater.WavePropagation2D {
onGrid = grid,
equation = equation,
-- one of no-limiter, zero, min-mod, superbee,
@@ -211,14 +202,9 @@
ionEqnSlvrY = createSlvr(fluidEqn, ionX, ionY, limiter, {1})
emfEqnSlvrY = createSlvr(emfEqn, emfX, emfY, limiter, {1})
-elcEqnSlvrZ = createSlvr(fluidEqn, elcY, elcZ, limiter, {2})
-ionEqnSlvrZ = createSlvr(fluidEqn, ionY, ionZ, limiter, {2})
-emfEqnSlvrZ = createSlvr(emfEqn, emfY, emfZ, limiter, {2})
-
slvrs = {
{elcEqnSlvrX, ionEqnSlvrX, emfEqnSlvrX},
{elcEqnSlvrY, ionEqnSlvrY, emfEqnSlvrY},
- {elcEqnSlvrZ, ionEqnSlvrZ, emfEqnSlvrZ},
}
elcEqnLaxSlvrX = createSlvr(fluidEqnLax, elc, elcX, "zero", {0})
@@ -228,20 +214,15 @@
elcEqnLaxSlvrY = createSlvr(fluidEqnLax, elcX, elcY, "zero", {1})
ionEqnLaxSlvrY = createSlvr(fluidEqnLax, ionX, ionY, "zero", {1})
emfEqnLaxSlvrY = createSlvr(emfEqn, emfX, emfY, "zero", {1})
-
-elcEqnLaxSlvrZ = createSlvr(fluidEqnLax, elcY, elcZ, "zero", {2})
-ionEqnLaxSlvrZ = createSlvr(fluidEqnLax, ionY, ionZ, "zero", {2})
-emfEqnLaxSlvrZ = createSlvr(emfEqn, emfY, emfZ, "zero", {2})
laxSlvrs = {
{elcEqnLaxSlvrX, ionEqnLaxSlvrX, emfEqnLaxSlvrX},
{elcEqnLaxSlvrY, ionEqnLaxSlvrY, emfEqnLaxSlvrY},
- {elcEqnLaxSlvrZ, ionEqnLaxSlvrZ, emfEqnLaxSlvrZ},
}
-qInList = {q, qX, qY}
-elcOutList = {elcX, elcY, elcZ}
-ionOutList = {ionX, ionY, ionZ}
+qInList = {q, qX}
+elcOutList = {elcX, elcY}
+ionOutList = {ionX, ionY}
----------------------------------
-- HYPERBOLIC EQUATION UPDATERS --
@@ -251,7 +232,7 @@
local myDtSuggested = 1e3 * math.abs(tEnd-tCurr)
local useLaxFlux = false
- for _, dir in ipairs({0, 1, 2}) do
+ for _, dir in ipairs({0, 1}) do
applyBc(qInList[dir + 1], tCurr, tEnd, dir)
for _, slvr in ipairs(slvrs[dir + 1]) do
slvr:setCurrTime(tCurr)
@@ -280,7 +261,7 @@
local myStatus = true
local myDtSuggested = 1e3*math.abs(tEnd-tCurr)
- for _, dir in ipairs({0, 1, 2}) do
+ for _, dir in ipairs({0, 1}) do
applyBc(qInList[dir + 1], tCurr, tEnd, dir)
for _, slvr in ipairs(laxSlvrs[dir + 1]) do
slvr:setCurrTime(tCurr)
@@ -300,7 +281,7 @@
---------------------
-- SOURCE UPDATERS --
---------------------
-srcSlvr = Updater.ImplicitFiveMomentSrc3D {
+srcSlvr = Updater.ImplicitFiveMomentSrc2D {
onGrid = grid,
numFluids = numFluids,
charge = charge,
@@ -453,8 +434,6 @@
if (tCurr > tNextFrame or tCurr >= tEnd) then
log(">>> Writing output %d at t = %g...", frame, tCurr)
- qNew:applyCopyBc(2, "lower")
- qNew:applyCopyBc(2, "upper")
writeData(qNew, frame, tCurr, "q")
frame = frame + 1
tNextFrame = tNextFrame + tFrame
Raw
gkeyll_fake_2d.lua
----------------------
-- HELPER FUNCTIONS --
----------------------
function log(...)
Lucee.logInfo(string.format(...))
end
function HERE()
local info = debug.getinfo(2)
local str = string.format("HERE: %d %s: %s",
info.currentline, info.source, tostring(info.name))
Lucee.logInfo(str)
end
function writeData(q, frame, tCurr, tag)
local time0 = os.time()
local fname = string.format("%s_%d.h5", tag, frame)
q:write(fname, tCurr)
log(" Writing %s took %gs", fname, os.time() - time0)
end
----------------
-- PARAMETERS --
----------------
-- physical constants
gasGamma = 5. / 3.
lightSpeed = 1.
mu0 = 1.
epsilon0 = 1. / mu0 / lightSpeed^2
-- kinetic parameters
numFluids = 2
mi = 1.
me = 1. / 25.
qi = 1.
qe = -1.
mass = {me, mi}
charge = {qe, qi}
massRatio = mi /me
-- problem specific parameters
wpe0_wce0 = 4.
n0 = 1.
betae0 = 0.5
betai0 = 0.5
wpe0 = math.sqrt(n0 * qe^2 / epsilon0 / me)
wce0 = wpe0 / wpe0_wce0
B0 = math.abs(wce0 * me / qe)
pmag0 = B0^2 / 2. / mu0
pe0 = pmag0 * betae0
pi0 = pmag0 * betai0
wpi0 = math.sqrt(n0 * qe^2 / epsilon0 / me)
wci0 = B0 * qi / mi
de0 = lightSpeed / wpe0
di0 = lightSpeed / wpi0
-- domain and grid parameters
xlo, xup, nx = -8. * di0, 8. * di0, 64
ylo, yup, ny = -16. * di0, 16. * di0, 128
zlo, zup, nz = -0.5 * di0, 0.5 * di0, 1
lower = {xlo, ylo, zlo}
upper = {xup, yup, zup}
cells = {nx, ny, nz}
-- computational options
cfl = 0.9
cflm = 1.1 * cfl
limiter = "monotonized-centered"
elcErrorSpeedFactor = 0
mgnErrorSpeedFactor = 1
elcErrorDampFactor = 0
mgnErrorDampFactor = 0
linearSolver = "partialPivLu"
alwaysUseLaxFlux = false
-- I/O control
tEnd = 1 / wci0
tFrame = 0.1 / wci0
Lucee.IsRestarting = false
Lucee.RestartFrame = -1
outputFwd = false
log("%30s = %g", "massRatio", massRatio)
log("%30s = %g", "di0", di0)
log("%30s = %g", "1/wci0", 1./wci0)
log("%30s = %g", "tFrame", tFrame)
-----------------------
-- INITIAL CONDITION --
-----------------------
function init(x, y, z)
z = 0.
local rho = n0 * (1. + 0.1 * math.cos(x / di0))
local vx, vy, vz = 0.1, 0.05, 0.02 * math.sin(x / di0)
local rho_e = rho / (1 + massRatio)
local rho_i = rho - rho_e
local p_e = pe0
local p_i = pi0
local rhovx_e = rho_e * vx
local rhovy_e = rho_e * vy
local rhovz_e = rho_e * vz
local u_e = p_e / (gasGamma - 1.) + 0.5 * (vx^2 + vy^2 + vz^2) * rho_e
local rhovx_i = rho_i * vx
local rhovy_i = rho_i * vy
local rhovz_i = rho_i * vz
local u_i = p_i / (gasGamma - 1.) + 0.5 * (vx^2 + vy^2 + vz^2) * rho_i
local Bx, By, Bz = B0, B0 * 2, B0 * 3
local Ex = - vy * Bz + vz * By
local Ey = - vz * Bx + vx * Bz
local Ez = - vx * By + vy * Bx
return 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, 0., 0.
end
----------------------------
-- DECOMPOSITION AND GRID --
----------------------------
decomposition = DecompRegionCalc3D.CartGeneral {}
grid = Grid.RectCart3D {
lower = lower,
upper = upper,
cells = cells,
decomposition = decomposition,
periodicDirs ={0, 1},
}
----------
-- DATA --
----------
createData = function(numComponents)
if not numComponents then
numComponents = 18
end
return DataStruct.Field3D {
onGrid = grid,
numComponents = numComponents,
ghost = {2, 2},
-- writeGhost = {2, 2},
}
end
q = createData()
qNew = createData()
qDup = createData()
-- reuse array for dimensional-splitting
qX = qNew
qY = q
qZ = qNew
getFields = function(q)
return q:alias(0, 5), q:alias(5, 10), q:alias(10, 18)
end
elc,ion,emf = getFields(q)
elcNew,ionNew,emfNew = getFields(qNew)
elcX,ionX,emfX = getFields(qX)
elcY,ionY,emfY = getFields(qY)
elcZ,ionZ,emfZ = getFields(qZ)
-------------------------
-- BOUNDARY CONDITIONS --
-------------------------
function applyBc(q, tCurr, tEnd, dir)
if (dir == 2) then
q:applyCopyBc(2, "lower")
q:applyCopyBc(2, "upper")
end
q:sync()
end
---------------------------------------
-- HYPERBOLIC EQUATIONS AND SOLVERS --
---------------------------------------
fluidEqn = HyperEquation.Euler {
gasGamma = gasGamma,
}
fluidEqnLax = HyperEquation.Euler {
gasGamma = gasGamma,
numericalFlux = "lax",
}
emfEqn = HyperEquation.PhMaxwell {
lightSpeed = lightSpeed,
elcErrorSpeedFactor = elcErrorSpeedFactor,
mgnErrorSpeedFactor = mgnErrorSpeedFactor,
}
createSlvr = function(equation, input, output, limiter, updateDirections)
local slvr = Updater.WavePropagation3D {
onGrid = grid,
equation = equation,
-- one of no-limiter, zero, min-mod, superbee,
-- van-leer, monotonized-centered, beam-warming
limiter = limiter,
cfl = cfl,
cflm = cflm,
updateDirections = updateDirections,
}
slvr:setIn( {input} )
slvr:setOut( {output} )
return slvr
end
elcEqnSlvrX = createSlvr(fluidEqn, elc, elcX, limiter, {0})
ionEqnSlvrX = createSlvr(fluidEqn, ion, ionX, limiter, {0})
emfEqnSlvrX = createSlvr(emfEqn, emf, emfX, limiter, {0})
elcEqnSlvrY = createSlvr(fluidEqn, elcX, elcY, limiter, {1})
ionEqnSlvrY = createSlvr(fluidEqn, ionX, ionY, limiter, {1})
emfEqnSlvrY = createSlvr(emfEqn, emfX, emfY, limiter, {1})
elcEqnSlvrZ = createSlvr(fluidEqn, elcY, elcZ, limiter, {2})
ionEqnSlvrZ = createSlvr(fluidEqn, ionY, ionZ, limiter, {2})
emfEqnSlvrZ = createSlvr(emfEqn, emfY, emfZ, limiter, {2})
slvrs = {
{elcEqnSlvrX, ionEqnSlvrX, emfEqnSlvrX},
{elcEqnSlvrY, ionEqnSlvrY, emfEqnSlvrY},
{elcEqnSlvrZ, ionEqnSlvrZ, emfEqnSlvrZ},
}
elcEqnLaxSlvrX = createSlvr(fluidEqnLax, elc, elcX, "zero", {0})
ionEqnLaxSlvrX = createSlvr(fluidEqnLax, ion, ionX, "zero", {0})
emfEqnLaxSlvrX = createSlvr(emfEqn, emf, emfX, "zero", {0})
elcEqnLaxSlvrY = createSlvr(fluidEqnLax, elcX, elcY, "zero", {1})
ionEqnLaxSlvrY = createSlvr(fluidEqnLax, ionX, ionY, "zero", {1})
emfEqnLaxSlvrY = createSlvr(emfEqn, emfX, emfY, "zero", {1})
elcEqnLaxSlvrZ = createSlvr(fluidEqnLax, elcY, elcZ, "zero", {2})
ionEqnLaxSlvrZ = createSlvr(fluidEqnLax, ionY, ionZ, "zero", {2})
emfEqnLaxSlvrZ = createSlvr(emfEqn, emfY, emfZ, "zero", {2})
laxSlvrs = {
{elcEqnLaxSlvrX, ionEqnLaxSlvrX, emfEqnLaxSlvrX},
{elcEqnLaxSlvrY, ionEqnLaxSlvrY, emfEqnLaxSlvrY},
{elcEqnLaxSlvrZ, ionEqnLaxSlvrZ, emfEqnLaxSlvrZ},
}
qInList = {q, qX, qY}
elcOutList = {elcX, elcY, elcZ}
ionOutList = {ionX, ionY, ionZ}
----------------------------------
-- HYPERBOLIC EQUATION UPDATERS --
----------------------------------
function updateHyperEqns(tCurr, tEnd, step)
local myStatus = true
local myDtSuggested = 1e3 * math.abs(tEnd-tCurr)
local useLaxFlux = false
for _, dir in ipairs({0, 1, 2}) do
applyBc(qInList[dir + 1], tCurr, tEnd, dir)
for _, slvr in ipairs(slvrs[dir + 1]) do
slvr:setCurrTime(tCurr)
local status, dtSuggested = slvr:advance(tEnd)
myStatus = status and myStatus
myDtSuggested = math.min(myDtSuggested, dtSuggested)
end
if (myStatus) then
if ((fluidEqn:checkInvariantDomain(elcOutList[dir + 1]) == false)
or (fluidEqn:checkInvariantDomain(ionOutList[dir + 1]) == false)) then
log("** Negative density/pressure after updateHyperEqns along %s", tostring(dir))
useLaxFlux = true
end
end
if (not myStatus) or useLaxFlux then
return myStatus, myDtSuggested, useLaxFlux
end
end
return myStatus, myDtSuggested, useLaxFlux
end
function updateHyperEqnsLax(tCurr, tEnd, step)
local myStatus = true
local myDtSuggested = 1e3*math.abs(tEnd-tCurr)
for _, dir in ipairs({0, 1, 2}) do
applyBc(qInList[dir + 1], tCurr, tEnd, dir)
for _, slvr in ipairs(laxSlvrs[dir + 1]) do
slvr:setCurrTime(tCurr)
local status, dtSuggested = slvr:advance(tEnd)
myStatus = status and myStatus
myDtSuggested = math.min(myDtSuggested, dtSuggested)
end
if not myStatus then
return myStatus, myDtSuggested
end
end
return myStatus, myDtSuggested
end
---------------------
-- SOURCE UPDATERS --
---------------------
srcSlvr = Updater.ImplicitFiveMomentSrc3D {
onGrid = grid,
numFluids = numFluids,
charge = charge,
mass = mass,
epsilon0 = epsilon0,
linearSolver = linearSolver,
elcErrorSpeedFactor = elcErrorSpeedFactor,
mgnErrorSpeedFactor = mgnErrorSpeedFactor,
elcErrorDampFactor = elcErrorDampFactor,
mgnErrorDampFactor = mgnErrorDampFactor,
}
function updateSource(elcIn, ionIn, emfIn, qIn, tCurr, tEnd)
srcSlvr:setOut( {elcIn, ionIn, emfIn} )
srcSlvr:setCurrTime(tCurr)
srcSlvr:advance(tEnd)
end
----------------------------------------
-- HYPERBOLIC-EQUATION SYSTEM SOLVERS --
----------------------------------------
function updateSystem(tCurr, tEnd, step)
local dthalf = 0.5 * (tEnd-tCurr)
updateSource(elc, ion, emf, q, tCurr, tCurr + dthalf)
local status, dtSuggested, useLaxFlux = updateHyperEqns(tCurr, tEnd, step)
updateSource(elcNew, ionNew, emfNew, qNew, tCurr, tCurr + dthalf)
if (status and not useLaxFlux) then
if ((fluidEqn:checkInvariantDomain(elcNew) == false)
or (fluidEqn:checkInvariantDomain(ionNew) == false)
or (qNew:hasNan())) then
log("** Negative density/pressure or Nan at end of updateSystem")
useLaxFlux = true
end
end
return status, dtSuggested, useLaxFlux
end
function updateSystemLax(tCurr, tEnd, step)
local dthalf = 0.5 * (tEnd - tCurr)
updateSource(elc, ion, emf, q, tCurr, tCurr + dthalf)
local status, dtSuggested = updateHyperEqnsLax(tCurr, tEnd, step)
updateSource(elcNew, ionNew, emfNew, qNew, tCurr, tCurr + dthalf)
return status, dtSuggested
end
----------
-- MAIN --
----------
function runSimulation(initDt)
local tCurr = 0.
local frame = 1
local tNextFrame = tCurr + tFrame
local stepTotal = 1
local step = 1
local stepInFrame = 1
local myDt = initDt
local dtSuggested = initDt
local status = true
local useLaxFlux = false
if (Lucee.IsRestarting) then
rFileName = "q_" .. Lucee.RestartFrame .. ".h5"
log("\nReading %s...", rFileName)
local time0 = os.time()
tCurr = q:read(rFileName)
local time = os.time() - time0
log("Reading %s...done. Took %gs", rFileName, time)
if not tCurr then
tCurr = tStart + Lucee.RestartFrame * tFrame
end
frame = Lucee.RestartFrame + 1
tNextFrame = tCurr + tFrame
log('\nRestarting from frame %d tCurr = %g\n', Lucee.RestartFrame, tCurr)
else
log("Initializing...")
q:set(init)
log("Initializing...done")
end
applyBc(q, 0, 0, 0)
applyBc(q, 0, 0, 1)
applyBc(q, 0, 0, 2)
if (not Lucee.IsRestarting) then
log(">>> Writing output 0 at t = %g...", tCurr)
writeData(q, 0, tCurr, "q")
end
log("Time-integration loop starting...")
while true do
if alwaysUseLaxFlux then
useLaxFlux = true
end
qDup:copy(q)
if (tCurr + myDt > tEnd) then
myDt = tEnd - tCurr
end
if useLaxFlux then
log("Taking step %5d (%4d in frame %3d; %5d in total); t = %10g, dt = %10g; using Lax fluxes",
step, stepInFrame, frame, stepTotal, tCurr, myDt)
status, dtSuggested = updateSystemLax(tCurr, tCurr+myDt, step)
if status then
useLaxFlux = false
end
else
log("Taking step %5d (%4d in frame %3d; %5d in total); t = %10g, dt = %10g",
step, stepInFrame, frame, stepTotal, tCurr, myDt)
status, dtSuggested, useLaxFlux = updateSystem(tCurr, tCurr+myDt, step)
end
stepTotal = stepTotal + 1
if not status then
log (" ** dt %g too large! Will retake step with dt %g; useLaxFlux = %s",
myDt, dtSuggested, tostring(useLaxFlux))
myDt = dtSuggested
q:copy(qDup)
elseif useLaxFlux then
log(" ** Negative pressure or density! Will retake step with Lax fluxes")
q:copy(qDup)
else
local badElc = fluidEqn:checkInvariantDomain(elcNew) == false
local badIon = fluidEqn:checkInvariantDomain(ionNew) == false
local badQ = qNew:hasNan()
if (badElc or badIon or badQ) then
log(" ** Error occurred!")
if badElc then log(" ** Negative electron pressure/density!") end
if badIon then log(" ** Negative ion pressure/density!") end
if badQ then log(" ** Nans occurred!") end
log(" ** Stopping simulation, dumping this and last states!")
writeData(qDup, -1, tCurr, "dead")
writeData(qNew, 0, tCurr+myDt, "dead")
break
end
q:copy(qNew)
tCurr = tCurr + myDt
if (tCurr > tNextFrame or tCurr >= tEnd) then
log(">>> Writing output %d at t = %g...", frame, tCurr)
qNew:applyCopyBc(2, "lower")
qNew:applyCopyBc(2, "upper")
writeData(qNew, frame, tCurr, "q")
frame = frame + 1
tNextFrame = tNextFrame + tFrame
stepInFrame = 0
end
myDt = dtSuggested
step = step + 1
stepInFrame = stepInFrame + 1
if (tCurr >= tEnd) then
break
end
end
end
end
------------------------
-- RUN THE SIMULATION --
------------------------
runSimulation(tFrame)
Raw
gkeyll_real_2d.lua
----------------------
-- HELPER FUNCTIONS --
----------------------
function log(...)
Lucee.logInfo(string.format(...))
end
function HERE()
local info = debug.getinfo(2)
local str = string.format("HERE: %d %s: %s",
info.currentline, info.source, tostring(info.name))
Lucee.logInfo(str)
end
function writeData(q, frame, tCurr, tag)
local time0 = os.time()
local fname = string.format("%s_%d.h5", tag, frame)
q:write(fname, tCurr)
log(" Writing %s took %gs", fname, os.time() - time0)
end
----------------
-- PARAMETERS --
----------------
-- physical constants
gasGamma = 5. / 3.
lightSpeed = 1.
mu0 = 1.
epsilon0 = 1. / mu0 / lightSpeed^2
-- kinetic parameters
numFluids = 2
mi = 1.
me = 1. / 25.
qi = 1.
qe = -1.
mass = {me, mi}
charge = {qe, qi}
massRatio = mi /me
-- problem specific parameters
wpe0_wce0 = 4.
n0 = 1.
betae0 = 0.5
betai0 = 0.5
wpe0 = math.sqrt(n0 * qe^2 / epsilon0 / me)
wce0 = wpe0 / wpe0_wce0
B0 = math.abs(wce0 * me / qe)
pmag0 = B0^2 / 2. / mu0
pe0 = pmag0 * betae0
pi0 = pmag0 * betai0
wpi0 = math.sqrt(n0 * qe^2 / epsilon0 / me)
wci0 = B0 * qi / mi
de0 = lightSpeed / wpe0
di0 = lightSpeed / wpi0
-- domain and grid parameters
xlo, xup, nx = -8. * di0, 8. * di0, 64
ylo, yup, ny = -16. * di0, 16. * di0, 128
lower = {xlo, ylo}
upper = {xup, yup}
cells = {nx, ny}
-- computational options
cfl = 0.9
cflm = 1.1 * cfl
limiter = "monotonized-centered"
elcErrorSpeedFactor = 0
mgnErrorSpeedFactor = 1
elcErrorDampFactor = 0
mgnErrorDampFactor = 0
linearSolver = "partialPivLu"
alwaysUseLaxFlux = false
-- I/O control
tEnd = 1 / wci0
tFrame = 0.1 / wci0
Lucee.IsRestarting = false
Lucee.RestartFrame = -1
outputFwd = false
log("%30s = %g", "massRatio", massRatio)
log("%30s = %g", "di0", di0)
log("%30s = %g", "1/wci0", 1./wci0)
log("%30s = %g", "tFrame", tFrame)
-----------------------
-- INITIAL CONDITION --
-----------------------
function init(x, y, z)
z = 0.
local rho = n0 * (1. + 0.1 * math.cos(x / di0))
local vx, vy, vz = 0.1, 0.05, 0.02 * math.sin(x / di0)
local rho_e = rho / (1 + massRatio)
local rho_i = rho - rho_e
local p_e = pe0
local p_i = pi0
local rhovx_e = rho_e * vx
local rhovy_e = rho_e * vy
local rhovz_e = rho_e * vz
local u_e = p_e / (gasGamma - 1.) + 0.5 * (vx^2 + vy^2 + vz^2) * rho_e
local rhovx_i = rho_i * vx
local rhovy_i = rho_i * vy
local rhovz_i = rho_i * vz
local u_i = p_i / (gasGamma - 1.) + 0.5 * (vx^2 + vy^2 + vz^2) * rho_i
local Bx, By, Bz = B0, B0 * 2, B0 * 3
local Ex = - vy * Bz + vz * By
local Ey = - vz * Bx + vx * Bz
local Ez = - vx * By + vy * Bx
return 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, 0., 0.
end
----------------------------
-- DECOMPOSITION AND GRID --
----------------------------
decomposition = DecompRegionCalc2D.CartGeneral {}
grid = Grid.RectCart2D {
lower = lower,
upper = upper,
cells = cells,
decomposition = decomposition,
periodicDirs ={0, 1},
}
----------
-- DATA --
----------
createData = function(numComponents)
if not numComponents then
numComponents = 18
end
return DataStruct.Field2D {
onGrid = grid,
numComponents = numComponents,
ghost = {2, 2},
--writeGhost = {2, 2},
}
end
q = createData()
qX = createData()
qNew = q
qY = qNew
qDup = createData()
getFields = function(q)
return q:alias(0, 5), q:alias(5, 10), q:alias(10, 18)
end
elc,ion,emf = getFields(q)
elcNew,ionNew,emfNew = getFields(qNew)
elcX,ionX,emfX = getFields(qX)
elcY,ionY,emfY = getFields(qY)
-------------------------
-- BOUNDARY CONDITIONS --
-------------------------
function applyBc(q, tCurr, tEnd, dir)
q:sync()
end
---------------------------------------
-- HYPERBOLIC EQUATIONS AND SOLVERS --
---------------------------------------
fluidEqn = HyperEquation.Euler {
gasGamma = gasGamma,
}
fluidEqnLax = HyperEquation.Euler {
gasGamma = gasGamma,
numericalFlux = "lax",
}
emfEqn = HyperEquation.PhMaxwell {
lightSpeed = lightSpeed,
elcErrorSpeedFactor = elcErrorSpeedFactor,
mgnErrorSpeedFactor = mgnErrorSpeedFactor,
}
createSlvr = function(equation, input, output, limiter, updateDirections)
local slvr = Updater.WavePropagation2D {
onGrid = grid,
equation = equation,
-- one of no-limiter, zero, min-mod, superbee,
-- van-leer, monotonized-centered, beam-warming
limiter = limiter,
cfl = cfl,
cflm = cflm,
updateDirections = updateDirections,
}
slvr:setIn( {input} )
slvr:setOut( {output} )
return slvr
end
elcEqnSlvrX = createSlvr(fluidEqn, elc, elcX, limiter, {0})
ionEqnSlvrX = createSlvr(fluidEqn, ion, ionX, limiter, {0})
emfEqnSlvrX = createSlvr(emfEqn, emf, emfX, limiter, {0})
elcEqnSlvrY = createSlvr(fluidEqn, elcX, elcY, limiter, {1})
ionEqnSlvrY = createSlvr(fluidEqn, ionX, ionY, limiter, {1})
emfEqnSlvrY = createSlvr(emfEqn, emfX, emfY, limiter, {1})
slvrs = {
{elcEqnSlvrX, ionEqnSlvrX, emfEqnSlvrX},
{elcEqnSlvrY, ionEqnSlvrY, emfEqnSlvrY},
}
elcEqnLaxSlvrX = createSlvr(fluidEqnLax, elc, elcX, "zero", {0})
ionEqnLaxSlvrX = createSlvr(fluidEqnLax, ion, ionX, "zero", {0})
emfEqnLaxSlvrX = createSlvr(emfEqn, emf, emfX, "zero", {0})
elcEqnLaxSlvrY = createSlvr(fluidEqnLax, elcX, elcY, "zero", {1})
ionEqnLaxSlvrY = createSlvr(fluidEqnLax, ionX, ionY, "zero", {1})
emfEqnLaxSlvrY = createSlvr(emfEqn, emfX, emfY, "zero", {1})
laxSlvrs = {
{elcEqnLaxSlvrX, ionEqnLaxSlvrX, emfEqnLaxSlvrX},
{elcEqnLaxSlvrY, ionEqnLaxSlvrY, emfEqnLaxSlvrY},
}
qInList = {q, qX}
elcOutList = {elcX, elcY}
ionOutList = {ionX, ionY}
----------------------------------
-- HYPERBOLIC EQUATION UPDATERS --
----------------------------------
function updateHyperEqns(tCurr, tEnd, step)
local myStatus = true
local myDtSuggested = 1e3 * math.abs(tEnd-tCurr)
local useLaxFlux = false
for _, dir in ipairs({0, 1}) do
applyBc(qInList[dir + 1], tCurr, tEnd, dir)
for _, slvr in ipairs(slvrs[dir + 1]) do
slvr:setCurrTime(tCurr)
local status, dtSuggested = slvr:advance(tEnd)
myStatus = status and myStatus
myDtSuggested = math.min(myDtSuggested, dtSuggested)
end
if (myStatus) then
if ((fluidEqn:checkInvariantDomain(elcOutList[dir + 1]) == false)
or (fluidEqn:checkInvariantDomain(ionOutList[dir + 1]) == false)) then
log("** Negative density/pressure after updateHyperEqns along %s", tostring(dir))
useLaxFlux = true
end
end
if (not myStatus) or useLaxFlux then
return myStatus, myDtSuggested, useLaxFlux
end
end
return myStatus, myDtSuggested, useLaxFlux
end
function updateHyperEqnsLax(tCurr, tEnd, step)
local myStatus = true
local myDtSuggested = 1e3*math.abs(tEnd-tCurr)
for _, dir in ipairs({0, 1}) do
applyBc(qInList[dir + 1], tCurr, tEnd, dir)
for _, slvr in ipairs(laxSlvrs[dir + 1]) do
slvr:setCurrTime(tCurr)
local status, dtSuggested = slvr:advance(tEnd)
myStatus = status and myStatus
myDtSuggested = math.min(myDtSuggested, dtSuggested)
end
if not myStatus then
return myStatus, myDtSuggested
end
end
return myStatus, myDtSuggested
end
---------------------
-- SOURCE UPDATERS --
---------------------
srcSlvr = Updater.ImplicitFiveMomentSrc2D {
onGrid = grid,
numFluids = numFluids,
charge = charge,
mass = mass,
epsilon0 = epsilon0,
linearSolver = linearSolver,
elcErrorSpeedFactor = elcErrorSpeedFactor,
mgnErrorSpeedFactor = mgnErrorSpeedFactor,
elcErrorDampFactor = elcErrorDampFactor,
mgnErrorDampFactor = mgnErrorDampFactor,
}
function updateSource(elcIn, ionIn, emfIn, qIn, tCurr, tEnd)
srcSlvr:setOut( {elcIn, ionIn, emfIn} )
srcSlvr:setCurrTime(tCurr)
srcSlvr:advance(tEnd)
end
----------------------------------------
-- HYPERBOLIC-EQUATION SYSTEM SOLVERS --
----------------------------------------
function updateSystem(tCurr, tEnd, step)
local dthalf = 0.5 * (tEnd-tCurr)
updateSource(elc, ion, emf, q, tCurr, tCurr + dthalf)
local status, dtSuggested, useLaxFlux = updateHyperEqns(tCurr, tEnd, step)
updateSource(elcNew, ionNew, emfNew, qNew, tCurr, tCurr + dthalf)
if (status and not useLaxFlux) then
if ((fluidEqn:checkInvariantDomain(elcNew) == false)
or (fluidEqn:checkInvariantDomain(ionNew) == false)
or (qNew:hasNan())) then
log("** Negative density/pressure or Nan at end of updateSystem")
useLaxFlux = true
end
end
return status, dtSuggested, useLaxFlux
end
function updateSystemLax(tCurr, tEnd, step)
local dthalf = 0.5 * (tEnd - tCurr)
updateSource(elc, ion, emf, q, tCurr, tCurr + dthalf)
local status, dtSuggested = updateHyperEqnsLax(tCurr, tEnd, step)
updateSource(elcNew, ionNew, emfNew, qNew, tCurr, tCurr + dthalf)
return status, dtSuggested
end
----------
-- MAIN --
----------
function runSimulation(initDt)
local tCurr = 0.
local frame = 1
local tNextFrame = tCurr + tFrame
local stepTotal = 1
local step = 1
local stepInFrame = 1
local myDt = initDt
local dtSuggested = initDt
local status = true
local useLaxFlux = false
if (Lucee.IsRestarting) then
rFileName = "q_" .. Lucee.RestartFrame .. ".h5"
log("\nReading %s...", rFileName)
local time0 = os.time()
tCurr = q:read(rFileName)
local time = os.time() - time0
log("Reading %s...done. Took %gs", rFileName, time)
if not tCurr then
tCurr = tStart + Lucee.RestartFrame * tFrame
end
frame = Lucee.RestartFrame + 1
tNextFrame = tCurr + tFrame
log('\nRestarting from frame %d tCurr = %g\n', Lucee.RestartFrame, tCurr)
else
log("Initializing...")
q:set(init)
log("Initializing...done")
end
applyBc(q, 0, 0, 0)
applyBc(q, 0, 0, 1)
applyBc(q, 0, 0, 2)
if (not Lucee.IsRestarting) then
log(">>> Writing output 0 at t = %g...", tCurr)
writeData(q, 0, tCurr, "q")
end
log("Time-integration loop starting...")
while true do
if alwaysUseLaxFlux then
useLaxFlux = true
end
qDup:copy(q)
if (tCurr + myDt > tEnd) then
myDt = tEnd - tCurr
end
if useLaxFlux then
log("Taking step %5d (%4d in frame %3d; %5d in total); t = %10g, dt = %10g; using Lax fluxes",
step, stepInFrame, frame, stepTotal, tCurr, myDt)
status, dtSuggested = updateSystemLax(tCurr, tCurr+myDt, step)
if status then
useLaxFlux = false
end
else
log("Taking step %5d (%4d in frame %3d; %5d in total); t = %10g, dt = %10g",
step, stepInFrame, frame, stepTotal, tCurr, myDt)
status, dtSuggested, useLaxFlux = updateSystem(tCurr, tCurr+myDt, step)
end
stepTotal = stepTotal + 1
if not status then
log (" ** dt %g too large! Will retake step with dt %g; useLaxFlux = %s",
myDt, dtSuggested, tostring(useLaxFlux))
myDt = dtSuggested
q:copy(qDup)
elseif useLaxFlux then
log(" ** Negative pressure or density! Will retake step with Lax fluxes")
q:copy(qDup)
else
local badElc = fluidEqn:checkInvariantDomain(elcNew) == false
local badIon = fluidEqn:checkInvariantDomain(ionNew) == false
local badQ = qNew:hasNan()
if (badElc or badIon or badQ) then
log(" ** Error occurred!")
if badElc then log(" ** Negative electron pressure/density!") end
if badIon then log(" ** Negative ion pressure/density!") end
if badQ then log(" ** Nans occurred!") end
log(" ** Stopping simulation, dumping this and last states!")
writeData(qDup, -1, tCurr, "dead")
writeData(qNew, 0, tCurr+myDt, "dead")
break
end
q:copy(qNew)
tCurr = tCurr + myDt
if (tCurr > tNextFrame or tCurr >= tEnd) then
log(">>> Writing output %d at t = %g...", frame, tCurr)
writeData(qNew, frame, tCurr, "q")
frame = frame + 1
tNextFrame = tNextFrame + tFrame
stepInFrame = 0
end
myDt = dtSuggested
step = step + 1
stepInFrame = stepInFrame + 1
if (tCurr >= tEnd) then
break
end
end
end
end
------------------------
-- RUN THE SIMULATION --
------------------------
runSimulation(tFrame)
@liangwang0734
Copy link
Author

liangwang0734 commented Oct 11, 2018 •
edited

  • applyCopyBc along the degenerate direction (periodic BC in gkeyll v1 does not work properly when there is only one cell)
  • on Trillian, the real 2d run took 44 s, the fake 2d run took 120 s.
    • path: /home/space/liang/scratch/gkeyll/2d_3d/

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment