Simulation

Main Simulation Functions

simulate_waveDM

The high-level entry point. It builds the initial conditions (using AstroIC.jl), creates a SimulationGrid, calls SPE3D_waveDM to evolve the SPE, then handles visualization and post-processing. Most of the keyword arguments listed below can be supplied directly to simulate_waveDM; the low-level-only ones are marked in the tables.

SPE3D_waveDM

The lower-level entry point. It expects you to provide the initial wavefunction IC and (optionally) the baryonic density/potential/forces yourself. Use simulate_waveDM unless you are integrating WaveDM.jl into a larger pipeline.

Simulation Workflow

flowchart LR
    A[Grid setup] --> B[Initial conditions]
    B --> C[Baryon setup<br/>mesh / static / dynamic]
    C --> D[IC potential]
    D --> E[Visualization init]
    E --> F[KDK main loop]
    F -->|every StepsBetweenSnapshots| E
    F --> G[Outputs<br/>CSV, JLD2, PNG, MP4]

1. Grid setup — uniform Cartesian grid (Nx × Ny × Nz, side lengths 2*Xmax etc.).
2. Initial conditions — the initial wavefunction (IC_vel at the high level, IC at the low level) and V (external potential) are evaluated on the grid.
3. Baryon setup — depending on baryon_mode, baryons are placed on the mesh, sampled as static N-body particles, or co-evolved as dynamic N-body particles.
4. IC potential — the gravitational potential $\Phi_{\rm total}$ at $t = 0$ is computed (FFT for periodic; an isolated tree solve for vacuum in the current driver).
5. Visualization init — observables for Realtime and the per-snapshot update are wired up.
6. KDK main loop — apply_kick_step! → apply_drift_step! for Nt steps.
7. Outputs — per-step and final files (CSV for properties, JLD2 for $\psi$/$\Phi$, PNG/MP4 for visualization).

Key Parameters

Grid and resolution

Initial conditions

Baryons

baryon_mode values:

• :ignored — no baryons. Use this explicitly for dark-matter-only halos and for presets that do not model baryons; the package default is :mesh.
• :mesh — baryons live on the same FFT grid; potential computed by the same FFT Poisson solver as the DM. Best for spatially extended baryons (e.g. dSph, ETG).
• :particles_static — baryons are sampled once as N-body particles and held fixed; their gravitational potential is added to $\Phi_{\rm total}$ at every step.
• :particles_dynamic — baryons are co-evolved with the DM; their potential is re-computed at every step and the DM acceleration is interpolated back to the particles. Most expensive, most realistic.

Tidal field

Boundary conditions

Timestep

Performance

Visualization

Output

Minimal Working Examples

1. Smoke test (CPU)

using WaveDM

simulate_waveDM(;
    model = :MW,
    baryon_mode = :ignored,
    Xmax = 2.0,
    Tmax = 0.001,
    Nt   = 3,
    Nx = Ny = Nz = 64,
    gpu = false,
    Realtime = false,
    unicode_plot = true,
    outputdir = mktempdir(),
)

2. GPU production run

using WaveDM
using Unitful

simulate_waveDM(;
    model = :MW,
    Xmax = 5.0,
    Tmax = 5.0u"Gyr",
    Nx = Ny = Nz = 256,
    Nt = 5000,
    autoset_timestep = true,
    gpu = true,
    Realtime = true,
    plot_virial = true,
    title = "MW_256_T5Gyr",
    outputdir = "results/MW_256_T5Gyr",
)

3. Crater II ultra-faint dwarf

using WaveDM
using Unitful

simulate_waveDM(;
    model       = :dwarf_UFDs,
    Galaxy_id   = 6,                 # Crater II
    baryon_mode = :ignored,          # Crater II is DM-dominated
    mₐ          = 50 * 1e-22 * 1.783e-36u"kg",  # m₂₂ = 50
    V           = (x, y, z, ψ) -> 0.0,
    Nx          = 384,
    Xmax        = 20u"kpc",
    Tmax        = 6.0u"Gyr",
    autoset_timestep = true,
    absorb_coeff = 10,
    bulk_size   = 4,                 # value used in the paper
    gpu         = true,
    Realtime    = true,
    title       = "CraterII",
)

