GitHub

Time-Dependent Simulation

This tutorial walks you through your first AURORA simulation—from setting up the model to running a short time-dependent simulation with a flickering input flux.

Loading AURORA

using AURORA

Step 1: Create the atmospheric model

An AuroraModel bundles the numerical grids (altitude, energy, pitch angle), the background ionosphere (neutral and electron densities from MSIS/IRI), and collision cross-section data.

# Find (or generate) MSIS and IRI data files
# Default parameters: VISIONS-2 rocket launch conditions
# (2018-12-07 11:15 UTC, 76°N 5°E)
msis_file = find_msis_file()
iri_file  = find_iri_file()

model = AuroraModel(
    [100, 500],    # altitude limits [km]
    180:-30:0,     # pitch-angle bin edges [°] → 6 beams
    1000,          # maximum energy [eV]
    msis_file,
    iri_file,
    13             # magnetic field angle to zenith [°]
)
AuroraModel (not initialized):
├── AltitudeGrid(100.0 - 500.0 km, 401 points)
├── EnergyGrid(2.05 - 998.39 eV, 164 bins)
├── PitchAngleGrid(6 beams)
├── ScatteringData: (not initialized)
├── Ionosphere(401 altitudes)
├── Species: N2, O2, O
└── B angle to zenith: 13.0°

The constructor performs several steps automatically:

  • Builds an AltitudeGrid with fine spacing (~150 m) at the lowest altitudes and coarser spacing above.
  • Builds an EnergyGrid with adaptive bin widths (finer at low energies where numerous elastic and inelastic collisions need to be properly resolved).
  • Builds a PitchAngleGrid from the given angle edges.
  • Loads the neutrals (N₂, O₂, O) densities and temperature from the MSIS file.
  • Loads electron densities and temperature from the IRI file.
  • Precomputes electron collision cross sections and scattering data.

See Discretization to read more about the grids.

Scattering and cascading transfer matrices are cached automatically under the package's internal_data/ directory. If you want to rebuild them or avoid writing cache files during an exploratory run, pass a CachePolicy to AuroraModel or use the cache keywords on initialize!.

Step 2: Define the precipitating electron flux

The precipitating electrons are specified by an InputFlux, which combines an energy spectrum with a temporal modulation.

For this first example we use a simple 5 Hz flickering input.

flux = InputFlux(
  FlatSpectrum(1e-3; E_min=100),
  SinusoidalFlickering(5.0);
  beams=1,
  z_source=3000.0,
)
InputFlux:
├── Spectrum:   FlatSpectrum(IeE_tot=0.001 W/m², E_min=100.0 eV)
├── Modulation: SinusoidalFlickering(f=5.0 Hz, amplitude=1.0)
├── Beams:      [1]
└── Source altitude: 3000.0 km

Energy flux and modulation:

  • FlatSpectrum(1e-3; E_min=100) specifies a total energy flux of 1 mW/m², uniform in energy above 100 eV.
  • SinusoidalFlickering(5.0) applies a 5 Hz sinusoidal modulation with default amplitude 1.0, so the flux oscillates between 0 and 1 mW/m².

See the Input flux tutorial for other spectrum types, modulation types, and more details.

Beam and launch arguments:

  • The beams=1 argument means that only the most field-aligned downward-going beam carries the precipitating flux.
  • The z_source argument sets the altitude from which the electrons are launched. This introduces an energy-dependent travel-time delay, so higher-energy electrons reach the top of the ionosphere before lower-energy ones.

Step 3: Create and run the simulation

Bundle the model and input flux into an AuroraSimulation. For a time-dependent simulation, pass a TimeDependentMode specifying the total duration and output cadence:

savedir = mkpath(joinpath("data", "my_first_simulation"))

sim = AuroraSimulation(
    model,
    flux,
    savedir;
    mode=TimeDependentMode(
        duration=0.5,           # total simulation time [s]
        dt=0.01,                # output time step [s] (save every 10 ms)
        CFL_number=256,         # CFL stability parameter (higher → coarser internal stepping)
        max_memory_gb=4.0       # memory budget [GB] — controls loop partitioning
    )
)
AuroraSimulation (Time-dependent):
├── Model:       AuroraModel(AltitudeGrid(100.0 - 500.0 km, 401 points), EnergyGrid(2.05 - 998.39 eV, 164 bins), PitchAngleGrid(6 beams))
├── Flux:        InputFlux(FlatSpectrum(IeE_tot=0.001 W/m², E_min=100.0 eV), SinusoidalFlickering(f=5.0 Hz, amplitude=1.0), beams=[1])
├── Mode:        TimeDependentMode(duration=0.5, dt=0.01)
├── Savedir:     data/my_first_simulation
├── duration:      0.5 s
├── dt:            0.01 s
├── dt (internal): 0.002 s
├── CFL factor:    5
├── n_loop:        1
└── Cache:       not initialized

Understanding the time grid

The RefinedTimeGrid is constructed automatically from the TimeDependentMode parameters. It determines two key quantities:

  • dt_internal: the internal time step, chosen to satisfy the CFL condition. This is typically much smaller than the requested output dt.
  • n_loop: the number of time chunks the simulation is split into. Each chunk is solved independently and its results are saved before moving to the next. This keeps memory usage within the specified max_memory_gb budget.

You can inspect these:

sim.time
RefinedTimeGrid:
├── Duration:        0.5 s
├── dt:              0.01 s
├── Internal dt:     0.002 s
├── CFL factor:      5
├── n_loop:          1
├── n_save_per_loop: 50
└── steps / loop:    251

Running

run!(sim)
[ Info: Starting time-dependent simulation...

Solving [loop 1/1]   7%|██▍                              |  ETA: 0:00:13
Solving [loop 1/1]  17%|█████▋                           |  ETA: 0:00:13
Solving [loop 1/1]  26%|████████▋                        |  ETA: 0:00:11
Solving [loop 1/1]  35%|███████████▌                     |  ETA: 0:00:09
Solving [loop 1/1]  44%|██████████████▌                  |  ETA: 0:00:07
Solving [loop 1/1]  54%|█████████████████▊               |  ETA: 0:00:06
Solving [loop 1/1]  63%|████████████████████▉            |  ETA: 0:00:05
Solving [loop 1/1]  74%|████████████████████████▍        |  ETA: 0:00:03
Solving [loop 1/1]  95%|███████████████████████████████▍ |  ETA: 0:00:01
Solving [loop 1/1] 100%|█████████████████████████████████| Time: 0:00:10

Internally, AURORA refines the time grid as needed to satisfy the CFL criterion and then sweeps through energies from high to low, advancing the Crank-Nicolson scheme for each energy. The solver loops over time chunks; within each chunk it solves the full energy sweep and saves the results before moving to the next.

Step 4: Post-process

AURORA provides several analysis functions that compute derived quantities from the raw electron flux and save them alongside the simulation output:

make_Ie_top_file(sim)              # boundary condition (input flux applied at top)
make_volume_excitation_file(sim)   # volumetric excitation rates for optical emissions
make_column_excitation_file(sim)   # column-integrated excitation rates
make_current_file(sim)             # field-aligned electron currents and energy fluxes
make_heating_rate_file(sim)        # electron heating rates
make_psd_file(sim)                 # electron phase-space density f(E, θ) and F(v∥)
Top flux saved in data/my_first_simulation/analysis/Ie_top.nc
Volume excitation rates saved in data/my_first_simulation/analysis/volume_excitation.nc
Column excitation rates saved in data/my_first_simulation/analysis/column_excitation.nc
Currents saved in data/my_first_simulation/analysis/currents.nc
Heating rates saved in data/my_first_simulation/analysis/heating_rate.nc
Converting Ie to PSD.
PSD saved in data/my_first_simulation/analysis/psd.nc

The output files are saved in data/my_first_simulation/:

readdir(savedir)
4-element Vector{String}:
 "analysis"
 "config.toml"
 "inputs"
 "simulation_data.nc"

Step 5: Visualize

AURORA.jl provides helper plotting functions to visualize results through a Makie extension. Install and load a Makie backend to access them (more information about this in Visualization).

Input flux

plot_input can be used to inspect the prescribed input spectrum directly from the simulation object:

using CairoMakie

fig = plot_input(sim)
Example block output

Volume excitation rate

plot_excitation! can be used to visualize the volume excitation rate for a given emission line over altitude and time. 

using CairoMakie
vol = load_volume_excitation(savedir)

fig = Figure()
ax  = Axis(fig[1, 1];
           xlabel = "Time (s)",
           ylabel = "Altitude (km)",
           title = "4278 Å")
hm  = plot_excitation!(ax, vol; field = :Q4278)
Colorbar(fig[1, 2], hm; label = "photons/m³/s")
fig
Example block output

Altitude profiles for specific time-steps can also be visualized using the time_index keyword argument.

using CairoMakie
vol = load_volume_excitation(savedir)

fig = Figure()
ax  = Axis(fig[1, 1];
           xlabel = "Excitation rate (photons/m³/s)",
           ylabel = "Altitude (km)",
           xscale = log10,
           title = "4278 Å")
plot_excitation!(ax, vol; field = :Q4278, time_index = 15)
fig
Example block output

Column emission intensity

plot_column_excitation! can be used to visualize the column-integrated emission intensities (as would be seen from the ground) for several emission lines (by default 4278 Å, 6730 Å, 7774 Å, 8446 Å, O1D and O1S).

using CairoMakie
col = load_column_excitation(savedir)

fig = Figure()
ax  = Axis(fig[1, 1]; xlabel = "Time (s)", ylabel = "Intensity (R)", yscale = log10)
plot_column_excitation!(ax, col)
Legend(fig[1, 2], ax)
fig
Example block output

Animation

animate_Ie_in_time produces an animation of the electron flux distribution plotted over altitude and energy, with different panels for different pitch-angle streams, stepping over time:

using CairoMakie
animate_Ie_in_time(savedir; framerate = 15)
The animation will be saved at animation.mp4.
(1/1) [12:45:04] Load coordinates... scan colorrange... create figure... animate.

Saving animation... done!