Construct a coarse mesh to overlap all particles for non-periodic boundary conditions

Compute forces using direct summation method if the particles run out of the non-periodic simulation box

Delete outbound particles if they run out of the non-periodic simulation box

struct DeviceConfig{Device, GA}

• type::Any: Supported: CPU, GPU

• GPU_NumThreads::Int64

• GPU_Algorithm::Any: ::GPUAlgorithm. Supported: AllPairs, Tiled. Default: AllPairs

DeviceConfig(
;
    type,
    GPU_NumThreads,
    GPU_Algorithm
) -> AstroNbodySim.DeviceConfig{CPU, AllPairs}

Direct Summation Method

Finite Differencing Method

Finite Element Method

Fast Fourier Transform

Finite Volume Method

GravityConfig(
;
    ...
) -> AstroNbodySim.GravityConfig{_A, Newton, DS, Vacuum} where _A
GravityConfig(
    units;
    ForceSofteningTable,
    MOND_nuIndex,
    GravityModel,
    OutboundLimiter,
    EnlargeMesh,
    BoundaryCondition,
    sparse,
    model,
    outbound,
    enlarge,
    boundary
) -> AstroNbodySim.GravityConfig{_A, Newton, DS, Vacuum} where _A

struct GravityConfig{Len, GM, OBL, BC}

• ForceSofteningTable::StaticArraysCore.MVector{6}: Gravitational softening lengths for different type of particles

• model::Any: Gravity model. Supported: Newton, MOND1983Milgrom, QUMOND

• MOND_nuIndex::Float64: Index of Mondian interpolation function

• outbound::Any: ::OutboundLimiter. Choose how to handle particles out of the non-periodic simulation box. Supported: - Delete: delete outbound particles - DS: use direct summation method to compute forces - CoarseMesh: Construct a coarse mesh to overlap all particles Default is DS.

• enlarge::Float64: Enlarge the simulation box by enlarge compared to the extent of system. Default is 3.0

• boundary::Any: Boundary condition. Supported: Vacuum, Periodic, Dirichlet. Default is Vacuum

• sparse::Bool: Use sparse matrix to improve performance

Example

GravityConfig(uAstro)
GravityConfig(nothing)
GravityConfig(;ForceSofteningTable=[0.01u"kpc" for i in 1:6])

struct LogInfo

• timers::Dict{String, Int64}: Timer enum names to access timing (continuously starting from 1)

• timing::Vector{UInt64}: Timings defined by timer enums

• analysers::Dict{String, Function}: Analyse on the whole simulation

Examples #TODO

Use uppercase letters to avoid

LogInfo(; timers, analysers) -> LogInfo

Magnetohydrodynamics

Machine Learning

`::OutboundLimiter`. Choose how to handle particles out of the non-periodic simulation box.
Supported:
    - `Delete`: delete outbound particles
    - `DS`: use direct summation method to compute forces 
    - `CoarseMesh`: Construct a coarse mesh to overlap all particles

struct OutputConfig{OT}

• dir::String: Output directory. Default: joinpath(pwd(), "output")

• func::Function: Output function. Default: AstroNbodySim.output

• type::Any: Output Type: traits of gadget2, hdf5, jld2 defined in AstroIO. Default: gadget2()

• format2::Bool: If true, write in format2 mode for Gadget2 format. Default: true

• acc::Bool: If true, output acceleration in snapshots. Default: false

• pot::Bool: If true, output potential in snapshots. Default: false

• SaveRestart::Bool: If true, save restart files regularly. Default: true

• SaveRestartFreq::Float64: Regularly after every fraction of the whole simulation time, save restart files. Default: 1.0

Notes on Keywords

• Keywords are compatible with SimConfig, and can also be overrided by field names of OutputConfig

Example

OutputConfig(; dir = "test", type = jld2(), format2 = false, acc = true, pot = true, SaveRestart = true, SaveRestartFreq = 0.33)

OutputConfig(
;
    OutputDir,
    OutputFunction,
    OutputType,
    format2,
    acc,
    pot,
    SaveRestart,
    SaveRestartFreq,
    dir,
    func,
    type
) -> AstroNbodySim.OutputConfig{gadget2}

mutable struct OutputInfo

• snapshotcount::Int64

• syncflag::Bool

OutputInfo(; snapshotcount) -> OutputInfo

mutable struct PhysicsInfo

• NumForceUpdateSinceLast::Int64

PhysicsInfo() -> PhysicsInfo

Smoothed Particle Hydrodynamics

struct SimConfig{F, U}

• name::String

• author::String

• daytime::Dates.DateTime

• floattype::Any: Numeric type of float numbers

• units::Any: Prefered units in simulation. Default is uAstro. See PhysicalParticles: uAstro, uSI, uGadget2, uCGS. To run without units, use nothing

• ZeroValues::Any: ZeroValue. Pre-constructed zero values for different types to simplify function arguments

• constants::Any: Constant. Physical constants

• loggingmode::Any: Choose how to display status of simulation. Supported: NormalMode, ProgressMode, SilentMode. Default is ProgressMode displaying progress bars

• time::Any: TimeConfig

• output::Any: OutputConfig

• solver::Any: SolverConfig

• grav::Any: GravityConfig

• device::Any: DeviceConfig

Examples

SimConfig(; GravitySolver = Tree())
SimConfig(; device = GPU(), TimeStep = 1.0e-5u"Gyr")
SimConfig(; units = uGadget2)
SimConfig(; units = nothing, ForceSofteningTable = [0.01 for i in 1:6])
SimConfig(; TimeEnd = 1.0u"Gyr", OutputDir = "Test/Dir")

SimConfig(
;
    name,
    author,
    daytime,
    floattype,
    units,
    ZeroValues,
    constants,
    loggingmode,
    TimeBegin,
    TimeEnd,
    TimeBetweenSnapshots,
    TimeStep,
    TimeBase,
    ErrTolTimestep,
    MinStep,
    MaxStep,
    TimeIntegrationAlgorithm,
    OutputDir,
    OutputFunction,
    OutputType,
    format2,
    acc,
    pot,
    SaveRestart,
    SaveRestartFreq,
    GravitySolver,
    ForceSofteningTable,
    MOND_nuIndex,
    GravityModel,
    device,
    GPU_NumThreads,
    GPU_Algorithm,
    OutboundLimiter,
    EnlargeMesh,
    BoundaryCondition,
    sparse
) -> SimConfig{DataType, Vector{Unitful.FreeUnits{N, D, nothing} where {N, D}}}

struct Simulation{D}

• config::Any

• id::Pair{Int64, Int64}

• pids::Vector{Int64}

• simdata::Any

• timeinfo::Any: ::TimeInfo

• outputinfo::Any: ::OutputInfo

• loginfo::Any: ::LogInfo

• physics::Any: ::PhysicsInfo

• stream::Any: ::StreamInfo

• visinfo::Any: ::VisualizationInfo

• buffer::Any: ::Buffer

• bgforce::Vector{Function}

• bgpotential::Vector{Function}

Simulation(
    d;
    floattype,
    units,
    pids,
    timers,
    analysers,
    Realtime,
    RenderTime,
    resolution,
    xlims,
    ylims,
    zlims,
    markersize,
    bgforce,
    bgpotential,
    TreeOpenAngle,
    ErrTolAcc,
    ToptreeAllocFactor,
    MaxTopnode,
    TopnodeFactor,
    TreeAllocFactor,
    MaxTreenode,
    ExtentMargin,
    PeanoBits3D,
    PeanoBits2D,
    epsilon,
    meshmode,
    assignment,
    Nx,
    Ny,
    Nz,
    NG,
    xMin,
    xMax,
    yMin,
    yMax,
    zMin,
    zMax,
    device,
    EnlargeMesh,
    BoundaryCondition,
    kw...
) -> Any

struct SolverConfig{SolverG}

• grav::Any: Gravitational force solver. Supported: DirectSum, Tree, FDM, FFT, ML. Default is DirectSum

SolverConfig(
;
    GravitySolver,
    grav
) -> AstroNbodySim.SolverConfig{DirectSum}

mutable struct StreamInfo

• loggingio::IOStream

• timerio::IOStream

• analyserio::IOStream

struct TimeConfig{Time, Step, I, TIA}

• Begin::Any: Physical time at the beginning of the simulation

• End::Any: Physical time at the end of the simulation

• BetweenSnapshots::Any: Time interval between snapshots

• BetweenSnapshotsInt::Any: Integer time interval between snapshots

• step::Any: Timestep config. Controled by keyword TimeStep. If TimeStep is zero, use AdaptiveTimestep; otherwise use ConstantTimestep

• algorithm::Any: Time integration algorithm. Supported: Euler, Leapfrog

• redshift::Float64: Redshift at start

• scalefactor::Float64: Scale factor at start

Notes on Keywords

• Keywords are compatible with SimConfig, and can also be overrided by field names of TimeConfig
• If timestep is non-zero, use ConstantTimestep; otherwise AdaptiveTimesteps
• TimeBase is the length of integer timeline. Default is 1<<30
• ErrTolTimestep controls the accuracy of adaptive time integration

Examples

TimeConfig(; TimeEnd = 1.0u"Gyr") # Keyword `End = TimeEnd` by default
TimeConfig(; End = 2.0u"Gyr")     # Override keyword from `SimConfig`

TimeConfig(
;
    TimeBegin,
    TimeEnd,
    TimeBetweenSnapshots,
    TimeStep,
    TimeBase,
    ErrTolTimestep,
    MinStep,
    MaxStep,
    TimeIntegrationAlgorithm,
    redshift,
    scalefactor,
    Begin,
    End,
    BetweenSnapshots,
    Step,
    algorithm
) -> Union{AstroNbodySim.TimeConfig{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(Gyr,), 𝐓, nothing}}, AdaptiveTimestep{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(Gyr,), 𝐓, nothing}}, Float64}, Int64, Leapfrog}, AstroNbodySim.TimeConfig{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(Gyr,), 𝐓, nothing}}, ConstantTimestep{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(Gyr,), 𝐓, nothing}}}, Int64, Leapfrog}}

mutable struct TimeInfo{T<:Number, I<:Integer}

• dt::Number: Time interval between neighbor time steps

• system_time_int::Integer

• system_time_float::Number

• last_system_time_int::Integer

• last_system_time_float::Number

• next_output_time_int::Integer

• next_output_time_float::Number

• redshift::Float64

• scalefactor::Float64

• min_endstep::Integer

• stepcount::Integer

TimeInfo(
    config::SimConfig
) -> TimeInfo{T, Int64} where T<:Number

Peano-Hilbert Space Filling Octree Method

struct TreeSimConfig

Controls accuracy of tree method

• TreeOpenAngle::Float64: If the view angle of tree node that relative to the sink point is too large, open this node. In radian unit. Default is 0.1

• ErrTolAcc::Float64: While opening the tree node, take the last acceleration into account. Improves accuracy at high redshift cosmology. Default is 0.025

mutable struct VisualizationInfo

• progress::ProgressMeter.Progress

• PlotData::Any

• resolution::Any

• fig::Any

• Realtime::Bool

• RenderTime::Float64

• last_plot_time::Float64

• xlims::Any

• ylims::Any

• zlims::Any

• markersize::Float64

VisualizationInfo(
;
    Realtime,
    RenderTime,
    resolution,
    xlims,
    ylims,
    zlims,
    markersize
) -> VisualizationInfo

QUMOND_PDM_density(m::MeshCartesianStatic, ACC0::Number)

Compute ρPDM on the RHS (right hand side) of QUMOND (QUasi-linear MOdified Newtonian Dynamics). Return ρPDM

QUMOND(m::MeshCartesianStatic, ACC0::Number, G::Number)

Apply QUMOND (QUasi-linear MOdified Newtonian Dynamics) formula to accelerations

QUMOND_acc(m::MeshCartesianStatic, ACC0::Number, G::Number)