4. Crater II with MW tidal field

using WaveDM
using Unitful

simulate_waveDM(;
    model       = :dwarf_UFDs,
    Galaxy_id   = 6,
    baryon_mode = :ignored,
    mₐ          = 50 * 1e-22 * 1.783e-36u"kg",
    V           = (x, y, z, ψ) -> 0.0,
    Nx          = 384,
    Xmax        = 20u"kpc",
    Tmax        = 6.0u"Gyr",
    autoset_timestep = true,
    absorb_coeff = 10,
    MW_tidal_field = true,
    df_traj     = df_traj,          # from the lookback step (see Examples §4)
    tidal_lookback_time = 6.0u"Gyr",
    LMC_tidal_field = true,         # include LMC perturbation
    gpu         = true,
    title       = "CraterII_MW_LMC",
)

5. Baryons as dynamic N-body particles

using WaveDM
using AstroIC
using PhysicalParticles, Unitful

baryons = generate(AstroIC.ExponentialDisk(; ...))  # see AstroIC.jl

simulate_waveDM(;
    Nx = 256,
    Xmax = 5.0,
    Tmax = 2.0u"Gyr",
    autoset_timestep = true,
    baryon_mode       = :particles_dynamic,
    baryon_particles  = baryons,
    GravitySolver     = Tree(),
    SofteningLength   = 0.1u"kpc",           # code default
    gpu = true,
)

6. Distributed run

using Distributed
addprocs(8)
@everywhere using WaveDM

simulate_waveDM(;
    Nx = Ny = Nz = 256,
    Tmax = 2.0u"Gyr",
    autoset_timestep = true,
    distributed_memory = true,         # use DArray for grid arrays
    gpu = false,                       # GPU and distributed are not orthogonal but mixing
                                      # them is an advanced use case
)

Return Values

ψ, fig, figMOND, chi2RC, dfProp, dfAcc, averaged_ψ2 = simulate_waveDM(...)

Output Files

The simulation writes several file types into outputdir:

• $(title), $(suffix) - Overview IC.png — initial snapshot
• $(title), $(suffix) - Overview.mp4 — animation of the entire run (if Realtime = true)
• $(title), $(suffix) - Overview Prop.png — final-state overview
• $(title), $(suffix) - RC.png / RAR.png — rotation curve / radial acceleration relation
• $(title), $(suffix) - Prop best fit.jld2 — best-fit snapshot (if extract_dwarf_granule = true)
• $(title), $(suffix) - acc.csv / acc best fit.csv — radial acceleration data
• $(title), $(suffix) - best_fit_time.csv — best-fit simulation time
• $(title), $(suffix) - timestep_*.jld2 — per-snapshot ψ / Φ (if average_snapshots = true)

The suffix is a human-readable summary of the run, e.g. Nx(384), Xmax(20.00), Nt(377), Tmax(6.00), DM_m(1.00).

Performance Considerations

• GPU vs CPU. gpu = true is most useful for grids $\gtrsim 256^3$. For tiny problems the kernel-launch overhead dominates and CPU is faster. On a 12 GB device, $N \gtrsim 512$ becomes memory-limited.
• Distributed memory. The distributed FFT Poisson solver is not implemented yet; use distributed workers for IC generation and N-body work only.
• Grid size. Powers of two (64, 128, 256, …) maximize FFT performance. The code also runs on non-power-of-two grids.
• Timestep. autoset_timestep = true keeps $\Delta t$ near the CFL limit without manual tuning.
• StepsBetweenSnapshots. Larger values skip more visualization updates per step, which is useful when the plot is the bottleneck.
• Memory. On GPU, $N^3$ complex arrays dominate. A $512^3$ single-precision grid is $\sim 2$ GiB; a $1024^3$ grid is $\sim 16$ GiB.
• KDK_flag. Leave at true. Setting it to false reduces the scheme to first-order Lie splitting and is provided only for code-experimentation purposes.

Related Modules

• Configs — configuration structures
• KDK — the kick–drift–kick implementation