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 other 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 can denote double, long double or float, depending on
the required numerical precision. This choice is such that sizeof(struct moonlet) is equal
to a line of processor cache (generally 64 bytes), when typ = double.
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.
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, if any, by running git pull in the same directory.
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.
Decimal points should not be removed or added to this file. That is, integers must stay that way, as well as floating point numbers.
Physical constants
The physical constants relative to the simulation are defined by these lines
#define Rearth 1.0 //Radius of the Earth (or central mass) is the unit of length
#define Mearth 1.0 //Mass of the Earth (or central mass) is the unit of mass
#define G 39.47841760435743 //Gravitational constant is set to 4*pi^2, so that the unit of time is the surface orbital period
#define density 0.1448 //The density of the moonlets (3344 kg/m^3) in Mearth/Rearth^3.
#define Tearth 3.4076 //Earth's sideral period in units of the surface orbital period. Must be larger than 1. Today value is 17.038
By default in Ncorpi\(\mathcal{O}\)N, the units of length, mass and time are the radius of the central body, its mass, and the orbital
period at its surface, respectively.
Initial conditions
The parameter file is also used to determine the initial conditions of the simulation.
#define N_max 10000 //The maximum number of moonlets the simulation can handle
#define N_0 1000 //The initial number of moonlets. Must be less than or equal to N_max.
#define M_0 0.02214 //Total (expected) moonlet mass at t=0
#define t_end 512.0 //Total simulation time (in surface orbital period)
#define time_step 0.0078125 //Timestep of the simulation.
#define output_step 128 //Output occurs every output_step timestep. Matters only if write_to_files_bool is 1
#define radius_stddev 0.57 //Standard deviation of moonlet's radii at t=0 (drawn uniformly), in units of the mean radius. Must be less than 1/sqrt(3) to prevent negative radius
#define eccentricity_min 0.0 //Minimal eccentricity for a moonlet at t=0
#define eccentricity_max 0.2 //Maximal eccentricity for a moonlet at t=0
#define sma_min 2.9 //Minimal semi-major axis for a moonlet at t=0
#define sma_max 14.0 //Maximal semi-major axis for a moonlet at t=0
#define inclination_min 0.0 //Minimal inclination (in radians) for a moonlet at t=0
#define inclination_max 0.174533 //Maximal inclination (in radians) for a moonlet at t=0
#define low_dumping_threshold 2.0 //Moonlets falling below this threshold (in Earth radii) are dumped from the simulation (collision with the Earth or disruption by tidal forces)
#define high_dumping_threshold 200.0 //Moonlets going beyond this threshold (in Earth radii) are dumped from the simulation (assumed unbounded)
#define max_ids_per_node 173 //The maximum number of ids in each node of the unrolled linked lists (chains). Choose such that sizeof(struct chain) be a multiple of the cache line
#define softening_parameter 0.0 //The softening parameter for mutual gravitational interations, in units of the sum of the radii.
#define seed 778345128 //The seed used for random number generation. Unimportant if seed_bool is 0.
#define switch_to_brute_force 512 //Threshold for N below which the program switches to the brute-force method for mutual interactions. Unimportant if brute_force_bool is 1
The orbital elements are drawn randomly from a uniform distribution whose bounds are defined by these lines.
The moonlet's radii are drawn uniformy from the range \(\left[\bar{R}\left(1-\sqrt{3}\sigma\right),\bar{R}\left(1+\sqrt{3}\sigma\right)\right]\),
where
- \(\bar{R}^3=\frac{3M_0}{4\pi\rho\left(1+3\sigma^2\right)}\)
- \(\sigma\) is
radius_stddev
- \(M_0\) is
M_0/N_0
- \(\rho\) is
density
Note that if you require non-random initial conditions, you can achieve that by modifying the function populate in
src/structure.c
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
#define seed_bool 0 //Determines if the seed for random number generation is chosen by the user. If seed_bool is 0, then the seed is the number of seconds since 01/01/1970
#define J2_bool 1 //Determines if the contribution from the symmetrical equatorial bulge is taken into account in the simulation
#define Sun_bool 0 //Determines if the perturbations from the Sun are taken into account in the simulation
#define collision_bool 1 //Determines if the moonlets are able to collide
#define mutual_bool 1 //Determines if there are mutual gravitational interactions between the moonlets.
/******** Boolean relative to conservation of the total momentum or total angular momentum ********/
#define tam_bool 0 //Determines if the total angular momentum should be conserved upon merging or fragmenting impact. If tam_bool is 0 then the total momentum is conserved
/******** Booleans relative to collision resolution. Exactly one of them must be 1 ********/
#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 and according to the fragmentation model of Ncorpi\(\mathcal{O}\)N.
/******** Booleans relative to the treatment of mutual interactions (collisions and self-gravity). Exactly one of them must be 1. ********/
#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. (Often the best speed/accuracy compromise for large N)
#define standard_tree_bool 0 //Determines if a standard tree code should be used for mutual interactions. (Significantly slower than falcON for the same accuracy).
#define mesh_bool 0 //Determines if the mesh algorithm should be used for mutual interactions. (Fastest but gravity with neighbours and three largest moonlets only).
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.
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
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.
Collisions
The behaviour of Ncorpi\(\mathcal{O}\)N's fragmentation model is determined by these lines of the parameter file
#define beta_slope 2.6470588235294118//Slope of the power law for fragment size distribution in Leinhardt and Stewart (2012). Must be such that N_tilde is integer.
#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
Running a simulation
If write_to_files_bool is 1, then Ncorpi\(\mathcal{O}\)N writes to output files, at a frequency determined by
output_step. The path towards the simulation files is determined by the string path of the parameter file, for
example
#define path "./ncorpiON/"
The file ./ncorpiON/e.txt would then contain the eccentricities. Each line corresponds to an output step and has as many
columns as the number \(N\) of particles the simulation contained when the output occured.
Compilation
Once the parameter file is filled, Ncorpi\(\mathcal{O}\)N is compiled by running
make clean && make && make clean
Run
The simulation is run with the command
./ncorpion
Making animations
Ncorpi\(\mathcal{O}\)N comes with an add-on implemented in python3 that can makes images out of the simulations files. Then,
ffmpeg can be used to animate these images.
Image creation
For this part, you will need to have python3 installed with the libraries matplotlib and numpy.
You can create the images by running
mkdir -p ./ncorpiON/gif && ./image_creation.sh number_of_images ./ncorpiON/
where number_of_images is t_end/(time_step*output_step), which is \(512\) in our case.
By default, image creation is parallelized on 8 processor cores, but this can be changed by updating the file src/image_creation.sh.
The steps compilation+run+image creation can be done all at once with the command
./ncorpion.sh number_of_images ./ncorpion/
Image assembly
Once the images are produced, you can make an animation out of them by running the command
./ncorpion_animation.sh ./ncorpion/
which calls ffmpeg. This will create the files ./ncorpion/gif/ncorpion.gif and
./ncorpion/gif/ncorpion.mp4.
In our example, we end up with this nice gif
