implode

IMPLODE stays for "planetesImal forMation via Pebble cLOuD collapsE"

This is a Fortran code to calculate collapse of a bound pebble cloud driven by inelastic collisions of pebbles and gas drag written by Joanna Drążkowska and Rico Visser.

The representative particle approach is used throughout the code. The motion of particles is integrated with a Runge-Kutta Fehlberg variable step scheme. A Monte Carlo algorithm is used to solve the collisional evolution.

Code usage

Setup

Compile with make

Run with ./implode parameters where parameters is a parameter file. A couple of example files are included in the repository, in particular virialtest.par is a setup with no gas and no collisions considered where the particles are initialized with random velocities corresponding to virial equilibrium, settlingtest.par is a test in which there is gas but no collisions and the particles are not given any initial velocity dispersion, and setup.par corresponds to one of the models presented in Visser et al. paper.

WARNING: the code is OPENMP parallel and thus will by default use all availble resources. To set the number of threads use set OMP_NUM_THREADS=n or export OMP_NUM_THREADS=n depending on your operating system.

Before re-compiling, clean the directory with make clean - WARNING: it will remove the output if you do not re-name it. When you just change the parameter file and not the source code, there is no need to re-compile. The simulation will not run if the output is overwritten.

Parameters that can be set in the parameter file and changed without re-compiling the code are:

gas is gas drag taken into account? (logical)
collisions are collisions taken into account? (logical)
nr_parts number of representative particles (integer)
nr_parts_per_zone minimum number of representative particles per zone (integer, note that doe to Monte Carlo method convergence it should not be below 200)
radial_distance_AU radial distance of the clump (float, [AU])
gas_surfde_cgs surface density of the gas (at the clump location, float [g cm-2])
final_radius_km final radius of planetesimal (float, [km])
max_time_yr maximum time of the simulation (float, [years])
delta_t_omegaK basic timestep for advection solver (float, [1/OmegaK the local Keplerian frequency])
softening_length_Rsolid softening length for the advection solver (float, [in the units of the final planetesimal radius])
error_tolerance relative error tolerance for the advection solver [float]
init_min_size minimum size of pebbles in the initial size distribution (float, [cm])
init_max_size maximum size of pebbles in the initial size distribution (float, [cm])
initial_size_distr is the initial size distribution considered? (logical, if false, all grains start at init_max_size)
size_distr_kappa slope of the initial size distribution (float, default is 0.1666666667 corresponding to the MRN distribution)
monomer_size_cm monomer size (sets the minimum grain size in the case of fragmentation, float, [cm])
material_density_cgs material denisty of pebbles and the final planetesimal (float, [g cm-3])
bouncing is bouncing included as a collision outcome? (logical)
COR coefficient of restitution for bouncing collisions (float)
fragmentation is fragmentation included as a collision outcome? (logical)
dmmax collisional acceleration parameter, minimum relative change of the target mass in a sticking collision (float)
initial_vel_disp_factor multiplier of the initial velocities (float, 1 corresponds to the virial equilibrium)
init_free_fall initialize the radial velocity component with free-fall/terminal speed velocity? (logical)
output_freq number of timesteps between outputs (integer)
output_dt_year maximum time between outputs (float [years])
screen_out write extended screen output? (logical)
full_output should the full output be saved? (in =false, only the most recent full output is saved in the output...dat file, see below)

Output

The full output is stored in output...dat file and it is written every output_freq timesteps, the frequency can be set in the parameter file. Additionally, there is a timestep criterion output_dt_year to avoid sparse output, the maximum timestep between outputs can also be set in the parameter file. The output format is:

The first row has 3 columns: these are time in years, physical mass represented by one superparticle (in gram, constant), initial radius of the cloud (corresponding to the Hill radius, in cm) The other rows have nr_parts (=total number of superparticles used) columns each and are:
x position [cm]
y position [cm]
z position [cm]
mass of the physical particle represented by each superparticle, normalized by the initial minimum pebble mass specified in the parameter file
x velocity component [cm/s]
y velocity component [cm/s]
z velocity component [cm/s]
number of collision the superparticle underwent (integer)
is the particle part of the core (not evolving anymore) (logical)
number of bouncing collisions the superparticle underwent (integer)
number of fragmenting collisions the superparticle underwent (integer)

Additionally, there is the energy...dat file, which, for every timestep, saves the 1. time [yrs], 2. total kinetic energy T [erg], 3. potential energy abs(U) [erg], 4. core mass fraction, 5. total number of collisions, 6. total number of bouncing collisions, 7. total number of fragmentation collisions, 8. optical depth of active particles containing the inner 1% of mass, 9. settling threshold radius location [Hill radii].

The sizedistr-init...dat contains the initial size distribution of particles and sizedistr...dat contains the final size distribution. The columns are 1. mass in the units of monomer mass (specified in the parameter file), 2. size [cm], 3. the distribution function f(m)m^2.

License

This code is licensed under the GNU General Public License v3. Among other things, this means that if you modified this code for a publication, the original code needs to be cited (see below), the changes need to be stated, and the new code needs to be open sourced under the same license.

Credits

Please cite Visser, Drążkowska, & Dominik (2021) and the Zenodo DOI.