1. Compute ρ_PDM
2. Solve modified potential on the mesh
3. Compute acceleration by finite differencing the potential

Return acceleration

QUMOND_phi(m::MeshCartesianStatic, ACC0::Number, G::Number)

First compute ρ_PDM, then solve QUMOND (QUasi-linear MOdified Newtonian Dynamics) equation on the mesh. Return modified potential

add_timer(
    sim::Simulation,
    name::String,
    t_start::UInt64,
    t_end::UInt64
) -> Any

Add time t_end - t_start to the named timer. This is useful for procedures distributed over several functions.

add_timer(sim::Simulation, name::String, t::UInt64) -> Any

Add time t to the named timer. This is useful for procedures distributed over several functions.

alter_line(filename::String, match::String, newline::String; all = false)

alter the first line containing match to newline

Keywords

• all::Bool = false : If true, alter all lines containing match to newline

alter_param(filename::String, param::String, value)

alter all lines containing param to param value

apply_background_force(f::Function, data::StructArray)

Apply f to all particles in data by calling data.Acc[i] += f(data[i]). f must return a PVector in the unit of physical acceleration.

begin_timer(sim::Simulation, name::String) -> UInt64

Begin the named timer, save time_ns()` to it

function clear(pids = procs())

Clear distributed memories in AstroNbodySim.registry, AstroNbodySim.PhysicalTrees.registry

compute_OldAcc(data::StructArray)

data.OldAcc[i] = norm(data.Acc[i])

compute_acc(
    m::MeshCartesianStatic,
    dim::Val{3},
    mode::VertexMode
)

Compute accelerations by second-order central differencing of potential on the mesh. Results are stored in acc and return nothing.

function computeforcepairwise(p1::AbstractParticle, p2::AbstractParticle, G::Number, h::Number)

return acc of p1 and p2 in a Tuple

compute_kinetic(p::AbstractParticle) -> Any

Compute kinetic energy of particle: 0.5 * p.Mass * p.Vel * p.Vel

crosstime(R, v)

t_cross = R / v

The typical time needed to cross the system.

dynamicaltime(ρ, G)

t_dyn = sqrt(3 * π / (16 * G * ρ))

The time required to travel halfway across the system.

eccentricity(
    G::Number,
    M::Number,
    v::Number,
    h::Number
) -> Any

Eccentricity of elliptic orbit

• G: gravitational constant
• M: mass of central object
• v: velocity at radius r
• h: specific relative angular momentum at radius r

eccentricity(
    G::Number,
    M::Number,
    r::AbstractPoint,
    v::AbstractPoint
) -> Any

Eccentricity vector of elliptic orbit

• G: gravitational constant
• M: mass of central object
• r: position at specific time
• v: velocity vector at position r

ellipticPeriod(G::Number, M::Number, a::Number) -> Any

Orbital period of elliptic orbit

• G: gravitational constant
• M: mass of central object
• a: length of semi-major axis

ellipticSemiMajor(
    G::Number,
    M::Number,
    r::Number,
    v::Number
) -> Any

Length of semi-major axis of elliptic orbit

• G: gravitational constant
• M: mass of central object
• r: orbit radius at specific time
• v: velocity at radius r

freefalltime(ρ, G)

t_ff = sqrt(3 * π / (32 * G * ρ))

The time it takes a sphere with zero pressure to collapse to a point.

Copy data to CPU

Copy data to CPU

Access data on GPU

hubbletime(H0)

t_hubble = 1/H0

The age of the universe.

interactiontime(R, v, N)

The typical time between two short-range interactions that cause a change in kinetic energy comparable to the intrinsic kinetic energy of the particle.

tinteraction = N * tcross

meandensity(M, R)

ρ = 3*M/(4πR^3)

Mean density of a spherical system.

mond_Milgrom1983(sim::Simulation, data::StructArray)

Apply Milgrom 1983 formula of MOND (MOdified Newtonian Dynamics) to accelerations

next_simulationid()

Produces an incrementing ID that will be used for simulations.

nu(y::Number, n::Number = 2)

Generalized MOND (MOdified Newtonian Dynamics) interpolation function. See also nu1, nu2

nu1(y::Number)

MOND (MOdified Newtonian Dynamics) interpolation function with index = 1. See also nu, nu2

nu1(y::Number)

MOND (MOdified Newtonian Dynamics) interpolation function with index = 2. See also nu, nu1

orbitaltime(ρ, G)

t_orb = sqrt(3π/(Gρ))

The time it takes to complete a (circular) orbit.

pseudoNewtonianAcc(G, M, c, R, n::AbstractPoint)

acc = - G * M / (R - rg)^2 * n

Pseudo-Newtonian acceleration around a compat central object. Diverge at gravity radius r_g.

pseudoNewtonianAcc(G, M, c, R, n::AbstractPoint)

acc = - G * M / (R - r_g(G, M, c))^2 * n

Pseudo-Newtonian acceleration around a compat central object. Diverge at gravity radius r_g.

pseudoNewtonianPotential(G, M, c, R)

pot = - G * M / (R - r_g(G, M, c))

Pseudo-Newtonian potential around a compat central object. Diverge at gravity radius r_g.

pseudoNewtonianPotential(G, M, rg)

pot = - G * M / (R - rg)

Pseudo-Newtonian potential around a compat central object. Diverge at gravity radius r_g.

r_g(G, M, c)

r_g = 2 * G * M / c^2

where G is the gravity constant, M is the mass of the central object, c is light speed.

Schwarzchild radius.

relaxtime(R, v, N)

trelax = N/(10 lnN) tcross

Relaxation time: The time over which the change in kinetic energy due to the long-range collisions has accumulated to a value that is comparable to the intrinsic kinetic energy of the particle.

set_softlen_local!(sim::Simulation, c::Collection, h::Number)

Set softening length of collection c to h.

set_softlen_local!(sim::Simulation, h::Number)

Set all softening lengths to h

softlen(c::Collection, table::MVector)
softlen(p::AbstractParticle, table::MVector)
softlen(p::AbstractParticle, sim::Simulation)
softlen(c::Collection, sim::Simulation)

Retrive softening length from table

suggest_softlen!(sim::Simulation)

Set all softening lengths to suggested value.

suggest_softlen(V::Number, N::Int64)
suggest_softlen(data::Union{Array, StructArray})
suggest_softlen(data::Union{Array, StructArray}, collection::Collection)
suggest_softlen(sim::Simulation)

return recommended softening length

See also suggest_softlen! and set_softlen!

tabstring(s::String, tab::String = "    ")

Retract all lines in string s by tab

total_angular_momentum(sim::Simulation) -> Any

Compute total angular momentum vector of the system

total_energy(sim::Simulation) -> Any

Total energy (kinetic + potential) of sim data

total_kinetic(sim::Simulation) -> Any

Compute kinetic energy of particles on workers and return the sum

total_kinetic(data::StructArray) -> Any

Sum kinetic energy: 0.5 * data.Mass[i] * data.Vel[i] * data.Vel[i]

total_kinetic_local(sim::Simulation) -> Any

Compute and sum kinetic energy of particles in data. Potentials need to be computed in advance.

Return nothing if empty.

total_momentum(sim::Simulation, axis::Symbol) -> Any

Compute momentum of the system in the direction of axis

total_momentum(sim::Simulation) -> Any

Compute total momentum vector of the system

total_potential(data::StructArray) -> Any

Sum potential energy of particles in data. Potentials need to be computed in advance.

Return nothing if empty.

typicalvelocity(G, M, R)

v = sqrt(G*M/R)

Typical velocity of a spherical system.

function writegadget2param(filename::String, ICfilename::String)

Params

• filename : name of param file
• ICfilename : name of initial conditions

write_timing(sim::Simulation) -> Any

Save the timings to file

run(sim::Simulation) -> Any

This function does all the work for you:

1. prepare sim environment and preprocess data
2. evaluate force, initialize timesteps, output
3. mainloop

Base.setproperty!(x::LogInfo, symbol::Symbol, d::Dict)

Modify the memories of Dict rather than modifing the pointer of Dict