Last active
October 5, 2018 01:00
-
-
Save liangwang0734/5f11c9e874ae868d2e596721348be82a to your computer and use it in GitHub Desktop.
test gkeyll
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import h5py | |
import numpy | |
numpy.set_printoptions(precision=16) | |
def printValues(fname): | |
fp = h5py.File(fname) | |
field = fp["StructGridField"][...] | |
print(field) | |
fp.close() | |
for fname in ["gem_qBeforeSrc.h5", "gem_qAfterSrc.h5", "gem_qSweep.h5"]: | |
print(fname) | |
printValues(fname) | |
print("\n") |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
-- GEM-challenge problem | |
Pi = Lucee.Pi | |
log = Lucee.logInfo | |
-- physical parameters | |
gasGamma = 5./3. | |
elcCharge = -1.0 | |
ionCharge = 1.0 | |
ionMass = 1.0 | |
elcMass = 0.25 | |
lightSpeed = 1.0 | |
epsilon0 = 1.0 | |
mgnErrorSpeedFactor = 1.0 | |
Lx = 25.6 | |
Ly = 12.8 | |
n0 = 1.0 | |
VAe = 0.5 | |
plasmaBeta = 1.0 | |
lambda = 0.5 | |
TiOverTe = 5.0 | |
nbOverN0 = 0.2 | |
pert = 0.1 | |
Valf = VAe*math.sqrt(elcMass/ionMass) | |
B0 = Valf*math.sqrt(n0*ionMass) | |
OmegaCi0 = ionCharge*B0/ionMass | |
psi0 = pert*B0 | |
-- resolution and time-stepping | |
NX = 64 | |
NY = 32 | |
NX = 2 | |
NY = 2 | |
cfl = 1. | |
tStart = 0.0 | |
tEnd = 1 | |
nFrames = 1 | |
log(string.format("elcMass/ionMass=1/%g",ionMass/elcMass)) | |
log(string.format("Lx=%gdi=%gde", Lx,Lx*math.sqrt(ionMass/elcMass))) | |
log(string.format("plasmaBeta=%g", plasmaBeta)) | |
log(string.format("Valf/c=%g", Valf)) | |
log(string.format("lambda/di=%g", lambda)) | |
log(string.format("TiOverTe=%g", TiOverTe)) | |
log(string.format("nbOverN0=%g", nbOverN0)) | |
log(string.format("pert=%g", pert)) | |
log(string.format("tEnd=%g, nFrames=%d",tEnd,nFrames)) | |
log(string.format("NX=%d,dx=%gdi=%gde,cfl=%g\n", NX,Lx/NX,Lx/NX*math.sqrt(ionMass/elcMass),cfl)) | |
------------------------------------------------ | |
-- COMPUTATIONAL DOMAIN, DATA STRUCTURE, ETC. -- | |
------------------------------------------------ | |
-- decomposition object | |
decomp = DecompRegionCalc2D.CartGeneral {} | |
-- computational domain | |
grid = Grid.RectCart2D { | |
lower = {0, 0}, | |
upper = {1., 1.}, | |
cells = {4, 1,}, | |
decomposition = decomp, | |
periodicDirs = {0,1}, | |
} | |
-- solution | |
q = DataStruct.Field2D { | |
onGrid = grid, | |
numComponents = 18, | |
ghost = {2, 2}, | |
} | |
-- solution after update along X (ds algorithm) | |
qX = DataStruct.Field2D { | |
onGrid = grid, | |
numComponents = 18, | |
ghost = {2, 2}, | |
} | |
-- final updated solution | |
qNew = DataStruct.Field2D { | |
onGrid = grid, | |
numComponents = 18, | |
ghost = {2, 2}, | |
} | |
-- duplicate copy in case we need to take the step again | |
qDup = DataStruct.Field2D { | |
onGrid = grid, | |
numComponents = 18, | |
ghost = {2, 2}, | |
} | |
qNewDup = DataStruct.Field2D { | |
onGrid = grid, | |
numComponents = 18, | |
ghost = {2, 2}, | |
} | |
-- aliases to various sub-systems | |
elcFluid = q:alias(0, 5) | |
ionFluid = q:alias(5, 10) | |
emField = q:alias(10, 18) | |
elcFluidX = qX:alias(0, 5) | |
ionFluidX = qX:alias(5, 10) | |
emFieldX = qX:alias(10, 18) | |
elcFluidNew = qNew:alias(0, 5) | |
ionFluidNew = qNew:alias(5, 10) | |
emFieldNew = qNew:alias(10, 18) | |
----------------------- | |
-- INITIAL CONDITION -- | |
----------------------- | |
-- initial conditions | |
function init(x,y,z) | |
local n = 1 + .3 * x | |
print(x, n) | |
local rhoe = 1 + .3 * x | |
local rhovxe = 2 + .3 * x | |
local rhovye = 3 + .3 * x | |
local rhovze = 4 + .3 * x | |
local pe = 5 + .3 * x | |
local ue = pe / (gasGamma-1) + 0.5 * (rhovxe^2 + rhovye^2 + rhovze^2) / rhoe | |
local n = 1.4 | |
local rhoi = 1.4 | |
local rhovxi = 2.4 | |
local rhovyi = 3.4 | |
local rhovzi = 4.4 | |
local pi = 5.4 | |
local ui = pi / (gasGamma-1) + 0.5 * (rhovxi^2 + rhovyi^2 + rhovzi^2) / rhoi | |
local Ex = 1 + .1 * x | |
local Ey = 1 + .2 * x | |
local Ez = 1 + .3 * x | |
local Bx = 2 + .1 * x | |
local By = 2 + .2 * x | |
local Bz = 2 + .3 * x | |
return rhoe, rhovxe, rhovye, rhovze, ue, | |
rhoi, rhovxi, rhovyi, rhovzi, ui, | |
Ex, Ey, Ez, Bx, By, Bz, 0., 0. | |
end | |
------------------------ | |
-- Boundary Condition -- | |
------------------------ | |
-- boundary applicator objects for fluids and fields | |
bcElcCopy = BoundaryCondition.Copy { components = {0, 4} } | |
bcElcWall = BoundaryCondition.ZeroNormal { components = {1, 2, 3} } | |
bcIonCopy = BoundaryCondition.Copy { components = {5, 9} } | |
bcIonWall = BoundaryCondition.ZeroNormal { components = {6, 7, 8} } | |
bcElcFld = BoundaryCondition.ZeroTangent { components = {10, 11, 12} } | |
bcMgnFld = BoundaryCondition.ZeroNormal { components = {13, 14, 15} } | |
bcPot = BoundaryCondition.Copy { components = {16, 17} } | |
--FIXME: fact in bcPot | |
-- create boundary condition object | |
function createBc(myDir, myEdge) | |
local bc = Updater.Bc2D { | |
onGrid = grid, | |
-- boundary conditions to apply | |
boundaryConditions = { | |
bcElcCopy, bcElcWall, | |
bcIonCopy, bcIonWall, | |
bcElcFld, bcMgnFld, bcPot, | |
}, | |
-- direction to apply | |
dir = myDir, | |
-- edge to apply on | |
edge = myEdge, | |
} | |
return bc | |
end | |
-- create updaters to apply boundary conditions | |
bcBottom = createBc(1, "lower") | |
bcTop = createBc(1, "upper") | |
-- function to apply boundary conditions to specified field | |
function applyBc(fld, tCurr, myDt) | |
-- for i,bc in ipairs({bcBottom, bcTop}) do | |
-- bc:setOut( {fld} ) | |
-- bc:advance(tCurr+myDt) | |
-- end | |
-- sync ghost cells | |
fld:sync() | |
end | |
---------------------- | |
-- EQUATION SOLVERS -- | |
---------------------- | |
-- regular Euler equations | |
elcEulerEqn = HyperEquation.Euler { | |
gasGamma = gasGamma, | |
} | |
ionEulerEqn = HyperEquation.Euler { | |
gasGamma = gasGamma, | |
} | |
-- (Lax equations are used to fix negative pressure/density) | |
elcEulerLaxEqn = HyperEquation.Euler { | |
gasGamma = gasGamma, | |
numericalFlux = "lax", | |
} | |
ionEulerLaxEqn = HyperEquation.Euler { | |
gasGamma = gasGamma, | |
numericalFlux = "lax", | |
} | |
maxwellEqn = HyperEquation.PhMaxwell { | |
lightSpeed = lightSpeed, | |
elcErrorSpeedFactor = 0.0, | |
mgnErrorSpeedFactor = mgnErrorSpeedFactor | |
} | |
-- ds solvers for regular Euler equations along X | |
elcFluidSlvrDir0 = Updater.WavePropagation2D { | |
onGrid = grid, | |
equation = elcEulerEqn, | |
-- one of no-limiter, min-mod, superbee, | |
-- van-leer, monotonized-centered, beam-warming | |
limiter = "van-leer", | |
cfl = cfl, | |
cflm = 1.1*cfl, | |
updateDirections = {0} -- directions to update | |
} | |
ionFluidSlvrDir0 = Updater.WavePropagation2D { | |
onGrid = grid, | |
equation = ionEulerEqn, | |
limiter = "van-leer", | |
cfl = cfl, | |
cflm = 1.1*cfl, | |
updateDirections = {0} | |
} | |
maxSlvrDir0 = Updater.WavePropagation2D { | |
onGrid = grid, | |
equation = maxwellEqn, | |
limiter = "van-leer", | |
cfl = cfl, | |
cflm = 1.1*cfl, | |
updateDirections = {0} | |
} | |
-- ds solvers for regular Euler equations along Y | |
elcFluidSlvrDir1 = Updater.WavePropagation2D { | |
onGrid = grid, | |
equation = elcEulerEqn, | |
limiter = "van-leer", | |
cfl = cfl, | |
cflm = 1.1*cfl, | |
updateDirections = {1} | |
} | |
ionFluidSlvrDir1 = Updater.WavePropagation2D { | |
onGrid = grid, | |
equation = ionEulerEqn, | |
limiter = "van-leer", | |
cfl = cfl, | |
cflm = 1.1*cfl, | |
updateDirections = {1} | |
} | |
maxSlvrDir1 = Updater.WavePropagation2D { | |
onGrid = grid, | |
equation = maxwellEqn, | |
limiter = "van-leer", | |
cfl = cfl, | |
cflm = 1.1*cfl, | |
updateDirections = {1} | |
} | |
-- ds solvers for Lax Euler equations along X | |
elcLaxSlvrDir0 = Updater.WavePropagation2D { | |
onGrid = grid, | |
equation = elcEulerLaxEqn, | |
limiter = "van-leer", | |
cfl = cfl, | |
cflm = 1.1*cfl, | |
updateDirections = {0} | |
} | |
ionLaxSlvrDir0 = Updater.WavePropagation2D { | |
onGrid = grid, | |
equation = ionEulerLaxEqn, | |
limiter = "van-leer", | |
cfl = cfl, | |
cflm = 1.1*cfl, | |
updateDirections = {0} | |
} | |
maxLaxSlvrDir0 = Updater.WavePropagation2D { | |
onGrid = grid, | |
equation = maxwellEqn, | |
limiter = "van-leer", | |
cfl = cfl, | |
cflm = 1.1*cfl, | |
updateDirections = {0} | |
} | |
-- ds solvers for Lax Euler equations along Y | |
elcLaxSlvrDir1 = Updater.WavePropagation2D { | |
onGrid = grid, | |
equation = elcEulerLaxEqn, | |
limiter = "van-leer", | |
cfl = cfl, | |
cflm = 1.1*cfl, | |
updateDirections = {1} | |
} | |
ionLaxSlvrDir1 = Updater.WavePropagation2D { | |
onGrid = grid, | |
equation = ionEulerLaxEqn, | |
limiter = "van-leer", | |
cfl = cfl, | |
cflm = 1.1*cfl, | |
updateDirections = {1} | |
} | |
maxLaxSlvrDir1 = Updater.WavePropagation2D { | |
onGrid = grid, | |
equation = maxwellEqn, | |
limiter = "van-leer", | |
cfl = cfl, | |
cflm = 1.1*cfl, | |
updateDirections = {1} | |
} | |
-- updater for source terms | |
sourceSlvr = Updater.ImplicitFiveMomentSrc2D { | |
onGrid = grid, | |
numFluids = 2, | |
charge = {elcCharge, ionCharge}, | |
mass = {elcMass, ionMass}, | |
epsilon0 = epsilon0, | |
-- linear solver to use: one of partialPivLu or colPivHouseholderQr | |
linearSolver = "partialPivLu", | |
hasStaticField = false, | |
} | |
-- function to update source terms | |
function updateSource(elcIn, ionIn, emIn, tCurr, t) | |
print("updateSource dt", t-tCurr) | |
sourceSlvr:setOut( {elcIn, ionIn, emIn} ) | |
sourceSlvr:setCurrTime(tCurr) | |
sourceSlvr:advance(t) | |
-- momentum drag relaxation | |
end | |
-- function to update the fluid and field using dimensional splitting | |
function updateFluidsAndField(tCurr, t) | |
local myStatus = true | |
local myDtSuggested = 1e3*math.abs(t-tCurr) | |
local useLaxSolver = False | |
-- X-direction updates | |
applyBc(q, tCurr, t-tCurr) | |
for i,slvr in ipairs({elcFluidSlvrDir0, ionFluidSlvrDir0, maxSlvrDir0}) do | |
slvr:setCurrTime(tCurr) | |
local status, dtSuggested = slvr:advance(t) | |
myStatus = status and myStatus | |
myDtSuggested = math.min(myDtSuggested, dtSuggested) | |
end | |
if ((elcEulerEqn:checkInvariantDomain(elcFluidX) == false) | |
or (ionEulerEqn:checkInvariantDomain(ionFluidX) == false)) then | |
useLaxSolver = true | |
end | |
if ((myStatus == false) or (useLaxSolver == true)) then | |
return myStatus, myDtSuggested, useLaxSolver | |
end | |
-- apply BCs to intermediate update after X sweep | |
applyBc(qX, tCurr, t-tCurr) | |
-- Y-direction updates | |
for i,slvr in ipairs({elcFluidSlvrDir1, ionFluidSlvrDir1, maxSlvrDir1}) do | |
slvr:setCurrTime(tCurr) | |
local status, dtSuggested = slvr:advance(t) | |
myStatus = status and myStatus | |
myDtSuggested = math.min(myDtSuggested, dtSuggested) | |
end | |
if ((elcEulerEqn:checkInvariantDomain(elcFluidNew) == false) | |
or (ionEulerEqn:checkInvariantDomain(ionFluidNew) == false)) then | |
useLaxSolver = true | |
end | |
return myStatus, myDtSuggested, useLaxSolver | |
end | |
write = true | |
-- function to take one time-step with Euler solver | |
function solveTwoFluidSystem(tCurr, t) | |
local dthalf = 0.5*(t-tCurr) | |
-- update source terms | |
if write then | |
q:write("qBeforeSrc.h5", 0) | |
end | |
updateSource(elcFluid, ionFluid, emField, tCurr, tCurr+dthalf) | |
applyBc(q, tCurr, t-tCurr) | |
if write then | |
q:write("qAfterSrc.h5", 0) | |
end | |
-- update fluids and fields | |
local status, dtSuggested, useLaxSolver = updateFluidsAndField(tCurr, t) | |
if write then | |
q:write("qSweep.h5", 0) | |
end | |
-- update source terms | |
updateSource(elcFluidNew, ionFluidNew, emFieldNew, tCurr, tCurr+dthalf) | |
applyBc(qNew, tCurr, t-tCurr) | |
if write then | |
write = false | |
end | |
os.exit() | |
return status, dtSuggested,useLaxSolver | |
end | |
-- function to update the fluid and field using dimensional splitting Lax scheme | |
function updateFluidsAndFieldLax(tCurr, t) | |
local myStatus = true | |
local myDtSuggested = 1e3*math.abs(t-tCurr) | |
applyBc(q, tCurr, t-tCurr) | |
for i,slvr in ipairs({elcLaxSlvrDir0, ionLaxSlvrDir0, maxLaxSlvrDir0}) do | |
slvr:setCurrTime(tCurr) | |
local status, dtSuggested = slvr:advance(t) | |
myStatus = status and myStatus | |
myDtSuggested = math.min(myDtSuggested, dtSuggested) | |
end | |
applyBc(qX, tCurr, t-tCurr) | |
-- Y-direction updates | |
for i,slvr in ipairs({elcLaxSlvrDir1, ionLaxSlvrDir1, maxLaxSlvrDir1}) do | |
slvr:setCurrTime(tCurr) | |
local status, dtSuggested = slvr:advance(t) | |
myStatus = status and myStatus | |
myDtSuggested = math.min(myDtSuggested, dtSuggested) | |
end | |
return myStatus, myDtSuggested | |
end | |
-- function to take one time-step with Lax Euler solver | |
function solveTwoFluidLaxSystem(tCurr, t) | |
local dthalf = 0.5*(t-tCurr) | |
-- update source terms | |
updateSource(elcFluid, ionFluid, emField, tCurr, tCurr+dthalf) | |
applyBc(q, tCurr, t-tCurr) | |
-- update fluids and fields | |
local status, dtSuggested = updateFluidsAndFieldLax(tCurr, t) | |
-- update source terms | |
updateSource(elcFluidNew, ionFluidNew, emFieldNew, tCurr, tCurr+dthalf) | |
applyBc(qNew, tCurr, t-tCurr) | |
return status, dtSuggested | |
end | |
---------------------------- | |
-- DIAGNOSIS AND DATA I/O -- | |
---------------------------- | |
-- dynvector to store integrated flux | |
byAlias = qNew:alias(14, 15) | |
byFlux = DataStruct.DynVector { numComponents = 1 } | |
byFluxCalc = Updater.IntegrateFieldAlongLine2D { | |
onGrid = grid, | |
-- start cell | |
startCell = {0, NY/2}, | |
-- direction to integrate in | |
dir = 0, | |
-- number of cells to integrate | |
numCells = NX, | |
-- integrand | |
integrand = function (by) | |
return math.abs(by) | |
end, | |
} | |
byFluxCalc:setIn( {byAlias} ) | |
byFluxCalc:setOut( {byFlux} ) | |
-- dynvector to store Ez at X-point | |
ezAlias = qNew:alias(12, 13) | |
xpointEz = DataStruct.DynVector { numComponents = 1 } | |
xpointEzRec = Updater.RecordFieldInCell2D { | |
onGrid = grid, | |
-- index of cell to record | |
cellIndex = {(NX-1)/2, (NY-1)/2}, | |
} | |
xpointEzRec:setIn( {ezAlias} ) | |
xpointEzRec:setOut( {xpointEz} ) | |
-- dynvector to store number density at X-point | |
neAlias = qNew:alias(0, 1) | |
xpointNe = DataStruct.DynVector { numComponents = 1 } | |
xpointNeRec = Updater.RecordFieldInCell2D { | |
onGrid = grid, | |
-- index of cell to record | |
cellIndex = {(NX-1)/2, (NY-1)/2}, | |
} | |
xpointNeRec:setIn( {neAlias} ) | |
xpointNeRec:setOut( {xpointNe} ) | |
-- dynvector to store electron uz at X-point | |
uzeAlias = qNew:alias(3, 4) | |
xpointUze = DataStruct.DynVector { numComponents = 1 } | |
xpointUzeRec = Updater.RecordFieldInCell2D { | |
onGrid = grid, | |
-- index of cell to record | |
cellIndex = {(NX-1)/2, (NY-1)/2}, | |
} | |
xpointUzeRec:setIn( {uzeAlias} ) | |
xpointUzeRec:setOut( {xpointUze} ) | |
-- dynvector to store ion uz at X-point | |
uziAlias = qNew:alias(8, 9) | |
xpointUzi = DataStruct.DynVector { numComponents = 1 } | |
xpointUziRec = Updater.RecordFieldInCell2D { | |
onGrid = grid, | |
-- index of cell to record | |
cellIndex = {(NX-1)/2, (NY-1)/2}, | |
} | |
xpointUziRec:setIn( {uziAlias} ) | |
xpointUziRec:setOut( {xpointUzi} ) | |
-- compute diagnostic | |
function calcDiagnostics(tCurr, dt) | |
for i,diag in ipairs({byFluxCalc, xpointEzRec, xpointNeRec, xpointUzeRec, xpointUziRec}) do | |
diag:setCurrTime(tCurr) | |
diag:advance(tCurr+dt) | |
end | |
end | |
-- write data to H5 files | |
function writeFields(frame, t) | |
qNew:write( string.format("q_%d.h5", frame), t ) | |
byFlux:write( string.format("byFlux_%d.h5", frame) ) | |
xpointEz:write(string.format("xpointEz_%d.h5", frame) ) | |
xpointNe:write(string.format("xpointNe_%d.h5", frame) ) | |
xpointUze:write(string.format("xpointUze_%d.h5", frame) ) | |
xpointUzi:write(string.format("xpointUzi_%d.h5", frame) ) | |
end | |
Fluids = qNew:alias(0, 10) | |
function writeFieldsBwd(frame, t) | |
Fluids:write( string.format("Fluids_%d_bwd.h5", frame), t ) | |
end | |
function writeFieldsFwd(frame, t) | |
Fluids:write( string.format("Fluids_%d_fwd.h5", frame), t ) | |
end | |
---------------------------- | |
-- TIME-STEPPING FUNCTION -- | |
---------------------------- | |
function runSimulation(tStart, tEnd, nFrames, initDt) | |
local frame = 1 | |
local tFrame = (tEnd-tStart)/nFrames | |
local nextIOt = tFrame | |
local step = 1 | |
local tCurr = tStart | |
local myDt = initDt | |
local status, dtSuggested | |
local useLaxSolver = false | |
local doWriteFieldsFwd = false | |
-- the grand loop | |
while true do | |
-- copy q and qNew in case we need to take this step again | |
qDup:copy(q) | |
qNewDup:copy(qNew) | |
-- if needed adjust dt to hit tEnd exactly | |
if (tCurr+myDt > tEnd) then | |
myDt = tEnd-tCurr | |
end | |
-- advance fluids and fields | |
if (useLaxSolver) then | |
-- call Lax solver if positivity violated | |
log (string.format(" Taking step %5d at time %6g with dt %g (using Lax solvers)", step, tCurr, myDt)) | |
status, dtSuggested = solveTwoFluidLaxSystem(tCurr, tCurr+myDt) | |
useLaxSolver = false | |
else | |
log (string.format(" Taking step %5d at time %6g with dt %g", step, tCurr, myDt)) | |
status, dtSuggested, useLaxSolver = solveTwoFluidSystem(tCurr, tCurr+myDt) | |
end | |
if (status == false) then | |
-- time-step too large | |
log (string.format(" ** Time step %g too large! Will retake with dt %g", myDt, dtSuggested)) | |
myDt = dtSuggested | |
qNew:copy(qNewDup) | |
q:copy(qDup) | |
elseif (useLaxSolver == true) then | |
-- negative density/pressure occured | |
log (string.format(" ** Negative pressure or density at %8g! Will retake step with Lax fluxes", tCurr+myDt)) | |
q:copy(qDup) | |
qNew:copy(qNewDup) | |
else | |
-- check if a nan occured | |
if (qNew:hasNan()) then | |
log (string.format(" ** NaN occured at %g! Stopping simulation", tCurr)) | |
break | |
end | |
-- compute diagnostics | |
calcDiagnostics(tCurr, myDt) | |
-- copy updated solution back | |
q:copy(qNew) | |
if (tCurr+myDt <= nextIOt and tCurr+myDt+dtSuggested > nextIOt) then | |
log (string.format(" Writing electron data (bwd) at time %g (frame %d) ...", tCurr+myDt, frame)) | |
writeFieldsBwd(frame, tCurr+myDt) | |
end | |
if doWriteFieldsFwd then | |
log (string.format(" Writing electron data (fwd) at time %g (frame %d) ...", tCurr+myDt, frame-1)) | |
writeFieldsFwd(frame-1, tCurr+myDt) | |
doWriteFieldsFwd = false | |
end | |
-- write out data | |
if (tCurr+myDt > nextIOt or tCurr+myDt >= tEnd) then | |
log (string.format(" Writing data at time %g (frame %d) ...\n", tCurr+myDt, frame)) | |
writeFields(frame, tCurr+myDt) | |
frame = frame + 1 | |
nextIOt = nextIOt + tFrame | |
step = 0 | |
doWriteFieldsFwd = true | |
end | |
tCurr = tCurr + myDt | |
myDt = dtSuggested | |
step = step + 1 | |
-- check if done | |
if (tCurr >= tEnd) then | |
break | |
end | |
end | |
end -- end of time-step loop | |
return dtSuggested | |
end | |
---------------------------- | |
-- RUNNING THE SIMULATION -- | |
---------------------------- | |
-- setup initial condition | |
q:set(init) | |
q:sync() | |
qNew:copy(q) | |
-- set input/output arrays for various solvers | |
elcFluidSlvrDir0:setIn( {elcFluid} ) | |
elcFluidSlvrDir0:setOut( {elcFluidX} ) | |
ionFluidSlvrDir0:setIn( {ionFluid} ) | |
ionFluidSlvrDir0:setOut( {ionFluidX} ) | |
maxSlvrDir0:setIn( {emField} ) | |
maxSlvrDir0:setOut( {emFieldX} ) | |
elcFluidSlvrDir1:setIn( {elcFluidX} ) | |
elcFluidSlvrDir1:setOut( {elcFluidNew} ) | |
ionFluidSlvrDir1:setIn( {ionFluidX} ) | |
ionFluidSlvrDir1:setOut( {ionFluidNew} ) | |
maxSlvrDir1:setIn( {emFieldX} ) | |
maxSlvrDir1:setOut( {emFieldNew} ) | |
elcLaxSlvrDir0:setIn( {elcFluid} ) | |
elcLaxSlvrDir0:setOut( {elcFluidX} ) | |
ionLaxSlvrDir0:setIn( {ionFluid} ) | |
ionLaxSlvrDir0:setOut( {ionFluidX} ) | |
maxLaxSlvrDir0:setIn( {emField} ) | |
maxLaxSlvrDir0:setOut( {emFieldX} ) | |
elcLaxSlvrDir1:setIn( {elcFluidX} ) | |
elcLaxSlvrDir1:setOut( {elcFluidNew} ) | |
ionLaxSlvrDir1:setIn( {ionFluidX} ) | |
ionLaxSlvrDir1:setOut( {ionFluidNew} ) | |
maxLaxSlvrDir1:setIn( {emFieldX} ) | |
maxLaxSlvrDir1:setOut( {emFieldNew} ) | |
-- apply BCs on initial conditions | |
-- applyBc(q, 0.0, 0.0) | |
-- applyBc(qNew, 0.0, 0.0) | |
-- write initial conditions | |
calcDiagnostics(0.0, 0.0) | |
writeFields(0, 0.0) | |
initDt = 100.0 | |
runSimulation(tStart, tEnd, nFrames, initDt) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment