Structure
Ncorpi\(\mathcal{O}\)N is entirely written in the C language. Since it was primarily developed to study the formation of the Moon, the particles are
referred to as moonlets, even though all types of systems can be efficiently simulated with Ncorpi\(\mathcal{O}\)N, such as protoplanetary disks or
planetary rings.
Moonlet structure
The particles in Ncorpi\(\mathcal{O}\)N are stored in an array of struct moonlet. This array of size N_max is allocated once and for all
at the beginning of the simulation. The integer N_max is chosen by the user and denotes the maximum number of moonlets that the simulation
can handle. An error is returned if the number of particles ever exceeds N_max. The moonlet structure contains the cartesian
coordinates of a particle and is defined as
struct moonlet {
typ x;
typ y;
typ z;
typ vx;
typ vy;
typ vz;
typ mass;
typ radius;
};
where the reserved keyword typ has been defined as double.
This choice is such that sizeof(struct moonlet) is equal to a line of processor cache (generally 64 bytes).
Mutual interactions
Ncorpi\(\mathcal{O}\)N can handle mutual interactions between the particles in four different ways :
- Brute-forcely : Running time proportional to \(N^2\)
- With the standard Barnes-Hut tree code : Running time proportional to \(N \ln N\)
- With a mesh-grid based algorithm : Running time proportional to \(N\)
- With the fast multipole method Falcon : Running time proportional to \(N\)
The brute-force method refers to the straightforward algorithm that goes over all \(N(N-1)/2\) pairs of particles. It does not require
any particular structure and is fastest for \(N < 100\). It however becomes prohibitively slow for large values of \(N\).
Mesh-based algorithm
This algorithm uses a mesh, or 3D-grid, to treat mutual interactions between the moonlets. Collisions are searched, and mutual gravity is
computed, only between particles that are neighbours (a neighbours of a particle is a particle in the same cell or in an adjacent cell).
The user must choose the expected number of neighbours per particle, and the mesh-size (sidelength of a cell) is updated at every timestep
to meet the user's requirement.
The cells of the mesh are stored in an array. With \(1000\) cells per dimension, this method requires an array
of size \(1000^3\). In order not to allocate a too large chunk of memory, the array contains pointers towards the cells,
instead of the cells themselves, as pointers are much lighter than cells. A cell is just a chain (also called unrolled linked list) of the ids
(in the moonlet array defined above) of the moonlets it contains. It is defined as the structure
struct chain {
int ids[max_ids_per_node]; //The array of ids in that node of the chain
struct chain * queue; //A pointer towards the next node of the chain
int how_many; //How many ids are in this node
};
The mesh is a hash table that remembers what ids are in a given cell. It is an array of pointers towards chains of ids,
and has the type struct chain **. Most cells of the mesh are empty, and in that case, the hash table contains a
NULL pointer.
FalcON method and Barnes-Hut tree code
Both the Barnes-Hut code and falcON fast multipole method uses an octree (renamed boxdot tree). It is defined in C as the structure
struct boxdot {
struct boxdot * oct[8]; //The eight octants, or children of that cell.
struct chain * dots; //A chain containing the ids of the moonlets in that cell.
typ corner[3]; //Top-left-front corner of the node.
typ sidelength; //Sidelength of the node (half that of its parent).
int id; //The unique id of that cell, determined by its Hilbert-Peano order : index in the flattree
int how_many; //The number of moonlets in that cell.
int rotation; //A rotation used for determination of the Hilbert-Peano digit
int level; //The level of that node in the boxdot tree
};
This data structure is convenient for the construction of the octree, but its poor locality (children are not
next to each other in memory) leads to frequent cache misses, and travelling in the tree requires multiple time-consuming pointer dereferences.
Therefore, it is used only for the construction, and falcON and the Barnes-Hut code are performed on a regular array that I name flattree.
The flattree is merely an array of nodes (or cells) of the tree. A node is defined by the data structure
struct node {
typ com[3]; //When treating self-gravity : Center of mass and expansion center. When treating collisions : Average position of the moonlets in the node
typ C1[3]; //Acceleration, or first order tensor field of interactions
typ center[3]; //Center of the node
typ sidelength; //Sidelength of the node
typ M0; //When treating self-gravity : Mass of the cell or zeroth multipole moment. When treating collisions : max(R_i + v_i*timestep)
typ r_max; //Convergence radius of the Taylor expansion
typ r_crit; //When treating self-gravity : r_max/theta. When treating collisions : r_max + M0
int * dots; //Array of ids of moonlets contained in that node
int idParent; //Index of the parent node
int idFirstChild; //Index of the first child (Child with smallest unique id)
int how_many_children; //The number of children that node has
int how_many_dots; //The number of moonlets that node contains
};
and contains relevant informations to compute mutual gravity and to search for collisions. The multipole moments \(\mathbf{M}^{(2)}\) to
\(\mathbf{M}^{(5)}\) and the tensor fields of interactions \(\mathbf{C}^{(2)}\) to \(\mathbf{C}^{(6)}\) are stored in separate arrays, indexed identically
to the flattree. The index of a node in the flattree
is given by its Hilbert-Peano order, derived from a 3D version of Hilbert's 1891 space filling curve. This ensures maximum locality.
See Ncorpi\(\mathcal{O}\)N's paper for details on the Hilbert-Peano order and explanations on how these two methods work.
When building the octree, a two to three speed-up factor is achieved by adding the moonlets
in the Hilbert order of the previous timestep. Initially, only the root cell exists and all of its children are NULL.
It is centered on the central mass
and has a sidelength root_sidelength. The user defines the quantity subdivision_threshold and I proceed in the
following manner to add a particle to a cell struct boxdot * Boxdot:
if (Boxdot -> how_many < subdivision_threshold){
if (Boxdot -> how_many == 0){
Allocate memory to the cell;
}
Boxdot -> how_many++;
}
else{
Add the id of the moonlet recursively to the cell's children;
if (Boxdot -> how_many == subdivision_threshold){
Add all the moonlet's ids of the cell to its children;
}
Boxdot -> how_many++;
}
Add the id of the moonlet to the chain Boxdot -> dots;
In order to build the octree, this sequence of instructions is repeated as long as there are moonlets left to add. It has
a time-complexity \(\mathcal{O}(\ln N)\), hence an overall complexity \(\mathcal{O}(N \ln N)\) for the whole building procedure.
In practise, FalcON turns out to be 10 to 15 times faster than the standard Barnes-Hut code for large N.
Set Up
Installation
Ncorpi\(\mathcal{O}\)N is an open-source software distributed with the GNU general public license. It can be obtained from the following
github repository. Navigate to the directory where you would like
Ncorpi\(\mathcal{O}\)N to be installed and run the command
git clone git@github.com:Jeremycouturier/NcorpiON.git
Later on, you can download updates of Ncorpi\(\mathcal{O}\)N by running git pull in the same directory.
Alternatively, if you don't have git installed and don't want to install it, you can simply download the repository by hand
from github.
Running a simulation
Once your simulation is ready, you can compile Ncorpi\(\mathcal{O}\)N and launch the simulation all with one command:
./ncorpion.sh
This command must be run from within the src directory of Ncorpi\(\mathcal{O}\)N.
Parameter file
The parameter file located at src/parameters.h contains all the data needed to run a particular simulation.
The user uninterested in modifying Ncorpi\(\mathcal{O}\)N won't need to touch any other file than this one and the makefile. This file is used to set up the simulation
and to define the behavior of Ncorpi\(\mathcal{O}\)N.
Decimal points should not be removed or added to this file. That is, integers must stay that way, as well as floating point numbers.
Booleans
Boolean values are set to determine how Ncorpi\(\mathcal{O}\)N should behave
#define write_to_files_bool 1 //Determines if the simulation writes to output files. Set to 0 to run speed tests, or if you are satisfied with what is displayed in the terminal
//You can also set this boolean to 0 if you only want to 3D visualize the simulation in your browser.
#define make_animation_bool 1 //Determines if animations of the simulation are produced. write_to_files_bool and write_elliptic_bool must both be set to 1
#define write_cartesian_bool 0 //Determines if the cartesian elements x, y, z, vx, vy, vz should be output. Unimportant if write_to_files_bool is 0. Output in simulation's units
#define write_elliptic_bool 1 //Determines if the elliptic elements a, lambda, k, h, q, p should be output. Unimportant if write_to_files_bool is 0. Output in simulation's units
//If write_to_files_bool is 1 but both write_cartesian_bool and write_elliptic_bool are 0, then only radius.txt, mass.txt and stat.txt are output
#define central_mass_bool 1 //Determines if there is a central mass. If 0, none of the bodies in the simulation play any particular role. If 1, the central body is initially
//at (x,y,z,vx,vy,vz) = {0.0} and treated independently. Should be set to 1 if one body is very massive and if brute_force_bool is 0, so gravity
//with that body is computed directly and not from a tree or a mesh. If 1, some physical effects can be considered as well(see below).
#define reduce_to_COM_bool 1 //Determines if the position and speed of the center of mass of the system are cancelled by a translation before simulating. If 0, then the center
//of mass will drift linearly with time, which could be a problem due to the limit of floating-point representation. Setting to 1 is recommended.
//Note that the position and speed of the central body (if any) will slightly deviate from 0 after the reduction of the center of mass.
#define random_initial_bool 1 //Determines if the initial conditions are drawn randomly between bounds defined below. If 0, the initial conditions are read from the file init.txt
//located at the path defined above, with one line per body and 8 columns, the first six for the initial conditions and the last two for mass and
//radius (in simulation unit). If central_mass_bool is 1 and random_initial_bool is 0, do not include the central body in init.txt.
#define initial_cartesian_bool 1 //Determines if the initial conditions are given in cartesian coordinates (X,Y,Z,vX,vY,vZ). If 0, then the initial conditions are (semi-major axis,
//eccentricity, inclination, true anomaly, argument of periapsis, longitude of ascending node). Unimportant if random_initial_bool is 1.
//The initial conditions in the file init.txt must be given in simulation's units and radians.
#define seed_bool 1 //Determines if the seed for random number generation is chosen by the user. If seed_bool is 0, the seed is the number of seconds since 01/01/1970
#define one_collision_only_bool 0 //Determines if bodies are only allowed to collide once per timestep. If 0, there is no restriction on the number of collisions a body can experience
//during a timestep. Setting first to 1 and then to 0 is a good way to know if the timestep is adapted to the bodies' mean free path.
#define openGL_bool 0 //Determines if a 3D real-time visualization of the simulation is enabled. You must match openGL_bool to the same value in the makefile
#define resume_simulation_bool 0 //Determines if, at the end of the simulation, NcorpiON generates a file named init.txt that can be used to resume the simulation. The file init.txt
//is stored at the path indicated above. To resume the simulation, you need to set random_initial_bool to 0, initial_cartesian_bool to 1, and N_0 to
//the number of lines of init.txt. If init.txt already exists in path pth, it will be overwritten. Simulation's variables should be updated.
/******** Booleans relative to interactions with the central mass or a distant object. Set to 0 if central_mass_bool is 0 ********/
#define J2_bool 1 //Determines if the contribution from the J2 is taken into account in the simulation. The (x,y) plane of the simulation must be the equatorial plane
#define Sun_bool 0 //Determines if the perturbations from a distant object that the system orbits (or is orbited by) are taken into account in the simulation
#define central_tides_bool 0 //Determines if orbiting bodies raise tides on the central body. The tidal model used by NcorpiON is the constant timelag model
#define inner_fluid_disk_bool 1 //Determines if there is an inner fluid disk (disk of liquid material below the Roche radius from which bodies spawn). See Salmon & Canup 2012
//If set to 1, its mass is added to that of the central body when computing gravitational interactions and when preserving the total momemtum
/******** Booleans relative to mutual interactions between the bodies ********/
#define collision_bool 1 //Determines if the bodies are able to collide. If 0, the bodies pass through each other and you should set a non-zero value for softening_parameter
#define mutual_bool 1 //Determines if there is mutual gravity between the bodies. If central_mass_bool is 1, the bodies interact with the central mass even if set to 0.
/******** Booleans relative to the treatment of mutual interactions (collisions and self-gravity). Exactly one of them must be 1. Unimportant if mutual_bool is 0 ********/
#define brute_force_bool 0 //Determines if a brute-force method should be used for mutual interactions.
#define falcON_bool 1 //Determines if falcON algorithm should be used for mutual interactions. (Best speed/accuracy compromise for large N + preserve total momentum)
#define standard_tree_bool 0 //Determines if a standard tree code should be used for mutual interactions. (Much slower than falcON for the same accuracy, momentum not preserved).
#define mesh_bool 0 //Determines if the mesh algorithm should be used for mutual interactions. (Fastest but gravity with neighbours and three largest bodies only).
/******** Booleans relative to collision resolution. Exactly one of them must be 1. Unimportant if collision_bool is 0 ********/
#define elastic_collision_bool 0 //Determines if the collisions are all elastic.
#define inelastic_collision_bool 0 //Determines if the collisions are all inelastic with inelasticity parameter collision_parameter.
#define instant_merger_bool 0 //Determines if the collisions all result in a merger.
#define fragmentation_bool 1 //Determines if the outcome of the collision should be either a merge, a partial fragmentation, a full fragmentation, or a catastrophic disruption,
//depending on the relative velocity. Set to 1 to use the fragmentation model of NcorpiON.
You can choose, among other things, how collisions should be resolved or how mutual interactions should be handled. Here we choose to
handle mutual interactions with FalcON (by far best) and to resolve collisions with the built-in fragmentation model of Ncorpi\(\mathcal{O}\)N.
If you want to 3D visualize the simulation, you need to set the boolean openGL_bool to 1. This boolean also needs to be
set to the same value in the makefile.
System of units
The system of units used by NcorpiON is defined by these three parameters
#define R_unit 1.0 //If central_mass_bool is 1, then radius of the central body. Otherwise, unimportant for the simulation. You define your own unit of length.
#define M_unit 1.0 //If central_mass_bool is 1, then mass of the central body. Otherwise this is the mass for conversions cartesian <-> elliptic
#define G 39.47841760435743 //The gravitational constant. It is here set to 4*pi^2, so that a body of semi-major axis 1 orbiting a mass 1 has a period 1
For example, if you simulate a protoplanetary disk, M_unit would be the mass of the star and R_unit its radius. If you set both of these to
1, this means that the unit of mass (resp. length) of you system is the mass (resp. radius) of the star. If furthermore you
set the gravitational constant to \(4\pi^2\), then the unit of time would be the orbital period of a circular orbit at the surface of the star.
Initial conditions
The initial conditions of the simulation can be defined in two ways. They are either chosen at random by Ncorpi\(\mathcal{O}\)N, in which case all you
have to do is to define their bounds in the parameter file, or you can give them yourself in a file. The file must be named init.txt and be placed
at the input/output path defined in the parameter file. This file has one line per body and 8 columns. Depending if you give your initial conditions in
cartesian or elliptic coordinates, theses columns can be \((x,y,z,v_x,v_y,v_z,m,R)\) or \((a,e,i,\nu,\omega,\Omega,m,R)\).
Tree
If a tree algorithm is used to handle mutual interactions (that is, if either falcON_bool or standard_tree_bool
is 1), then these lines of the parameter file must be updated
#define expansion_order 3 //The order p of the Taylor expansions. This is p = 3 in Dehnen (2002). NcorpiON allows up to p = 6. Minimum is 1 since order 0 yields no acceleration
#define theta_min 0.5 //Minimal value of the tolerance parameter theta. Must be strictly less than 1. Sensible values are 0.2 < theta_min < 0.7
//The precision (and computational time) of the mutual gravity computed by Ncorpion increases with increasing p and decreasing theta_min.
#define subdivision_threshold 17 //A cubic cell is not divided as long as it contains at most that many moonlets. Called s in Dehnen (2002). Must be > 0
//The computational time depends a lot on this threshold. Suggested values are s ~ 26 if p = 3 and s ~ 80 if p = 6, but should be tweaked by the user.
#define root_sidelength 80.0 //Sidelength of the root cell. Must be machine representable. Particles outside of the root cell only feel Earth's gravity and do not affect others
#define level_max 25 //The maximum allowed number of levels in the tree. Root is at level 0 and a cell at level level_max - 1 is never divided.
#define child_multipole_threshold 1 //If number of particles/number of children is at most this threshold then the multipole moments of a cell are computed directly from the particles. Otherwise, they are computed from that of the children
The following lines are specifically needed when falcON is used
#define N_cc_pre 8 //For two cells with N1 and N2 moonlets, if N1N2 < N_cc_pre, then the interaction is computed brute-forcely, no matter if they are well-separated or not
#define N_cc_post 64 //For two cells with N1 and N2 moonlets, if N1N2 < N_cc_post and they are not well-separated, then the interaction is computed brute-forcely
#define N_cs 64 //If N1 < N_cs, the self-interaction of a cell containing N1 moonlets is computed brute-forcely
#define N_cc_collision 16 //This is N_cc_post. For two non well-separated nodes with Na and Nb moonlets, if NaNb < N_cc_collision, then collisions are searched brute-forcely, otherwise, the largest node is subdivised. N_cc_pre is 0 when falcON is used for collision detection.
#define N_cs_collision 12 //For a node with Na moonlets, if Na < N_cs_collision, then collisions are searched brute-forcely, otherwise, the node is subdivised
Whereas these lines are related to the standard Barnes-Hut tree code
#define N_cb_pre 6 //If N1 < N_cb_pre, the interaction between a moonlet and a cell with N1 moonlets is performed brute-forcely, no matter if they are well-separated or not
#define N_cb_post 16 //If N1 < N_cb_post and they are not well-separated, the interaction between a moonlet and a cell with N1 moonlets is performed brute-forcely
#define N_cb_collision 16 //This is N_cb_post. For a moonlet not well-separated from a node with Na moonlets, if Na < N_cb_collision, then collisions are searched brute-forcely, otherwise, the node is subdivised. N_cb_pre is 0 when the standard tree code is used for collision detection.
Some of theses values are not straightforward to understand, but details are given in functions TreeWalk (for falcON) and StandardTree of Sects. 3.3.4 and 3.3.5 of Ncorpi\(\mathcal{O}\)N's paper.
Mesh
If the mesh \(\mathcal{O}(N)\) algorithm is used to handle mutual interactions, then theses lines of the parameter file must be updated
#define collision_cube_min 80.0 //Minimal sidelength of the mesh-grid centered on the Earth. The mesh-size will never be less than collision_cube_min/collision_cube_cells
#define collision_cube_cells 1024 //The number of mesh cells per dimension of the collision cube. Must be even. For 16+ GiB of RAM (resp. 8 or 4 GiB), choose ~1000 (resp. ~800 or ~500)
#define how_many_neighbours 16.0 //The desired expected number of neighbours for a moonlet
Collisions
The behaviour of Ncorpi\(\mathcal{O}\)N's fragmentation model is determined by these lines of the parameter file
#define N_tilde 15 //2*beta_slope/(3 - beta_slope). This is the ratio between the ejected mass and the mass of the second largest fragment -> N° of fragments in the tail
#define mu_parameter 0.55 //The exponent of the impact velocity in the coupling parameter. See Table 3 of Housen & Holsapple (2011)
#define C1_parameter 1.5 //A dimensionless parameter of impact theories. See Table 3 of Housen & Holsapple (2011)
#define k_parameter 0.2 //A dimensionless parameter of impact theories. See Table 3 of Housen & Holsapple (2011)
#define frag_threshold 0.000000005 //Ejected mass threshold below which the collision results in a merger. If the ejected mass is above that threshold but the mass of the second largest fragment is below, then the tail is reunited into a single moonlet. If the mass of the second largest fragment is above that threshold then the collision results in a full fragmentation.
#define pq_min_max {-1,3,-1,1} //Extremal integer values for p_k and q_k to determine the position of the tail fragments with respect to the largest fragment. Must define a rectangle containing exactly N_tilde points with integer coordinates. More precisely, if pq_min_max = {a,b,c,d}, then we must have N_tilde = (b-a+1)(d-c+1). See NcorpiON's paper (Sect. 5.3.5) for details
Getting started on Ncorpi\(\mathcal{O}\)N
If you don't know where to begin, an idea is to simply run ./ncorpion.sh from the src directory before doing anything else. This
will launch the default simulation of Ncorpi\(\mathcal{O}\)N, which is just a collisional and fragmenting disk around a central mass.
The simulation will simply display the progress in terminal, but will not output to files, nor produce animations or launch a 3D visualization to your browser.
As such, it is reasonably likely that it will work out of the box on your system.
From there, you can set write_to_files_bool to 1 and define an input/output path so Ncorpi\(\mathcal{O}\)N will actually save the
simulation on files. You will be able to control the frequency of outputs with the parameter output_step. Then you can set make_animation_bool
to 1 so that animations will be produced. Here is the kind of animations that Ncorpi\(\mathcal{O}\)N will output.
Finally, you can set the boolean openGL_bool to 1 (in both the parameter file and the makefile), so that a 3D animation of the simulation
will be happening on your browser at the address http://localhost:1234/ while the simulation is running. Here is what you should
have if you visit this address.
The 3D real-time visualization will need MPI to be installed on your system. Similarly, you will need ffmpeg and python for the 2D animations. Look at the
github repository for installation instructions.
