Running multistar

A typical run may look as follows:

from multistar import multi as M
from physconst import YR
m = M('-eight')
m.rund(1*YR)
m.plot3()

rund

dynamic (time step) driver

Arguments

dt

total time to run (seconds)

If not set, it can be computed from (tbegin, and) tend

Optional arguments

nstep

maximum number of integrator steps to run

0: no limit (default)

In case of backups, e.g., too large time steps, these are computed, but not saved as data.

dtd

minimum time step for recording output (s)

tsave

maximum time to save (discard older)

truncate

True: last time step of run is truncated to requested dt (default)
False: code will stop when time >= dt is reached.

runp

Poincare driver (find and output sections)

Based on rund.

from physconst import YR
from multistar import multi as M
m = M('-DI_Herculis')
m.runp(YR,'periapsis',dtmin=1.e-6)
n = M('-DI_Herculis')
n.rund(YR)
plot(n.t, n.rno[0])
plot(m.t, m.rno[0], marker='*', ls='none')

Arguments

Supports same arguments as rund and additionally:

section (str)

section method to use

SLICE

simple coordinate slice (default)

Parameters:

coord (int)

coordinate index in y array to use for section

default: 0

First are the x, y, z coordinates of first orbit (indices 0-2), then of next orbit, …, followed by x, y, z velocities of first orbit, etc., followed by angular momentum components for each object (if present), followed by the orientation of each object (if present).

PERIAPSIS

crossings of radial velocity

Parameters:

orbit (int)

index of orbit

default: 0

sign (float)

sign for crossing

+1: periapsis (default)
-1: apoapsis

PLANE

crossing of arbitrary plane for 3-coordinate

Note

Currently does not work for orientation in quaternion space.

It also seems a bit unclear what exactly to do in spin vector space.

Parameters:

coord (int)

index of coordinate vector in y array

first are norb orbits, then norb velocities, then nstar angular momenta, then orientations

normal (float(3))

normal vector of plane

default: [1, 0, 0]

distance (float)

distance of plane from origin in normal vector direction

default: 0

SYNOD

orbits cross common plane (in same plane as orbital angular momentum vector)

Parameters:

orbit1 (int)

first orbit

default: 0

orbit2 (int)

second orbit

default: 1

mode (int)

determine what angular momentum vector to use

0: average specific angular momentum (default)
1: first orbit
2: second orbit

CUSTOM

select custom section

User code source file can be set using the MULTISTAR_CUSTOM_SECTIONS environment variable. The filename needs to be an absolute path or use ~ for the home directory. When first setting this up or changing the file location, it may be necessary to delete the multistar build directory so meson creates the correct file dependence. See the file fortran/custom_sections.f90 for a template and fortran/sections.f90 for examples.

custom_data

data to be sent to custom section

The data is converted byte-wise to a np.int64 array using np.asarray. For better control, the use should do this conversion/preparation already.

In your custom section, you may convert the data back into any format you wish using Fortran’s transfer function.

cross (int)

+1: cross to negative to positive (default)

0: either sign change of section function

-1: cross to positive to negative

saveall (bool)

save all converged data rund outputs, not just crossing data

default: False

Changed arguments

nstep: number of sections to compute

runs

Step driver (fixed time step of output, may still do sub-cycling)

Arguments

dt

steps size (seconds)

if not set, it can be set from (tbegin,) tend, and nstep

nstep

number of steps (default: 1000)

Common parameters

tmax

maximum wallclock runtime (s)

nsave

maximum number of data to save (discard older)

cont

False

run from start

True

continue, keep data

Ellipis

continue, discard old data

None

(default)

True: if model was run before,
False: otherwise

<number>

time from which to start

(discard old data if an imaginary number is used)

discard:

Ellipsis

discard initial data when continuing

True

discard all data but the last point

implies nsave = 1

False

keep data (default)

dt0

initial time step (s).

default for maximum time step.

dtmin

minimum time step (s).

0: set to dt0
< 0: set to 1. (FORTRAN DEFAULT)

dtmax

maximum time step (s).

0: set to dt0
< 0: set to infinity (1.e99, FORTRAN DEFAULT)

flags

‘run’: to continue using same flags
None: use default flags (from self._flags)
<int>: use these flags

eps

convergence accuracy

current default is 1e-13

token

set some run identifier (used by genetic code)

cutoff

set cutoff matrix for Jacobi coordinates

interact

set object interaction matrix

etadelay

time between ETA outputs (default 10 s)

tbegin

start time of simulation (s)

This is use for sync with star data file

tend

end time of simulation (s)

silent

set verbosity. (default: self.silent set on object generation, with default None)

verb

binary verbosity flags (bits) for Fortran driver

bit constants are defined in module multistar.constants

constant	bit	description
IVERB_INTERACT	0	interact
IVERB_DRIVER	1	driver
IVERB_INTEGRATE	2	integrate
IVERB_ETA	3	eta
IVERB_SHADOW	4	shadow

Default is to disable all flags if silent is True and to enable all flags if silent is False or None. Bits need to be added as 2**bit.

In module multistar.constants, the following values are defined

constant	value	numerical
VERB_ALL	-1	-1
VERB_NONE	0	0
VERB_INTERACT	2**IVERB_INTERACT	1
VERB_DRIVER	2**IVERB_DRIVER	2
VERB_INTEGRATE	2**IVERB_INTEGRATE	4
VERB_ETA	2**IVERB_ETA	8
VERB_SHADOW	2**IVERB_SHADOW	16

Run status

status

numerical run status flag

constant	value	name
STATUS_UNKNOWN	-1	UNKNOWN
STATUS_OK	0	OK
STATUS_COLLIDE	1	COLLIDE
STATUS_ESCAPE	2	ESCAPE
STATUS_INTERACT	3	INTERACT
STATUS_DIVERGE	4	DIVERGE

property Driver.status_name: str: name of status as string

status_stars

array of stars causing termination

If no such exception occurred, the value should be [0,0].

status_component

component index that caused DIVERGE exception

If no such exception occurred, the value should be -1.

Run statistics

property Driver.runtime: total wallclock time (s)

property Driver.steps: total number of steps computed

property Driver.evals: total number of drivative function evaluations

Run API

Driver.rund(dt=None, *, dt0=None, dtmin=None, dtmax=None, dtd=None, nstep=None, nsave=None, tsave=None, eps=None, tmax=None, silent=None, flags='run', token=None, cont=None, cutoff=None, interact=None, truncate=True, iext=0, etadelay=None, discard=False, tbegin=None, tend=None, verb=None)

Fortran dynamic driver

Parameters:

dt – total time aiming to run

if complex, run to imag(dt)

if negative, run backward in time
nstep – maximum number of steps to run

0 == no limit (default)
dtd – minimum time step for recording output
dt0 – initial time step.

Also used as default value for maximum time step.
dtmin – minimum time step.

if 0, set to dt0

if < 0, set to 1. (FORTRAN DEFAULT)
dtmax – maximum time step.

if 0, set to dt0

if < 0, set to infinity (1.e99, FORTRAN DEFAULT)
tmax – maximum wallclock runtime (s)
cont –

False
run from start

True
continue, keep data

Ellipis
continue, discard old data

None (default)

True
if model was run before,

False
otherwise

real value
continue from that time, keep data

imag value
continue from that time, discard old data

Same as using real value and setting discard to True.

Note

Times are specificed realtibe to t array, not relative to absolute time, t + time.
discard –

Ellipsis
discard initial data when continuing

True

discard all data but the last point
implies nsave = 1

False
keep data (default)
nsave – maximum number of data to save
tsave – maximum time to save
flags –

"run"
continue using same flags (default)

None
use default flags (from self._flags)

int
use these flags

str
interpret usual flag text string
verb – verbose flags for driver

Driver.runs(dt=None, nstep=None, *, nsave=None, dt0=None, dtmin=None, dtmax=None, eps=None, tmax=None, silent=None, token=None, flags='run', cont=None, cutoff=None, interact=None, etadelay=None, discard=False, tbegin=None, tend=None, verb=None)

Fortran step driver

Parameters:

tmax – maximum wallclock runtime (s)
cont –

False
run from start

True
ontinue, keep data

Ellipis
continue, discard data

real value
continue from that time, keep data

imag value
continue from that time, discard old data

Same as using real value and setting discard to True.

Note

Times are specificed realtibe to t array, not relative to absolute time, t + time.
discard –

Ellipsis
discard initial data when continuing

True
discard all data but the last point

implies nsave = 1

False
do nothing (default)
nsave – maximum number of data to save

None:
save all nstep models
flags –

"run"
to continue using same flags

None
use default flags (from self._flags)

int
use these flags

Todo

Documnetaion incomplete.

Driver.runp(dt=None, section='SLICE', cross=1, *, dt0=None, dtmin=None, dtmax=None, dtd=None, nstep=None, nsave=None, tsave=None, eps=None, tmax=None, silent=None, flags='run', token=None, cont=None, cutoff=None, interact=None, truncate=True, iext=0, etadelay=None, discard=False, tbegin=None, tend=None, verb=None, saveall=False, **section_args)

Fortran poincare driver

Parameters:

dt – total time aiming to run

if complex, run to imag(dt)

if negative, run backward in time
nstep – maximum number of steps to run

0 == no limit (default)
dtd – minimum time step for recording output
dt0 – initial time step.

Also used as default value for maximum time step.
dtmin – minimum time step.

if 0, set to dt0

if < 0, set to 1. (FORTRAN DEFAULT)
dtmax – maximum time step.

if 0, set to dt0

if < 0, set to infinity (1.e99, FORTRAN DEFAULT)
tmax – maximum wallclock runtime (s)
cont –

False
run from start

True
continue, keep data

Ellipis
continue, discard old data

None (default)

True
if model was run before,

False
otherwise

real value
continue from that time, keep data

imag value
continue from that time, discard old data

Same as using real value and setting discard to True.

Note

Times are specificed realtibe to t array, not relative to absolute time, t + time.
discard –

Ellipsis
discard initial data when continuing

True

discard all data but the last point
implies nsave = 1

False
keep data (default)
nsave – maximum number of data to save
tsave – maximum time to save
flags –

"run"
continue using same flags (default)

None
use default flags (from self._flags)

int
use these flags

str
interpret usual flag text string
verb – verbose flags for driver

Flags

Flags are internally represented as bit field represented in (64 bit) integer number. The Driver module provides a user-friendly interface.

Current flags

property Driver.flags: current flags as string

Driver.get_flags(numeric=False)

Driver.set_flags(flags=None, /, silent=None, return_old_flags=True, **kwargs): set run flags, returns old flags

Driver.get_flag(flag)

Driver.get_flag_bit(bit)

Driver.set_flag(flag=None, /, silent=None, value=True, **kwargs)

Driver.unset_flag(*args, **kwargs)

Flags of last run

For flags used in recent run (what should be used in analysis routines):

property Driver.run_flags: current run flags as string

Driver.get_run_flags(numeric=False)

Driver.get_run_flag(flag)

Driver.get_run_flag_bit(bit)

Metaflags

Access to historic switches and meta flags.

Currently, GR is the only such flag.

Driver.get_meta_flag(flag)

Driver.get_run_meta_flag(flag)

Star Data

Get the star data for the data points in the solution vector.

Driver.get_stardata(silent=None)

return star data for run points

Cached.

Driver.get_stardata_(i=None, t=None, silent=None)

return star data for select points only

Not chached.

To update data file paths after loading a run where data was discarded on save for efficiency:

Driver.update_stardata_path(path: str | Path)

update path for data files, keeping file names

Parameters:: path (star | pathlib.Path) – new path for data files

Driver.load_stardata(path=None)

load star data files, if neccessary (e.g., after loading model)

Parameters:: path (star | pathlib.Path) – new path for data files (optional)

You may also manually overwrite the stardatafiles list of filenames.

When data is loaded on first use, the data hash is checked.

Derivative data

There are two levels of derivative data, simply the derivative vector as used by the integrator, the \(y^{\prime}\), and some more extended data on detailed partial contributions, stored separately. The latter can be quite memory-intensive.

These derivatives are computed from the converged solution data after each run. The full derivative routing is called, but only and exactly for the converged solution vectors resulting from the Burlisch-Stoer algorithm. That is, these exact values may not have not been used in any call to the derivative function during the integration itself.

Note

There is a raw version and a regular version of both cached and uncached functions. Currently the regular version returns the same as the “raw” version.

Todo

In the future, the non-raw version may return a different data format and layout; the raw API should be stable. The issue is that the dimension of some components may vary, e.g., vector components as compared to quaternions.

Simple derivative vector

Driver.get_prime_raw(silent=None, refresh=False, flags=None)

return star prime data for run points

Cached.

Driver.get_prime_raw_(i=None, silent=None, flags=None)

return star prime data for select run points

Not Cached.

Todo

check whether cached data is available

Driver.get_prime(*args, **kwargs): return prime vectors, from cached data

Driver.get_prime_(*args, **kwargs): return prime vectors, not cached

Extended derivative data

Driver.get_derivative_raw(silent=None, refresh=False, flags=None)

return star derivative data for run points

Cached.

Driver.get_derivative_raw_(i=None, silent=None, flags=None)

return star derivative data for select run points

Not Cached.

Todo

check whether cached data is available

Driver.get_derivative(*args, **kwargs): return derivative vectors

Driver.get_derivative_(*args, **kwargs): return derivative vectors

Gravitational wave data

Driver.get_gwdata(silent=None)

return gw data for run points

Cached.

Driver.get_gwdata_(i=None, t=None, silent=None)

return gw data for select points only

Not chached.

Driver.get_gwavedata(silent=None)

return gwave data for run points (QP)

Cached.

Driver.get_gwavedata_(i=None, silent=None)

return gwave data for select points only (QP)

not Cached

Driver.get_gwavefardata(silent=None)

return gwave farfield data for run points (QP)

cached

Driver.get_gwavefardata_(i=None, silent=None)

return gwave farfield data for select points only (QP)

not Cached

Driver.get_gwaveneardata(silent=None)

return gwave nearfield data for run points (QP)

cached

Driver.get_gwaveneardata_(i=None, silent=None)

return gwave nearfield data for select points only (QP)

not Cached

Driver.get_gwaveqqdata(silent=None)

return gwave qq data for run points (QP)

cached

Driver.get_gwaveqqdata_(i=None, silent=None)

return gwave qq data for select points only (QP)

not Cached

History data for GW interaction

Driver.clear_history(): clear history data

Shadows

When shadows have been set up, one can use the index operator to extract regular runs.

m = M("-circle_shadow.toml", update={'scales.div':1})
m.rund(15*YR)
m.plot3(fig=gcf())
m[8].plot3(ax=gca(),ls='--')

status_shadow

shadow image (1-based index) causing termination

If termination is not cause by shadow, the values is 0 to indicate base run. This is also the natural default when not using shadows.

Driver.__getitem__(index): return shadow by index, starting with 0 for original

Driver.__len__()

total number of shadow runs including original run.

Raises an Exception if there is no shadow.

property Driver.shape

shape of shadow runs.

() if no shadow is defined (result is “scalar”) and length-one tuple otherwise (one plus number of shadows).

For Numpy interface compatibility.

property Driver.ndim

dimension of shadow runs.

\(0\) if no shadow is defined (result is “scalar”) and \(1\) otherwise.

For Numpy interface compatibility.

property Driver.shadow: shadow configuration namespace object

property Driver.nshadow: return number of shadows

property Driver.shadow_object: int

return orbit or star that caused divergence

zero-based

for variables r and v, this referes to the corresponding Jacobi coordinate, for j, o, and f to the corresponding star.

property Driver.shadow_variable: str

return variable that caused divergence

May be r, v, j, o, or f.

r: radius
v: velocity
j: angular momentum
o: orientation vectore
f: orientation quaternion (in 4-vector form wxyz)

property Driver.shadow_component: int

return vector component that caused divergence

zero-based

property Driver.shadow_info: dict: return shadow information as dictionary

Run Management

Driver.zerotime(final=False, zero_offset=False)

reset time coordinate to current run time

TODO - fix for use with history

Driver.truncate(t=None, i=None): truncate end of model

Driver.chop(t=None, i=None): chop beginning of model

Driver.reset_time(): reset time to 0 for start of run

Driver.select(t=None, i=None)

keep specific model and set as starting value

Does not reset self.time

Todo

\(dt < 0\)

Driver.slice(t=None, i=None): slice model data by providing ‘slice’ objects for t or i

Driver.copy()

return cleaned copy of self

‘copy’ uses __getstate__, which excudes items not suited for pickling, e.g., plots and figure, as well as chached data.

Driver.get(i=None, t=None): return copy with - data if i or t is slice - just new setup with start index

classmethod Driver.from_result(model, t=None, i=None)

allows to continue run from time t relative to start time

t:: time relative to offset None - from first model (=copy), == 0 == t[0] Ellipsis - from last model (=continue) == t[-1]
i:: Ellipsis: last model
default:: copy from first model

Todo

fix for dt < 0

Driver.reset(): delete results for a fresh start

Driver.new(): return a copy of setup w/o result

Driver.free_cache(): delete cached data

property Driver.version: driver version number

TODO

Todo

re-design ‘flags’ interface, make self-contained object