diff --git a/README.md b/README.md index 75b22f0bc..5db472fec 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# Swiftest +# Swiftest User Manual ### The Purdue University Swiftest Team #### Carlisle Wishard, David Minton, Jennifer Pouplin, Jake Elliott, & Dana Singh @@ -145,7 +145,7 @@ sim.get_parameter(**kwargs) # View the defa sim.set_parameter(**kwargs) # Set any desired simulation parameters sim.write_param(**kwargs) # Write simulation parameters to the param.in sim.save(**kwargs) # Save the simulation initial conditions to init_cond.nc -sim.run(**kwargs) # Run the simulation (leave off if running from the terminal) +sim.run(**kwargs) # Run the simulation (leave off if running from the executable) ``` To read in a set of Swiftest output files using the Swiftest Python package, follow the general script below. For more details on the output files and user options, continue reading this section. @@ -195,11 +195,11 @@ The key word arguments available to the user for the ```add_body``` method are d Once all desired bodies have been added to the Swiftest simulation, the simulation parameters can be accessed and changed using the ```get_parameter``` and ```set_parameter``` methods. The key word arguments available to the user for the ```get_parameter``` and ```set_parameter``` are the same as those described in [simulation_kwargs](README_tables/simulation_kwargs.md). -After all desired parameters have been set, the parameters can be saved to the **param.in** using the ```write_param``` method. The key word arguments available to the user for the ```write_param``` method are described in [write_param_kwargs](README_tables/write_param_kwargs.md). +After all desired parameters have been set, the parameters can be saved to the **param.in** using the ```write_param``` method. The key word arguments available to the user for the ```write_param``` method are described in [write_param_kwargs](README_tables/write_param_kwargs.md). The state of the system can be saved to the initial conditions NetCDF file, **init_cond.nc**, using the ```save``` method. The key word arguments available to the user for the ```save``` method are described in [save_kwargs](README_tables/save_kwargs.md). -Finally, a simulation can be run from the same script in which it is created (or a separate Python script) using the ```run``` method. This is optional as the simulation can also be run from the terminal. More details on running a Swiftest simulation can be found in the section **Running a Swiftest Simulation**. The key word arguments available to the user for the ```run``` method are the same as those described in [simulation_kwargs](README_tables/simulation_kwargs.md). +Finally, a simulation can be run from the same script in which it is created (or a separate Python script) using the ```run``` method. This is optional as the simulation can also be run from an executable. More details on running a Swiftest simulation can be found in the section **Running a Swiftest Simulation**. The key word arguments available to the user for the ```run``` method are the same as those described in [simulation_kwargs](README_tables/simulation_kwargs.md). **ASCII Input Files** Swiftest accepts 4 ASCII input files. All four input files are necessary, however the structure of each input file varies slightly depending on the features and capabilities of the integrator selected. The four input files are as follows: @@ -209,7 +209,7 @@ Swiftest accepts 4 ASCII input files. All four input files are necessary, howeve - **pl.in** - The massive body input file. - **tp.in** - The test particle input file. -The parameter options used in the parameter input file are as described in [param_options](README_tables/param_options.md). +The parameter options used in the parameter input file are as described in [simulation_kwargs](README_tables/simulation_kwargs.md). The **cb.in** includes all central body initial conditions. The structure of the **cb.in** is as follows: @@ -257,7 +257,7 @@ Note that the ID numbers of the test particles are a continuation of the ID numb **Running a Swiftest Simulation** -The input files necessary to successfully run Swiftest should now be generated in the simulation directory. The user is now faced with a second choice: to run a Swiftest simulation from a Python environment or to run it directly from the terminal. Either option is possible with NetCDF format input files, however ASCII input files must be run directly from the terminal. +The input files necessary to successfully run Swiftest should now be generated in the simulation directory. The user is now faced with a second choice: to run a Swiftest simulation from a Python environment or to run it directly from an executable. Either option is possible with NetCDF format input files, however ASCII input files must be run directly from an executable. **Running via Python** @@ -274,9 +274,9 @@ sim = swiftest.Simulation(simdir = "directory_name", read_param=True) sim.run() ``` -**Running via a Terminal** +**Running via an Executable** -To run a Swiftest simulation through the terminal, create a symbolic link to the Swiftest driver from your current directory. +To run a Swiftest simulation through an executable, create a symbolic link to the Swiftest driver from your current directory. ``` $ ln -s ~/PATH/TO/swiftest/bin/swiftest_driver . @@ -295,9 +295,9 @@ Where ```INTEGRATOR``` is your integrator of choice, either ```whm```, ```rmvs`` The number and type of output files generated by Swiftest depends on the input parameters selected and the method through which Swiftest was run. The standard output files are as follows: - **data.nc** - Always generated, the output file containing the information for every body in the system, recorded every ```ISTEP_OUT``` timesteps and written every ```DUMP_CADENCE```. This file can be analyzed using the Swiftest Python package (```sim.data```). - **collisions.log** - The log containing the record of each fragmentation event, including the collisional regime, and the number of the fragments created, only if ```FRAGMENTATION``` is ```YES```, Swiftest SyMBA only. -- **swiftest.log** - A log containing a brief update on the status of the run. Only generated if Swiftest is run through the Python package or through a shell script. If Swiftest is run directly through the terminal, these updates are output directly to the terminal. -- **collisions.nc** - The details of each collision that occurs in a simulation are recorded in a NetCDF file. Only if ```CHK_CLOSE```/```close_encounter_check``` is ```YES```/```True```. This file can be analyzed using the Swiftest Python package (```sim.collisions```). -- **encounters.nc** - The details of each close encounter that occurs in a simulation are recorded in a NetCDF file. Only if ```CHK_CLOSE```/```close_encounter_check``` is ```YES```/```True```. This file can be analyzed using the Swiftest Python package (```sim.encounters```). +- **swiftest.log** - A log containing a brief update on the status of the run. Only generated if Swiftest is run through the Python package or through a shell script. If Swiftest is run through an executable, these updates are output directly to the terminal. +- **collisions.nc** - The details of each collision that occurs in a simulation are recorded in a NetCDF file. Only if ```CHK_CLOSE```/```close_encounter_check``` is ```YES```/```True```. This file can be analyzed using the Swiftest Python package (```sim.collisions```). +- **encounters.nc** - The details of each close encounter that occurs in a simulation are recorded in a NetCDF file. Only if ```CHK_CLOSE```/```close_encounter_check``` is ```YES```/```True```. This file can be analyzed using the Swiftest Python package (```sim.encounters```). - **init_cond.nc** - The initial conditions used to run the simulation. This file can be analyzed using the Swiftest Python package (```sim.init_cond```). - **encounter_check_plpl_timer.log** - The log containing the encounter check timer for each massive body/massive body encounter, only if ```CHK_CLOSE```/```close_encounter_check``` is ```YES```/```True``` and ```ENCOUNTER_CHECK```/```encounter_check_loops``` is ```ADAPTIVE```. - **encounter_check_pltp_time.log** - The log containing the encounter check timer for each massive body/test particle encounter, only if ```CHK_CLOSE```/```close_encounter_check``` is ```YES```/```True``` and ```ENCOUNTER_CHECK```/```encounter_check_loops``` is ```ADAPTIVE```. @@ -324,7 +324,7 @@ The first line includes the simulation time, the fraction of the simulation that **Restarting a Simulation From t $\neq$ 0** -Just like Swiftest allows the user to run a simulation through the terminal or through Python, Swiftest also allows the user to restart a simulation from t $\neq$ 0 in the same two manners. This can be useful in the case of an accidental termination of a simulation, such as through a power outage or computer failure. In many cases, it is also necessary to run a simulation to a new end point, past the original ```TSTOP```. +Just like Swiftest allows the user to run a simulation through an executable or through Python, Swiftest also allows the user to restart a simulation from t $\neq$ 0 in the same two manners. This can be useful in the case of an accidental termination of a simulation, such as through a power outage or computer failure. In many cases, it is also necessary to run a simulation to a new end point, past the original ```TSTOP```. **Restarting via Python** @@ -340,9 +340,9 @@ sim.run() Note that Swiftest will look in the ```/simdata``` subdirectory for the initial conditions by default. You may set a new path to the initial conditions using the ```param_file``` keyword argument. -**Restarting via a Terminal** +**Restarting via an Executable** -Every ```DUMP_CADENCE``` X ```ISTEP_OUT``` timesteps, Swiftest writes all simulation information from memory to the output files. At the same time, Swiftest also writes all simulation information to a new parameter file, titled **param.XXXXXXXXXXXXXXXXXX.in**. To restart a run from a previous parameter file, simply follow the instructions detailed in the **Running via a Terminal** section, replacing ```param.in``` with the name of the parameter file from which you wish to restart. +Every ```DUMP_CADENCE``` X ```ISTEP_OUT``` timesteps, Swiftest writes all simulation information from memory to the output files. At the same time, Swiftest also writes all simulation information to a new parameter file, titled **param.XXXXXXXXXXXXXXXXXX.in**. To restart a run from a previous parameter file, simply follow the instructions detailed in the **Running via an Executable** section, replacing ```param.in``` with the name of the parameter file from which you wish to restart. --- @@ -474,13 +474,19 @@ To process the output file, run the script titled **scattermovie.py **Fragmentation** -This example highlights the functionality of the Fraggle algorithm. It can be found in the ```/swiftest/examples/Fragmentation``` directory. It is intended to be run using the SyMBA integrator. It contains three pre-built collisional test cases: +This example highlights the functionality of the Fraggle algorithm. It can be found in the ```/swiftest/examples/Fragmentation``` directory. It is intended to be run using the SyMBA integrator. It contains 9 pre-built collisional test cases: - A Head-On Disruptive Collision +- An Off-Axis Disruptive Collision +- A Head-On Super-Catastrophic Disruptive Collision - An Off-Axis Super-Catastrophic Disruptive Collision -- A Disruptive Hit and Run Collision +- A Disruptive Hit and Run Collision +- A Pure Hit and Run Collision +- A Merger +- A Merger Crossing the Spin Barrier +- All of the Above -To generate a movie depicting the collision and results of each of the test cases, run the Python script titled **Fragmentation_Movie.py**. +To generate, run, and create a movie depicting the collision, run the Python script titled **Fragmentation_Movie.py**. Please note that this example requires a large amount of memory. For reference, this example was created and run using 4 nodes, each with 256 GB of memory. This amount of computational memory is necessary to generate a smooth movie. In this example, the trajectories of all bodies involved in the collision are saved at every point in the simulation. This is extremely expensive and should only be used to study a particular collisional event in detail. **helio_gr_test** @@ -522,11 +528,11 @@ Semi-interacting bodies are useful because the integrator is not required to cal **What should minimum fragment mass should I use (**```MIN_GMFRAG``` **or** ```MIN_MFRAG```**)?** -This mass threshold is necessary to ensure that Swiftest SyMBA does not generate huge amounts of very small fragments, grinding the model to a halt. While this value is largely empirical and dependent on each specific set of initial conditions, a good place to start is to set the minimum fragment mass threshold to be one tenth the size of the smallest body in your simulation. +This mass threshold is necessary to ensure that Swiftest SyMBA does not generate huge amounts of very small fragments, grinding the model to a halt. While this value is largely empirical and dependent on each specific set of initial conditions, a good place to start is to set the minimum fragment mass threshold to be one tenth the size of the smallest body in your simulation. You can also adjust ```FRAG_REDUCTION``` to keep the number of fragments within a reasonable range. **What are the limits of Swiftest SyMBA?** -While Swifest SyMBA is a powerful tool for modeling gravitational interactions between massive bodies, it does have its limits. While Swiftest SyMBA is capable of modeling systems containing thousands of massive bodies, the code does slow down significantly. For this reason, Swiftest SyMBA is best used for systems containing tens to hundreds of fully-interacting massive bodies. It is also best used for timescales on the order of a few hundred million years or less. While it is possible to model systems on a billion year timescale, the computational power required may be beyond what is available to the average user. In these cases, it is recommended that the user consider modeling with test particles instead of massive bodies. For systems that contain mainly test particles, with few to no close encounters between massive bodies, Swiftest RMVS is likely a more appropriate tool. +While Swifest SyMBA is a powerful tool for modeling gravitational interactions between massive bodies, it does have its limits. Swiftest SyMBA is best used for systems containing tens to hundreds of fully-interacting massive bodies. It is also best used for timescales on the order of a few hundred million years or less. While it is possible to model systems on a billion year timescale, the computational power required may be beyond what is available to the average user. In these cases, it is recommended that the user consider modeling with test particles instead of massive bodies. For systems that contain mainly test particles, with few to no close encounters between massive bodies, Swiftest RMVS is likely a more appropriate tool. To get a sense of the scope of your desired simulation, it is recommended that you run your initial conditions and parameters for a just few steps. Make sure that you set ```ISTEP_OUT``` and ```DUMP_CADENCE``` to output only once the simulation is complete, not between steps. Because writing to the output files and memory takes a significant amount of computational time compared to integrating the step, we want to avoid counting writing time in our diagnostic information. The terminal output contains information about the total wall time and the wall time per integration step. To get a sense of how long your run will take to complete your desired ```tmax```, simply scale up the wall time per integration step to the number of steps necessary for ```tmax``` to be reached. Remember that writing to the output files will take a considerable amount of time. Adjust your intitial conditions and parameters accordingly. @@ -564,4 +570,4 @@ Swiftest is free software: you can redistribute it and/or modify it under the te Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -You should have received a copy of the GNU General Public License along with Swiftest. If not, see . \ No newline at end of file +You should have received a copy of the GNU General Public License along with Swiftest. If not, see . diff --git a/README_tables/param_options.md b/README_tables/param_options.md deleted file mode 100644 index 6125a282b..000000000 --- a/README_tables/param_options.md +++ /dev/null @@ -1,57 +0,0 @@ -# **param.in** options -| Parameter Name | Parameter Description | Input Format | Compatible Integrators | -|-------------------------|------------------------------------------------------------------------------------------------------------------------------|---------------------------------------------------------------------------------------------------------|------------------------| -| ```T0``` | The reference time for the start of the simulation in time units. | floating point (ex. ```0.0```) | all -| ```TSTART``` | Simulation start time for a restarted run in time units. | floating point (ex. ```0.0```) | all -| ```TSTOP``` | Simulation end time in time units. Must be greater than ```TSTART```. | floating point (ex. ```100.0```) | all -| ```DT``` | Simulation step size in time units. Must be less than or equal to ```TSTOP```-```TSTART```. | floating point (ex. ```0.005```) | all -| ```ISTEP_OUT``` | The number of time steps (```DT```) between output saves to memory. | integer (ex. ```200```) | all -| ```DUMP_CADENCE``` | The number of output steps between when data saved to memory is written to file. Setting to ```0``` results in writing data to file only at the completion of the simulation. | integer (ex. ```100```) | all -| ```IN_TYPE``` | Input file format. | ```ASCII```, ```NETCDF_FLOAT```, and ```NETCDF_DOUBLE``` | all -| ```NC_IN``` | NetCDF input file name. Only if ```IN_TYPE``` is set to ```NETCDF_FLOAT``` or ```NETCDF_DOUBLE```. | string (ex. ```init_cond.nc```) | all -| ```PL_IN``` | Massive body input file name. Only if ```IN_TYPE``` is set to ```ASCII```. | string (ex. ```pl.in```) | all -| ```TP_IN``` | Test particle input file name. Only if ```IN_TYPE``` is set to ```ASCII```. | string (ex. ```tp.in```) | all -| ```CB_IN ``` | Central body input file name. Only if ```IN_TYPE``` is set to ```ASCII```. | string (ex. ```cb.in```) | all -| ```IN_FORM``` | Input format. | ```EL```, ```XV``` | all -| ```OUT_TYPE``` | Output file format. | ```NETCDF_FLOAT```, ```NETCDF_DOUBLE``` | all -| ```BIN_OUT``` | Output file name. | string (ex. ```bin.nc```) | all -| ```OUT_FORM``` | Output format. | ```XV```, ```XVEL``` | all -| ```OUT_STAT``` | Output status. | ```NEW```, ```APPEND```, ```REPLACE```, ```UNKNOWN``` | all -| ```CHK_QMIN``` | Pericenter distance at which a test particle is too close to the pericenter of the system in distance units. | floating point, turn off using ```-1.0``` | all -| ```CHK_RMIN``` | Heliocentric distance at which a test particle is considered merged with the central body in distance units. | floating point, turn off using ```-1.0``` | all -| ```CHK_RMAX``` | Heliocentric distance at which a test particle is too distant from the central body in distance units. | floating point (ex. ```1000.0```) | all -| ```CHK_EJECT``` | Heliocentric distance at which an unbound test particle is too distant from the central body in distance units. | floating point (ex. ```1000.0```) | all -| ```CHK_QMIN_COORD``` | Coordinate frame used to check for minimum pericenter distance. | ```HELIO```, ```BARY``` | all -| ```CHK_QMIN_RANGE``` | Upper and lower bounds of the semimajor axis range used to check the pericenter distance. | two floating points, turn off using ```-1.0 -1.0``` | all -| ```EXTRA_FORCE``` | Additional user defined force routines provided. | ```YES```, ```NO``` | all -| ```CHK_CLOSE``` | Check for close encounters. Requires radius of massive bodies to be provided in initial conditions. | ```YES```, ```NO``` | all -| ```INTERACTION_LOOPS``` | Method for checking for interactions between bodies. | ```TRIANGULAR```, ```FLAT```, ```ADAPTIVE``` | all -| ```ENCOUNTER_CHECK``` | Method for checking for close encounters between bodies. | ```TRIANGULAR```, ```SORTSWEEP```, ```ADAPTIVE``` | all -| ```MU2KG``` | Mass units to kilogram conversion factor. | floating point (ex. ```1.988409870698051e+30```) | all -| ```TU2S``` | Time units to seconds conversion factor. | floating point (ex. ```31557600.0```) | all -| ```DU2M``` | Distance units to meters conversion factor. | floating point (ex. ```149597870700.0```) | all -| ```BIG_DISCARD``` | Include data for all fully-interacting bodies (above GMTINY) in each discard. Swifter only. | ```YES```, ```NO``` | all -| ```GR``` | General relativity. | ```YES```, ```NO``` | all -| ```RHILL_PRESENT``` | Hill Radius present in massive body input file. | ```YES```, ```NO``` | SyMBA -| ```ENERGY``` | Track and report the total energy, angular momentum, and mass of the system. | ```YES```, ```NO``` | SyMBA -| ```COLLISION_MODEL``` | Resolve collisions. | ```MERGE```, ```BOUNCE```, ```FRAGGLE``` | SyMBA -| ```ROTATION``` | Rotation of massive bodies. Requires rotation vectors, radii, and moments of inertia to be provided in initial conditions. | ```YES```, ```NO``` | SyMBA -| ```GMTINY``` | Mass cutoff between fully and semi-interacting massive bodies in gravitational mass units. | floating point (ex. ```4e-06```) | SyMBA -| ```MIN_GMFRAG``` | Minimum fragment mass in gravitational mass units. | floating point (ex. ```1e-09```) | SyMBA -| ```TIDES``` | Tidal dissipation model. | ```YES```, ```NO``` | *(under development)* -| ```YORP``` | YORP effect. | ```YES```, ```NO``` | *(under development)* -| ```YARKOVSKY``` | Yarkovsky effect. | ```YES```, ```NO``` | *(under development)* - -In the above list, the following are defined as: -- ```HELIO``` - Use the heliocentric coordinate frame for ```CHK_QMIN``` -- ```BARY``` - Use the barycentric coordinate frame for ```CHK_QMIN``` -- ```XV``` - Heliocentric position and velocity components for ```IN_FORM``` and/or ```OUT_FORM``` -- ```EL``` - Osculating orbital elements for ```IN_FORM``` and/or ```OUT_FORM``` -- ```XVEL``` - Heliocentric position and velocity components and osculating orbital elements for ```OUT_FORM``` -- ```NETCDF_FLOAT``` - Single precision NetCDF format for ```OUT_TYPE``` -- ```NETCDF_DOUBLE``` - Double precision NetCDF format for ```OUT_TYPE``` -- ```MERGE``` - Perfectly conserve the mass of all colliding bodies into a single resultant body -- ```BOUNCE``` - Perfectly-elastic collision in which all bodies reverse trajectory after the collision -- ```FRAGGLE``` - Collisional fragments are generated as a result of a collision - -For more details on the ```INTERACTION_LOOPS``` and ```ENCOUNTER_CHECK``` options, see the **Updates to Swifter SyMBA** section below. \ No newline at end of file diff --git a/README_tables/simulation_kwargs.md b/README_tables/simulation_kwargs.md index 2f0b3f256..b2a7db5dd 100644 --- a/README_tables/simulation_kwargs.md +++ b/README_tables/simulation_kwargs.md @@ -1,53 +1,65 @@ # swiftest.Simulation(**kwargs) +The following key word arguments are specific to the ```swiftest.Simulation``` method. | Key Word Name | Key Word Description | Options | Compatible Integrators | |---------------------------------|----------------------------------------------------------------------------------------------------|----------------------------------------------------------------------------------------------------|------------------------| -|```simdir``` | Path to subdirectory in which to store data. Default is ```/simdir```. | pathlike string (ex. ```path/to/directory```) | all +|```simdir``` | Path to subdirectory in which to store data. Default is ```/simdir```. | pathlike string (ex. ```path/to/directory```) | all |```read_param``` | Read in a pre-existing parameter input file. Default is ```False```. | ```True```, ```False``` | all |```param_file``` | Name of the pre-existing parameter input file. Only used if ```read_param``` is set to ```True```. | string (ex. ```param.in```) | all -|```read_data``` | Read in a pre-existing simulation output file. Default is ```False```. | ```True```, ```False``` | all +|```read_data``` | Read in a pre-existing simulation output file. Default is ```False```. | ```True```, ```False``` | all |```codename``` | Name of the N-body code to use. Default is ```Swiftest```. | ```Swiftest```, ```Swifter```, ```Swift``` | all |```integrator``` | Name of the N-body integrator to use. Default is ```symba```. | ```symba```, ```helio```, ```rmvs```, ```whm``` | all -|```t0``` | The reference time for the start of the simulation in time units. Default is ```0.0```. | floating point (ex. ```0.0```) | all -|```tstart``` | Simulation start time for a restarted run in time units. Default is ```0.0```. | floating point (ex. ```0.0```) | all -|```tstop``` | Simulation end time in time units. Must be greater than ```tstart```. | floating point (ex. ```100.0```) | all -|```dt``` | Simulation step size in time units. Must be less than or equal to ```tstop```-```tstart```. | floating point (ex. ```0.005```) | all -|```istep_out``` | The number of time steps (```dt```) between output saves to memory. Either ```istep_out``` **OR** ```tstep_out``` may be set. Default is ```1```. | integer (ex. ```200```) | all -|```dump_cadence``` | The number of output steps between when data saved to memory is written to file. Setting to ```0``` results in writing data to file only at the completion of the simulation. Default is ```10```. | integer (ex. ```10```) | all -|```tstep_out``` | The approximate time between when outputs saved in memory are written to file in time units. Either ```istep_out``` **OR** ```tstep_out``` may be set. | floating point (ex. ```10.0```) that calculates ```istep_out = floor(tstep_out / dt)``` | all -|```init_cond_file_type``` | Input file format. Default is ```NETCDF_DOUBLE```. | ```NETCDF_DOUBLE```, ```NETCDF_FLOAT```, ```ASCII``` | all -|```init_cond_file_name``` | Input file name(s). If ```init_cond_file_type``` is set to ```NETCDF_DOUBLE``` or ```NETCDF_FLOAT```, default is ```init_cond.nc```. If ```init_cond_file_type``` is set to ```ASCII```, default is a dictionary: ```{"CB" : "cb.in", "PL" : "pl.in", "TP" : "tp.in"}```. | string (ex. ```my_init_cond.nc```) or dictionary ```{"CB" : "mycb.in", "PL" : "mypl.in", "TP" : "mytp.in"}``` | all -|```init_cond_format``` | Input format. Default is ```EL```. | ```EL```, ```XV``` | all -|```output_file_type``` | Output file format. Default is ```NETCDF_DOUBLE```. ```REAL4```, ```REAL8```, ```XDR4```, ```XDR8``` Swifter/Swift only. | ```NETCDF_DOUBLE```, ```NETCDF_FLOAT```, ```REAL4```, ```REAL8```, ```XDR4```, ```XDR8``` | all -|```output_file_name``` | Output file name. Default is ```bin.nc```. | string (ex. ```mydata.nc```) | all -|```output_format``` | Output format. Default is ```XVEL```. | ```XV```, ```XVEL``` | all -|```MU``` | Mass unit system to use in the simulation. Default is ```Msun```. | ```Msun```, ```Mearth```, ```kg```, ```g``` (case-insensitive) | all -|```DU``` | Distance unit system to use in the simulation. Default is ```AU```. | ```AU```, ```Rearth```, ```m```, ```cm``` (case-insensitive) | all -|```TU``` | Time unit system to use in the simulation. Default is ```Y```. | ```Y```, ```YR```, ```DAY``` (Julian day), ```d``` (Julian day), ```JD``` (Julian day), ```s``` (case-insensitive) | all -|```MU_name``` | The name of the mass unit. Overrides ```MU```. | string (ex. ```kilograms```) | all -|```DU_name``` | The name of the distance unit. Overrides ```DU```. | string (ex. ```meters```) | all -|```TU_name``` | The name of the time unit. Overrides ```TU```. | string (ex. ```seconds```) | all -|```MU2KG``` | Mass units to kilogram conversion factor. Overrides ```MU```. | floating point (ex. ```1.988409870698051e+30```) | all -|```DU2M``` | Distance units to meters conversion factor. Overrides ```DU```. | floating point (ex. ```31557600.0```) | all -|```TU2S``` | Time units to seconds conversion factor. Overrides ```TU```. | floating point (ex. ```149597870700.0```) | all -|```rmin``` | Heliocentric distance at which a test particle is considered merged with the central body in distance units. Default is the radius of the central body in system units. | floating point (ex. ```0.3```) | all -|```rmax``` | Heliocentric distance at which a test particle is too distant from the central body in distance units. Default is ```10000.0 AU```. | floating point (ex. ```10000.0```) | all -|```qmin_coord``` | Coordinate frame used to check for minimum pericenter distance. Default is ```HELIO```. | ```HELIO```, ```BARY``` | all -|```mtiny``` | Mass cutoff between fully and semi-interacting massive bodies in mass units. Either ```mtiny``` **OR** ```gmtiny``` may be set. | floating point (ex. ```1e23```) | all -|```gmtiny``` | Mass cutoff between fully and semi-interacting massive bodies in gravitational mass units. Default is ```0.0```. Either ```mtiny``` **OR** ```gmtiny``` may be set. | floating point (ex. ```4e-6```) | all -|```close_encounter_check``` | Check for close encounters. Default is ```True```. Requires radius of massive bodies to be provided in initial conditions. | ```True```, ```False``` | all -|```encounter_save``` | Save data for each close encounter to a file. Warning! This can generate very large files! | ```TRAJECTORY```, ```CLOSEST```, ```BOTH``` | SyMBA -|```general_relativity``` | General relativity. Default is ```True```. | ```True```, ```False``` | all -|```collision_model``` | Resolve collisions. Default is ```MERGE```. | ```MERGE```, ```BOUNCE```, ```FRAGGLE``` | SyMBA -|```minimum_fragment_gmass``` | Minimum fragment mass in gravitational mass units. Default is ```0.0```. Either ```minimum_fragment_gmass``` **OR** ```minimum_fragment_mass``` may be set. | floating point (ex. ```1e-9```) | SyMBA -|```minimum_fragment_mass``` | Minimum fragment mass in mass units. Either ```minimum_fragment_gmass``` **OR** ```minimum_fragment_mass``` may be set. | floating point (ex. ```1e20```) | SyMBA -|```rotation``` | Rotation of massive bodies. Requires rotation vectors, radii, and moments of inertia to be provided in initial conditions. Default is ```False```. | ```True```, ```False``` | SyMBA -|```compute_conservation_values```| Track and report the total energy, angular momentum, and mass of the system. Default is ```False```. | ```True```, ```False``` | SyMBA -|```rhill_present``` | Hill Radius present in massive body input file. Default is ```False```. | ```True```, ```False``` | SyMBA -|```extra_force``` | Additional user defined force routines provided. Default is ```False```. | ```True```, ```False``` | all -|```big_discard``` | Include data for all fully-interacting bodies (above GMTINY) in each discard. Swifter only. Default is ```False```. | ```True```, ```False``` | all -|```restart``` | If ```True```, the simulation given by ```output_file_name``` will be restarted from ```t0```. Default is ```False```. | ```True```, ```False``` | all -|```interaction_loops``` | Method for checking for interactions between bodies. Default is ```TRIANGULAR```. | ```TRIANGULAR```, ```FLAT```, ```ADAPTIVE``` | all -|```encounter_check_loops``` | Method for checking for close encounters between bodies. Default is ```TRIANGULAR```. | ```TRIANGULAR```, ```SORTSWEEP```, ```ADAPTIVE``` | all + +# Parameters +These parameters can be entered as key word arguments into the ```swiftest.Simulation``` method, the ```swiftest.get_parameter``` method, the ```swiftest.set_parameter``` method, the ```swiftest.run``` method, or through editing the ASCII **param.in** directly, as described in the ```Usage``` column. If the key word argument or parameter's name changes between usage as a Python key word argument and as a parameter in the ASCII input files, the name of the key word argument is listed first, followed by the name of the ASCII parameter in the ```Key Word / Parameter Name``` column. Many parameters retain their Swift / Swifter names in the ASCII for backwards compatibility, but change names in the Python for readability. + +| Key Word / Parameter Name | Key Word / Parameter Description | Options | Compatible Integrators | Usage | +|---------------------------------|----------------------------------------------------------------------------------------------------|----------------------------------------------------------------------------------------------------|------------------------|-------| +|```t0``` | The reference time for the start of the simulation in time units. Default is ```0.0```. | floating point (ex. ```0.0```) | all | Both | +|```tstart``` | Simulation start time for a restarted run in time units. Default is ```0.0```. | floating point (ex. ```0.0```) | all | Both | +|```tstop``` | Simulation end time in time units. Must be greater than ```tstart```. | floating point (ex. ```100.0```) | all | Both | +|```dt``` | Simulation step size in time units. Must be less than or equal to ```tstop```-```tstart```. | floating point (ex. ```0.005```) | all | Both | +|```istep_out``` | The number of time steps (```dt```) between output saves to memory. Either ```istep_out``` **OR** ```tstep_out``` may be set. Default is ```1```. | integer (ex. ```200```) | all | Both | +|```tstep_out``` | The approximate time between when outputs saved in memory are written to file in time units. Either ```istep_out``` **OR** ```tstep_out``` may be set. Calculates ```istep_out = floor(tstep_out / dt)``` | floating point (ex. ```10.0```) | all | Both +|```dump_cadence``` | The number of output steps between when data saved to memory is written to file. Setting to ```0``` results in writing data to file only at the completion of the simulation. Default is ```10```. | integer (ex. ```10```) | all | Both | +|```init_cond_file_type``` / ```IN_TYPE``` | Input file type. Default is ```NETCDF_DOUBLE```. | ```NETCDF_DOUBLE```, ```NETCDF_FLOAT```, ```ASCII``` | all | Both | +|```init_cond_file_name``` / ```PL_IN```, ```TP_IN```, ```CB_IN```, ```NC_IN``` | Input file name(s). If using Python and ```init_cond_file_type``` is set to ```NETCDF_DOUBLE``` or ```NETCDF_FLOAT```, default is ```init_cond.nc```. If using Python and ```init_cond_file_type``` is set to ```ASCII```, default is a dictionary: ```{"CB" : "cb.in", "PL" : "pl.in", "TP" : "tp.in"}```. If editing the ASCII **param.in**, each input file name is a separate parameter. If ```IN_TYPE``` is set to ```NETCDF_DOUBLE``` or ```NETCDF_FLOAT```, default is ```init_cond.nc```. If ```IN_TYPE``` is set to ```ASCII```, defaults are ```CB.IN```, ```PL.IN```, and ```TP.in```. | string (ex. ```my_init_cond.nc```), dictionary ```{"CB" : "mycb.in", "PL" : "mypl.in", "TP" : "mytp.in"}``` | all | Both +|```init_cond_format``` / ```IN_FORM``` | Input format. Default is ```EL```. | ```EL```, ```XV``` | all | Both | +|```output_file_type``` / ```OUT_TYPE``` | Output file type. Default is ```NETCDF_DOUBLE```. | ```NETCDF_DOUBLE```, ```NETCDF_FLOAT```, Swift / Swifter Only : ```REAL4```, ```REAL8```, ```XDR4```, ```XDR8``` | all | Both | +|```output_file_name``` / ```BIN_OUT``` | Output file name. Default is ```bin.nc```. | string (ex. ```mydata.nc```) | all | Both | +|```output_format``` / ```OUT_FORM``` | Output format. Default is ```XVEL```. | ```XV```, ```XVEL``` | all | Both | +| ```OUT_STAT``` | Output status. Default is ```REPLACE```. | ```NEW```, ```APPEND```, ```REPLACE```, ```UNKNOWN``` | all | ASCII only | +|```MU``` | Mass unit system to use in the simulation. Default is ```Msun```. | ```Msun```, ```Mearth```, ```kg```, ```g``` (case-insensitive) | all | Python only | +|```DU``` | Distance unit system to use in the simulation. Default is ```AU```. | ```AU```, ```Rearth```, ```m```, ```cm``` (case-insensitive) | all | Python only | +|```TU``` | Time unit system to use in the simulation. Default is ```Y```. | ```Y```, ```YR```, ```DAY``` (Julian day), ```d``` (Julian day), ```JD``` (Julian day), ```s``` (case-insensitive) | all | Python only | +|```MU_name``` | The name of the mass unit. Overrides ```MU```. | string (ex. ```kilograms```) | all | Python only | +|```DU_name``` | The name of the distance unit. Overrides ```DU```. | string (ex. ```meters```) | all | Python only | +|```TU_name``` | The name of the time unit. Overrides ```TU```. | string (ex. ```seconds```) | all | Python only | +|```MU2KG``` | Mass units to kilogram conversion factor. Overrides ```MU```. | floating point (ex. ```1.988409870698051e+30```) | all | Both | +|```DU2M``` | Distance units to meters conversion factor. Overrides ```DU```. | floating point (ex. ```31557600.0```) | all | Both | +|```TU2S``` | Time units to seconds conversion factor. Overrides ```TU```. | floating point (ex. ```149597870700.0```) | all | Both | +|```rmin``` / ```CHK_RMIN``` | Heliocentric distance at which a test particle is considered merged with the central body in distance units. Default is the radius of the central body in system units. | floating point (ex. ```0.3```) | all | Both | +|```rmax``` / ```CHK_RMAX``` | Heliocentric distance at which a test particle is too distant from the central body in distance units. Default is ```10000.0 AU```. | floating point (ex. ```10000.0```) | all | Both | +|```qmin_coord``` / ```CHK_QMIN_COORD``` | Coordinate frame used to check for minimum pericenter distance. Default is ```HELIO```. | ```HELIO```, ```BARY``` | all | Both | +| ```CHK_QMIN``` | Pericenter distance at which a test particle is too close to the pericenter of the system in distance units. | floating point, turn off using ```-1.0``` | all | ASCII only | +| ```CHK_QMIN_RANGE``` | Upper and lower bounds of the semimajor axis range used to check the pericenter distance. | two floating points, turn off using ```-1.0 -1.0``` | all | ASCII only | +| ```CHK_EJECT``` | Heliocentric distance at which an unbound test particle is too distant from the central body in distance units. | floating point (ex. ```1000.0```) | all | ASCII only | +|```extra_force``` / ```EXTRA_FORCE``` | Additional user defined force routines provided. Default is ```False``` / ```NO```. | ```True```, ```False``` / ```YES```, ```NO``` | all | Both | +|```close_encounter_check``` / ```CHK_CLOSE``` | Check for close encounters. Default is ```True``` / ```YES```. Requires radius of massive bodies to be provided in initial conditions. | ```True```, ```False``` / ```YES```, ```NO``` | all | Both | +|```interaction_loops``` | Method for checking for interactions between bodies. Default is ```TRIANGULAR```. | ```TRIANGULAR```, ```FLAT```, ```ADAPTIVE``` | all | Both | +|```encounter_check_loops``` / ```ENCOUNTER_CHECK``` | Method for checking for close encounters between bodies. Default is ```TRIANGULAR```. | ```TRIANGULAR```, ```SORTSWEEP```, ```ADAPTIVE``` | all | Both | +|```encounter_save``` | Save data for each close encounter to a file. Warning! This can generate very large files! | ```TRAJECTORY```, ```CLOSEST```, ```BOTH``` | SyMBA | Python only | +|```collision_model``` | Resolve collisions. Default is ```MERGE```. | ```MERGE```, ```BOUNCE```, ```FRAGGLE``` | SyMBA | Both | +|```nfrag_reduction``` | Factor to reduce the number of fragments generated by a collision. See below for a more detailed explanation. | floating point, default is ```30.0```, minumum is ```1.0``` | SyMBA | Both | +|```minimum_fragment_gmass``` / ```MIN_GMFRAG``` | Minimum fragment mass in gravitational mass units. Default is ```0.0```. Either ```minimum_fragment_gmass``` **OR** ```minimum_fragment_mass``` may be set. | floating point (ex. ```1e-9```) | SyMBA | Both | +|```minimum_fragment_mass``` | Minimum fragment mass in mass units. Either ```minimum_fragment_gmass``` **OR** ```minimum_fragment_mass``` may be set. | floating point (ex. ```1e20```) | SyMBA | Python only | +|```rotation``` | Rotation of massive bodies. Requires rotation vectors, radii, and moments of inertia to be provided in initial conditions. Default is ```False```. | ```True```, ```False``` / ```YES```, ```NO``` | SyMBA | Both | +|```rhill_present``` | Hill Radius present in massive body input file. Default is ```False```. | ```True```, ```False``` / ```YES```, ```NO``` | SyMBA | Both | +|```compute_conservation_values``` / ```ENERGY``` | Track and report the total energy, angular momentum, and mass of the system. Default is ```False```. | ```True```, ```False``` / ```YES```, ```NO``` | SyMBA | Both | +|```general_relativity``` / ```GR``` | General relativity. Default is ```True```. | ```True```, ```False``` / ```YES```, ```NO``` | all | Both | +|```big_discard``` | Include data for all fully-interacting bodies (above GMTINY) in each discard. Swifter only. Default is ```False```. | ```True```, ```False``` / ```YES```, ```NO``` | all | Both | +|```restart``` | If ```True```, the simulation given by ```output_file_name``` will be restarted from ```t0```. Default is ```False```. | ```True```, ```False``` | all | Python only | +|```gmtiny``` | Mass cutoff between fully and semi-interacting massive bodies in gravitational mass units. Default is ```0.0```. Either ```mtiny``` **OR** ```gmtiny``` may be set. | floating point (ex. ```4e-6```) | all | Both | +|```mtiny``` | Mass cutoff between fully and semi-interacting massive bodies in mass units. Either ```mtiny``` **OR** ```gmtiny``` may be set. | floating point (ex. ```1e23```) | all | Python only | In the above list, the following are defined as: - ```HELIO``` - Use the heliocentric coordinate frame. @@ -65,4 +77,7 @@ In the above list, the following are defined as: - ```CLOSEST``` - Save the position and velocity vectors of all bodies involved in the close encounter at the point of closest approach - ```MERGE``` - Perfectly conserve the mass of all colliding bodies into a single resultant body - ```BOUNCE``` - Perfectly-elastic collision in which all bodies reverse trajectory after the collision -- ```FRAGGLE``` - Collisional fragments are generated as a result of a collision \ No newline at end of file +- ```FRAGGLE``` - Collisional fragments are generated as a result of a collision +- ```NFRAG_REDUCTION``` - Based on the SFD defined in Leinhardt & Stewart (2012), Fraggle dictates that ```N``` fragments be generated between the mass of the second largest remnant and ```MIN_GMFRAG```, with the maximum value of ```N``` being ```3000```. To improve performance, this number can be reduced by ```N / NFRAG_REDUCTION```. + +For more details on the ```INTERACTION_LOOPS``` and ```ENCOUNTER_CHECK``` options, see the **Updates to Swifter SyMBA** section. \ No newline at end of file diff --git a/examples/Chambers2013/init_cond.py b/examples/Chambers2013/init_cond.py index 1b2437b6c..d5c9f9f19 100755 --- a/examples/Chambers2013/init_cond.py +++ b/examples/Chambers2013/init_cond.py @@ -43,7 +43,7 @@ mtiny = 1e-2 * Ms minimum_fragment_mass = 1e-5 * Ms -nfrag_reduction = 10.0 +nfrag_reduction = 30.0 rng = default_rng(seed=3031179) runname = "Chambers (2013)" diff --git a/examples/Chambers2013/scattermovie.py b/examples/Chambers2013/scattermovie.py index a5e99de4d..98b0645ef 100755 --- a/examples/Chambers2013/scattermovie.py +++ b/examples/Chambers2013/scattermovie.py @@ -152,7 +152,7 @@ def update(self,frame): return self.artists -sim = swiftest.Simulation(read_data=True) +sim = swiftest.Simulation(read_data=True, read_collisions=False, dask=True) for plot_style in valid_plot_styles: animation_file = f"Chambers2013-{plot_style}.mp4" print('Making animation') diff --git a/examples/Fragmentation/.gitignore b/examples/Fragmentation/.gitignore index 2828073f0..fadcdbd6c 100644 --- a/examples/Fragmentation/.gitignore +++ b/examples/Fragmentation/.gitignore @@ -1,5 +1,4 @@ * !.gitignore !Fragmentation_Movie.py -!swiftest_fragmentation.py !README.txt diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 527c02a2d..c44fe0b56 100755 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -11,16 +11,26 @@ """ """ -Generates a movie of a fragmentation event from set of Swiftest output files. +Generates and runs a set of Swiftest input files from initial conditions with the SyMBA integrator. All simulation +outputs are stored in the subdirectory named after their collisional regime. Inputs _______ -param.in : ASCII Swiftest parameter input file. -data.nc : A NetCDF file containing the simulation output. +None. + +Output +------ +collisions.log : An ASCII file containing the information of any collisional events that occured. +collisions.nc : A NetCDF file containing the collision output. +data.nc : A NetCDF file containing the simulation output. +encounters.nc : A NetCDF file containing the encounter output. +init_cond.nc : A NetCDF file containing the initial conditions for the simulation. +param.00...0.in : A series of parameter input files containing the parameters for the simulation at every output stage. +param.in : An ASCII file containing the inital parameters for the simulation. +param.restart.in : An ASCII file containing the parameters for the simulation at the last output. +swiftest.log : An ASCII file containing the information on the status of the simulation as it runs. +collision.mp4 : A movie file named after the collisional regime depicting the collision. -Returns -------- -fragmentation.mp4 : A .mp4 file of a fragmentation event. """ import swiftest @@ -156,7 +166,7 @@ def encounter_combiner(sim): enc = enc.sel(time=tgood) # The following will combine the two datasets along the time dimension, sort the time dimension, and then fill in any time gaps with interpolation - ds = xr.combine_nested([data,enc],concat_dim='time').sortby("time").interpolate_na(dim="time", method="akima") + ds = xr.combine_nested([data,enc],concat_dim='time').sortby("time").interpolate_na(dim="time") # Rename the merged Target body so that their data can be combined tname=[n for n in ds['name'].data if names[0] in n] diff --git a/examples/Fragmentation/README.txt b/examples/Fragmentation/README.txt index c37f52296..103e01b18 100644 --- a/examples/Fragmentation/README.txt +++ b/examples/Fragmentation/README.txt @@ -16,7 +16,6 @@ Date : December 6, 2022 Included in the Fragmentation example directory are the following files: - README.txt : This file - - swiftest_fragmentation.py : A Python Script that generates and runs three sets of initial conditions. - - Fragmentation_Movie.py : A Python Script that processes a data.nc file and generates a movie (.mp4) of a collisional event. + - Fragmentation_Movie.py : A Python Script that generates and runs initial conditions of collisional events, and then generates a movie (.mp4) of the collisional event. This example is intended to be run with Swiftest SyMBA. For details on how to generate, run, and analyze this example, see the Swiftest User Manual. \ No newline at end of file diff --git a/examples/Fragmentation/swiftest_fragmentation.py b/examples/Fragmentation/swiftest_fragmentation.py deleted file mode 100644 index 4c2368111..000000000 --- a/examples/Fragmentation/swiftest_fragmentation.py +++ /dev/null @@ -1,79 +0,0 @@ -#!/usr/bin/env python3 - -""" - Copyright 2023 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh - This file is part of Swiftest. - Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License - as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. - Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty - of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. - You should have received a copy of the GNU General Public License along with Swiftest. - If not, see: https://www.gnu.org/licenses. -""" - -""" -Generates and runs a set of Swiftest input files from initial conditions with the SyMBA integrator. All simulation -outputs for the disruption case are stored in the /disruption subdirectory. All simulation outputs for the hit and run -case are stored in the /hitandrun subdirectory. All simulation outputs for the super-catastrophic disruption case are -stored in the /supercat subdirectory. - -Input ------- -None. - -Output ------- -Three subdirectories: -disruption/ -hitandrun/ -supercat/ - -Each subdirectory contains: -data.nc : A NetCDF file containing the simulation output. -init_cond.nc : A NetCDF file containing the initial conditions for the simulation. -collision_000001.nc : A NetCDF file containing the data for the collisions. -encounter_000001.nc : A NetCDF file containing the data for the close encounters. -collisions.log : An ASCII file containing the information of any collisional events that occured. -param.in : An ASCII file containing the parameters for the simulation. -swiftest.log : An ASCII file containing the information on the status of the simulation as it runs. - -""" -import swiftest -import numpy as np -from numpy.random import default_rng - -run_arguments = {"tstart": 0.0, "tstop":1e-5, "dt": 1e-5, "istep_out": 1, "fragmentation":True, "encounter_save":"both", "compute_conservation_values":True, - "minimum_fragment_gmass":1.0e-11, "gmtiny":1.0e-11, "output_format":"XVEL", "init_cond_format":"XV"} - -# Initialize the simulation object as a variable with arguments. -sim_disruption = swiftest.Simulation(simdir="disruption", **run_arguments) -# Add the Sun using the JPL Horizons Database. -sim_disruption.add_solar_system_body(["Sun"]) -# Add a user-defined target body. -sim_disruption.add_body(name="Target", rh=[1.0, -1.807993e-05, 0.0], vh=[-2.562596e-04, 6.280005, 0.0], Gmass=1e-7, radius=7e-6) -# Add a user-defined projectile body. -sim_disruption.add_body(name="Projectile", rh=[1.0, 1.807993e-05, 0.0], vh=[-2.562596e-04, -6.280005, 0.0], Gmass=7e-10, radius=3.25e-6) -# Display the run configuration parameters. -sim_disruption.get_parameter() -# Write the parameters to the param.in -sim_disruption.write_param() -# Run the simulation. -sim_disruption.run() - -# Do the same as above for the hit and run case. -sim_hitandrun = swiftest.Simulation(simdir="hitandrun", **run_arguments) -sim_hitandrun.add_solar_system_body(["Sun"]) -sim_hitandrun.add_body(name="Target", rh=[1.0, -4.2e-05, 0.0], vh=[0.0, 6.28, 0.0], Gmass=1e-7, radius=7e-6, rot=[0.0, 0.0, 6.0e4]) -sim_hitandrun.add_body(name="Projectile", rh=[1.0, 4.2e-05, 0.0], vh=[-1.5, -6.28, 0.0], Gmass=7e-10, radius=3.25e-6, rot=[0.0, 0.0, 1.0e5]) -sim_hitandrun.get_parameter() -sim_hitandrun.write_param() -sim_hitandrun.run() - -# Do the same as above for the super-catastrophic disruption case. -sim_supercat = swiftest.Simulation(simdir="supercat", **run_arguments) -sim_supercat.add_solar_system_body(["Sun"]) -sim_supercat.add_body(name="Target", rh=[1.0, -4.2e-05, 0.0], vh=[0.0, 6.28, 0.0], Gmass=1e-7, radius=7e-6, rot=[0.0, 0.0, -6.0e4]) -sim_supercat.add_body(name="Projectile", rh=[1.0, 4.2e-05, 0.0], vh=[1.0, -6.28, 0.0], Gmass=1e-8, radius=3.25e-6, rot=[0.0, 0.0, 1.0e5]) -sim_supercat.get_parameter() -sim_supercat.write_param() -sim_supercat.run() diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 98e6c154a..6a376d387 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -825,7 +825,7 @@ def process_netcdf_input(ds, param): return ds -def swiftest2xr(param, verbose=True): +def swiftest2xr(param, verbose=True, dask=False): """ Converts a Swiftest binary data file into an xarray DataSet. @@ -833,6 +833,8 @@ def swiftest2xr(param, verbose=True): ---------- param : dict Swiftest parameters + dask : bool, default False + Use Dask to lazily load data (useful for very large datasets) Returns ------- @@ -842,7 +844,11 @@ def swiftest2xr(param, verbose=True): if ((param['OUT_TYPE'] == 'NETCDF_DOUBLE') or (param['OUT_TYPE'] == 'NETCDF_FLOAT')): if verbose: print('\nCreating Dataset from NetCDF file') - ds = xr.open_dataset(param['BIN_OUT'], mask_and_scale=False) + if dask: + ds = xr.open_mfdataset(param['BIN_OUT'], engine='h5netcdf', mask_and_scale=False) + else: + ds = xr.open_dataset(param['BIN_OUT'], mask_and_scale=False) + ds = process_netcdf_input(ds, param) else: print(f"Error encountered. OUT_TYPE {param['OUT_TYPE']} not recognized.") diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 1e851b0fe..10ebba59e 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -49,6 +49,7 @@ def __init__(self,read_param: bool = False, read_collisions: bool | None = None, read_encounters: bool | None = None, simdir: os.PathLike | str = "simdata", + dask: bool = False, **kwargs: Any): """ @@ -300,6 +301,8 @@ def __init__(self,read_param: bool = False, time calculation. The exact floating-point results of the interaction will be different between the two algorithm types. Parameter input file equivalent: `ENCOUNTER_CHECK` + dask : bool, default False + Use Dask to lazily load data (useful for very large datasets) verbose : bool, default True If set to True, then more information is printed by Simulation methods as they are executed. Setting to False suppresses most messages other than errors. @@ -367,7 +370,7 @@ def __init__(self,read_param: bool = False, self.read_encounters = read_encounters if read_param or read_data: - if self.read_param(read_init_cond = True): + if self.read_param(read_init_cond = True, dask=dask): # We will add the parameter file to the kwarg list. This will keep the set_parameter method from # overriding everything with defaults when there are no arguments passed to Simulation() kwargs['param_file'] = self.param_file @@ -390,7 +393,7 @@ def __init__(self,read_param: bool = False, if read_data: binpath = os.path.join(self.simdir, self.param['BIN_OUT']) if os.path.exists(binpath): - self.read_output_file() + self.read_output_file(dask=dask) else: raise FileNotFoundError(f"BIN_OUT file {binpath} not found.") return @@ -472,7 +475,7 @@ def _type_scrub(output_data): pbar.close() return - def run(self,**kwargs): + def run(self,dask: bool = False, **kwargs): """ Runs a Swiftest integration. Uses the parameters set by the `param` dictionary unless overridden by keyword arguments. Accepts any keyword arguments that can be passed to `set_parameter`. @@ -515,7 +518,7 @@ def run(self,**kwargs): # Read in new data self.read_encounters = True self.read_collisions = True - self.read_output_file() + self.read_output_file(dask=dask) return @@ -2592,7 +2595,8 @@ def read_param(self, param_file : os.PathLike | str | None = None, codename: Literal["Swiftest", "Swifter", "Swift"] | None = None, read_init_cond : bool | None = None, - verbose: bool | None = None): + verbose: bool | None = None, + dask: bool = False): """ Reads in an input parameter file and stores the values in the param dictionary. @@ -2607,6 +2611,9 @@ def read_param(self, verbose : bool, default is the value of the Simulation object's internal `verbose` If set to True, then more information is printed by Simulation methods as they are executed. Setting to False suppresses most messages other than errors. + dask : bool, default False + Use Dask to lazily load data (useful for very large datasets) + Returns ------- True if the parameter file exists and is read correctly. False otherwise. @@ -2632,7 +2639,7 @@ def read_param(self, if os.path.exists(init_cond_file): param_tmp = self.param.copy() param_tmp['BIN_OUT'] = init_cond_file - self.data = io.swiftest2xr(param_tmp, verbose=self.verbose) + self.data = io.swiftest2xr(param_tmp, verbose=self.verbose, dask=dask) self.init_cond = self.data.copy(deep=True) else: warnings.warn(f"Initial conditions file file {init_cond_file} not found.", stacklevel=2) @@ -2715,24 +2722,26 @@ def write_param(self, return def convert(self, param_file, newcodename="Swiftest", plname="pl.swiftest.in", tpname="tp.swiftest.in", - cbname="cb.swiftest.in", conversion_questions={}): + cbname="cb.swiftest.in", conversion_questions={}, dask = False): """ Converts simulation input files from one format to another (Swift, Swifter, or Swiftest). Parameters ---------- param_file : string - File name of the input parameter file + File name of the input parameter file newcodename : string - Name of the desired format (Swift/Swifter/Swiftest) + Name of the desired format (Swift/Swifter/Swiftest) plname : string - File name of the massive body input file + File name of the massive body input file tpname : string - File name of the test particle input file + File name of the test particle input file cbname : string - File name of the central body input file + File name of the central body input file conversion_questions : dictronary - Dictionary of additional parameters required to convert between formats + Dictionary of additional parameters required to convert between formats + dask : bool, default False + Use Dask to lazily load data (useful for very large datasets) Returns ------- @@ -2768,15 +2777,16 @@ def convert(self, param_file, newcodename="Swiftest", plname="pl.swiftest.in", t warnings.warn(f"Conversion from {self.codename} to {newcodename} is not supported.",stacklevel=2) return oldparam - def read_output_file(self,read_init_cond : bool = True): + def read_output_file(self,read_init_cond : bool = True, dask : bool = False): """ Reads in simulation data from an output file and stores it as an Xarray Dataset in the `data` instance variable. Parameters ---------- - read_init_cond : bool - Read in an initial conditions file along with the output file. Default is True - + read_init_cond : bool, default True + Read in an initial conditions file along with the output file. Default is True + dask : bool, default False + Use Dask to lazily load data (useful for very large datasets) Returns ------- self.data : xarray dataset @@ -2789,7 +2799,7 @@ def read_output_file(self,read_init_cond : bool = True): param_tmp = self.param.copy() param_tmp['BIN_OUT'] = os.path.join(self.simdir, self.param['BIN_OUT']) if self.codename == "Swiftest": - self.data = io.swiftest2xr(param_tmp, verbose=self.verbose) + self.data = io.swiftest2xr(param_tmp, verbose=self.verbose, dask=dask) if self.verbose: print('Swiftest simulation data stored as xarray DataSet .data') if read_init_cond: @@ -2797,14 +2807,14 @@ def read_output_file(self,read_init_cond : bool = True): print("Reading initial conditions file as .init_cond") if "NETCDF" in self.param['IN_TYPE']: param_tmp['BIN_OUT'] = self.simdir / self.param['NC_IN'] - self.init_cond = io.swiftest2xr(param_tmp, verbose=False) + self.init_cond = io.swiftest2xr(param_tmp, verbose=False, dask=dask) else: self.init_cond = self.data.isel(time=0) if self.read_encounters: - self.read_encounter_file() + self.read_encounter_file(dask=dask) if self.read_collisions: - self.read_collision_file() + self.read_collision_file(dask=dask) if self.verbose: print("Finished reading Swiftest dataset files.") @@ -2817,12 +2827,15 @@ def read_output_file(self,read_init_cond : bool = True): warnings.warn('Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".',stacklevel=2) return - def read_encounter_file(self): + def read_encounter_file(self, dask=False): enc_file = self.simdir / "encounters.nc" if not os.path.exists(enc_file): return - self.encounters = xr.open_dataset(enc_file) + if dask: + self.encounters = xr.open_mfdataset(enc_file,engine='h5netcdf', mask_and_scale=False) + else: + self.encounters = xr.open_dataset(enc_file, mask_and_scale=False) self.encounters = io.process_netcdf_input(self.encounters, self.param) # Remove any overlapping time values @@ -2834,7 +2847,7 @@ def read_encounter_file(self): return - def read_collision_file(self): + def read_collision_file(self, dask=False): col_file = self.simdir / "collisions.nc" if not os.path.exists(col_file): @@ -2843,12 +2856,16 @@ def read_collision_file(self): if self.verbose: print("Reading collisions history file as .collisions") - self.collisions = xr.open_dataset(col_file) + if dask: + self.collisions = xr.open_mfdataset(col_file,engine='h5netcdf', mask_and_scale=False) + else: + self.collisions = xr.open_dataset(col_file, mask_and_scale=False) + self.collisions = io.process_netcdf_input(self.collisions, self.param) return - def follow(self, codestyle="Swifter"): + def follow(self, codestyle="Swifter", dask=False): """ An implementation of the Swift tool_follow algorithm. Under development. Currently only for Swift simulations. @@ -2862,7 +2879,7 @@ def follow(self, codestyle="Swifter"): fol : xarray dataset """ if self.data is None: - self.read_output_file() + self.read_output_file(dask=dask) if codestyle == "Swift": try: with open('follow.in', 'r') as f: diff --git a/python/swiftest/swiftest/tool.py b/python/swiftest/swiftest/tool.py index 7fe2c8884..fdcd5585f 100644 --- a/python/swiftest/swiftest/tool.py +++ b/python/swiftest/swiftest/tool.py @@ -18,7 +18,7 @@ def magnitude(ds,x): dim = "space" ord = None return xr.apply_ufunc( - np.linalg.norm, ds[x], input_core_dims=[[dim]], kwargs={"ord": ord, "axis": -1} + np.linalg.norm, ds[x].where(~np.isnan(ds[x])), input_core_dims=[[dim]], kwargs={"ord": ord, "axis": -1}, dask="allowed" ) def wrap_angle(angle): diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 2957dffaf..042c67af7 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -67,7 +67,8 @@ module collision real(DP), dimension(NDIM,2) :: L_spin !! Two-body equivalent spin angular momentum vectors of the collider bodies prior to collision real(DP), dimension(NDIM,2) :: L_orbit !! Two-body equivalent orbital angular momentum vectors of the collider bodies prior to collision real(DP), dimension(2) :: ke_orbit !! Orbital kinetic energy of each individual impactor - real(DP), dimension(2) :: ke_spin !! Spin kinetic energy of each individual impactors + real(DP), dimension(2) :: ke_spin !! Spin kinetic energy of each individual impactor + real(DP), dimension(2) :: be !! Binding energy of each individual impactor real(DP), dimension(NDIM,2) :: Ip !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision real(DP), dimension(2) :: Gmass !! Two-body equivalent G*mass of the collider bodies prior to the collision real(DP), dimension(2) :: mass !! Two-body equivalent mass of the collider bodies prior to the collision @@ -404,15 +405,6 @@ module subroutine collision_util_bounce_one(r,v,rcom,vcom,radius) real(DP), intent(in) :: radius end subroutine collision_util_bounce_one - module subroutine collision_util_construct_constraint_system(collider, nbody_system, param, constraint_system, phase) - implicit none - class(collision_basic), intent(inout) :: collider !! Collision system object - class(base_nbody_system), intent(in) :: nbody_system !! Original Swiftest nbody system object - class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters - class(base_nbody_system), allocatable, intent(out) :: constraint_system !! Output temporary Swiftest nbody system object - character(len=*), intent(in) :: phase !! One of "before" or "after", indicating which phase of the calculation this needs to be done - end subroutine collision_util_construct_constraint_system - module subroutine collision_util_dealloc_fragments(self) implicit none class(collision_fragments), intent(inout) :: self diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index e5c4a86b6..dadde1dfd 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -33,8 +33,8 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa integer(I4B), dimension(2) :: nchild integer(I4B) :: i, j, nimpactors, idx_child real(DP), dimension(2) :: volume, density - real(DP) :: mchild, volchild, rrel_mag, rlim, dt, mtot - real(DP), dimension(NDIM) :: xc, vc, xcom, vcom, xchild, vchild, xcrossv, rrel, rrel_unit, dr + real(DP) :: mchild, volchild, rrel_mag, rlim, mtot, vdotr + real(DP), dimension(NDIM) :: xc, vc, xcom, vcom, xchild, vchild, xcrossv, rrel, vrel, rrel_unit, vrel_unit, dr real(DP), dimension(NDIM,2) :: mxc, vcc select type(nbody_system) @@ -151,7 +151,11 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa rrel_mag = .mag. rrel if (rrel_mag < rlim) then rrel_unit = .unit.rrel - dr(:) = (1.0_DP + 2*epsilon(1.0_DP)) * (rlim - rrel_mag) * rrel_unit(:) + vrel = impactors%vb(:,2) - impactors%vb(:,1) + vrel_unit = .unit.vrel + vdotr = dot_product(vrel_unit, rrel) + dr(:) = -(vdotr - sign(1.0_DP, vdotr) * sqrt(rlim**2 - rrel_mag**2 + vdotr**2)) * vrel_unit(:) + dr(:) = (1.0_DP + 2*epsilon(1.0_DP)) * dr(:) impactors%rb(:,1) = impactors%rb(:,1) - dr(:) * impactors%mass(2) / mtot impactors%rb(:,2) = impactors%rb(:,2) + dr(:) * impactors%mass(1) / mtot rrel = impactors%rb(:,2) - impactors%rb(:,1) @@ -684,4 +688,4 @@ module subroutine collision_resolve_pltp(self, nbody_system, param, t, dt, irec) return end subroutine collision_resolve_pltp -end submodule s_collision_resolve \ No newline at end of file +end submodule s_collision_resolve diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 2b468c6f7..19a688c8a 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -84,82 +84,6 @@ module subroutine collision_util_bounce_one(r,v,rcom,vcom,radius) return end subroutine collision_util_bounce_one - module subroutine collision_util_construct_constraint_system(collider, nbody_system, param, constraint_system, phase) - !! Author: David A. Minton - !! - !! Constructs a temporary system that is used to evaluate the convergence on energy and angular momentum constraints. - !! The rotations of all bodies other than those involved in the collision are set to 0 so that the changes in spin kinetic - !! energy and momentum are isolated to the collision system. - implicit none - ! Arguments - class(collision_basic), intent(inout) :: collider !! Collision system object - class(base_nbody_system), intent(in) :: nbody_system !! Original Swiftest nbody system object - class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters - class(base_nbody_system), allocatable, intent(out) :: constraint_system !! Output temporary Swiftest nbody system object - character(len=*), intent(in) :: phase !! One of "before" or "after", indicating which phase of the calculation this needs to be done - ! Internals - integer(I4B) :: i, status - class(swiftest_nbody_system), allocatable :: tmpsys - - select type(nbody_system) - class is (swiftest_nbody_system) - select type(param) - class is (swiftest_parameters) - associate(pl => nbody_system%pl, npl => nbody_system%pl%nbody, cb => nbody_system%cb, nfrag => collider%fragments%nbody) - - ! Set up a new temporary system based on the original - allocate(tmpsys, mold=nbody_system) - allocate(tmpsys%cb, source=cb) - allocate(tmpsys%pl, source=pl) - allocate(tmpsys%collider, source=collider) - call tmpsys%collider%set_original_scale() - - ! Remove spins and velocities from all bodies other than the new fragments so that we can isolate the kinetic energy and momentum of the collision system, but still be able to compute - ! the potential energy correctly - tmpsys%cb%Gmass = 0.0_DP - tmpsys%cb%mass = 0.0_DP - tmpsys%cb%rot(:) = 0.0_DP - tmpsys%pl%rot(:,:) = 0.0_DP - tmpsys%pl%vb(:,:) = 0.0_DP - - if (phase == "before") then - ! Put back the spins and velocities of the colliding bodies to compute pre-impact KE and L - tmpsys%pl%rot(:,collider%impactors%id(:)) = pl%rot(:,collider%impactors%id(:)) - tmpsys%pl%vb(:,collider%impactors%id(:)) = pl%vb(:,collider%impactors%id(:)) - else if (phase == "after") then - allocate(tmpsys%pl_adds, mold=pl) - allocate(tmpsys%pl_discards, mold=pl) - associate(impactors => tmpsys%collider%impactors, fragments => tmpsys%collider%fragments) ! Be sure to select the temporary version because its unit system has been updated - ! Update barycentric vector values - do concurrent(i = 1:nfrag) - fragments%rb(:,i) = fragments%rc(:,i) + impactors%rbcom(:) - fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:) - end do - - select case(impactors%regime) - case(COLLRESOLVE_REGIME_DISRUPTION) - status = DISRUPTED - case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - status = SUPERCATASTROPHIC - case(COLLRESOLVE_REGIME_HIT_AND_RUN) - status = HIT_AND_RUN_DISRUPT - end select - end associate - - call collision_resolve_mergeaddsub(tmpsys, param, nbody_system%t, status) - deallocate(tmpsys%collider%before) - deallocate(tmpsys%collider%after) - call tmpsys%pl%rearray(tmpsys, param) - end if - call move_alloc(tmpsys, constraint_system) - - end associate - end select - end select - - return - end subroutine collision_util_construct_constraint_system - module subroutine collision_util_get_idvalues_snapshot(self, idvals) !! author: David A. Minton @@ -218,9 +142,6 @@ module subroutine collision_util_get_energy_and_momentum(self, nbody_system, par !! Author: David A. Minton !! !! Calculates total system energy in either the pre-collision outcome state (phase = "before") or the post-collision outcome state (lbefore = .false.) - !! This subrourtine works by building a temporary internal massive body object out of the non-excluded bodies and optionally with fragments appended. - !! This will get passed to the energy calculation subroutine so that energy is computed exactly the same way is it is in the main program. - !! This will temporarily expand the massive body object in a temporary system object called constraint_system to feed it into symba_energy implicit none ! Arguments class(collision_basic), intent(inout) :: self !! Encounter collision system object @@ -228,7 +149,6 @@ module subroutine collision_util_get_energy_and_momentum(self, nbody_system, par class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters character(len=*), intent(in) :: phase !! One of "before" or "after", indicating which phase of the calculation this needs to be done ! Internals - class(base_nbody_system), allocatable :: constraint_system integer(I4B) :: i, phase_val select case(phase) @@ -246,40 +166,44 @@ module subroutine collision_util_get_energy_and_momentum(self, nbody_system, par select type(param) class is (swiftest_parameters) associate(fragments => self%fragments, impactors => self%impactors, nfrag => self%fragments%nbody, pl => nbody_system%pl, cb => nbody_system%cb) - call collision_util_construct_constraint_system(self, nbody_system, param, constraint_system, phase) - select type(constraint_system) - class is (swiftest_nbody_system) - call constraint_system%get_energy_and_momentum(param) - self%L_orbit(:,phase_val) = constraint_system%L_orbit(:) / self%Lscale - self%L_spin(:,phase_val) = constraint_system%L_spin(:) / self%Lscale - self%L_total(:,phase_val) = constraint_system%L_total(:) / self%Lscale - self%ke_orbit(phase_val) = constraint_system%ke_orbit / self%Escale - self%ke_spin(phase_val) = constraint_system%ke_spin / self%Escale - self%pe(phase_val) = constraint_system%pe / self%Escale - self%be(phase_val) = constraint_system%be / self%Escale - self%te(phase_val) = constraint_system%te / self%Escale - if (phase_val == 1) then - do concurrent(i = 1:2) - impactors%ke_orbit(i) = 0.5_DP * impactors%mass(i) * dot_product(impactors%vc(:,i), impactors%vc(:,i)) - impactors%ke_spin(i) = 0.5_DP * impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i) * dot_product(impactors%rot(:,i), impactors%rot(:,i)) - impactors%L_orbit(:,i) = impactors%mass(i) * impactors%rc(:,i) .cross. impactors%vc(:,i) - impactors%L_spin(:,i) = impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i) * impactors%rot(:,i) - end do - else if (phase_val == 2) then - do concurrent(i = 1:nfrag) - fragments%ke_orbit(i) = 0.5_DP * fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i)) - fragments%ke_spin(i) = 0.5_DP * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * dot_product(fragments%rot(:,i), fragments%rot(:,i)) - fragments%L_orbit(:,i) = fragments%mass(i) * fragments%rc(:,i) .cross. fragments%vc(:,i) - fragments%L_spin(:,i) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * fragments%rot(:,i) - end do - call swiftest_util_get_potential_energy(nfrag, [(.true., i = 1, nfrag)], constraint_system%cb%Gmass, fragments%Gmass, fragments%mass, fragments%rb, fragments%pe) - fragments%be = sum(-3*fragments%Gmass(1:nfrag)*fragments%mass(1:nfrag)/(5*fragments%radius(1:nfrag))) - fragments%L_orbit_tot(:) = sum(fragments%L_orbit(:,1:nfrag),dim=2) - fragments%L_spin_tot(:) = sum(fragments%L_spin(:,1:nfrag),dim=2) - fragments%ke_orbit_tot = sum(fragments%ke_orbit(1:nfrag)) - fragments%ke_spin_tot = sum(fragments%ke_spin(1:nfrag)) - end if - end select + if (phase_val == 1) then + do concurrent(i = 1:2) + impactors%ke_orbit(i) = 0.5_DP * impactors%mass(i) * dot_product(impactors%vc(:,i), impactors%vc(:,i)) + impactors%ke_spin(i) = 0.5_DP * impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i) * dot_product(impactors%rot(:,i), impactors%rot(:,i)) + impactors%be(i) = -3 * impactors%Gmass(i) * impactors%mass(i) / (5 * impactors%radius(i)) + impactors%L_orbit(:,i) = impactors%mass(i) * impactors%rc(:,i) .cross. impactors%vc(:,i) + impactors%L_spin(:,i) = impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i) * impactors%rot(:,i) + end do + self%L_orbit(:,phase_val) = sum(impactors%L_orbit(:,1:2),dim=2) + self%L_spin(:,phase_val) = sum(impactors%L_spin(:,1:2),dim=2) + self%L_total(:,phase_val) = self%L_orbit(:,phase_val) + self%L_spin(:,phase_val) + self%ke_orbit(phase_val) = sum(impactors%ke_orbit(1:2)) + self%ke_spin(phase_val) = sum(impactors%ke_spin(1:2)) + self%be(phase_val) = sum(impactors%be(1:2)) + call swiftest_util_get_potential_energy(2, [(.true., i = 1, 2)], 0.0_DP, impactors%Gmass, impactors%mass, impactors%rb, self%pe(phase_val)) + self%te(phase_val) = self%ke_orbit(phase_val) + self%ke_spin(phase_val) + self%be(phase_val) + self%pe(phase_val) + else if (phase_val == 2) then + do concurrent(i = 1:nfrag) + fragments%ke_orbit(i) = 0.5_DP * fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i)) + fragments%ke_spin(i) = 0.5_DP * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * dot_product(fragments%rot(:,i), fragments%rot(:,i)) + fragments%L_orbit(:,i) = fragments%mass(i) * fragments%rc(:,i) .cross. fragments%vc(:,i) + fragments%L_spin(:,i) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * fragments%rot(:,i) + end do + call swiftest_util_get_potential_energy(nfrag, [(.true., i = 1, nfrag)], 0.0_DP, fragments%Gmass, fragments%mass, fragments%rb, fragments%pe) + fragments%be = sum(-3*fragments%Gmass(1:nfrag)*fragments%mass(1:nfrag)/(5*fragments%radius(1:nfrag))) + fragments%L_orbit_tot(:) = sum(fragments%L_orbit(:,1:nfrag),dim=2) + fragments%L_spin_tot(:) = sum(fragments%L_spin(:,1:nfrag),dim=2) + fragments%ke_orbit_tot = sum(fragments%ke_orbit(1:nfrag)) + fragments%ke_spin_tot = sum(fragments%ke_spin(1:nfrag)) + self%L_orbit(:,phase_val) = fragments%L_orbit_tot(:) + self%L_spin(:,phase_val) = fragments%L_spin_tot(:) + self%L_total(:,phase_val) = self%L_orbit(:,phase_val) + self%L_spin(:,phase_val) + self%ke_orbit(phase_val) = fragments%ke_orbit_tot + self%ke_spin(phase_val) = fragments%ke_spin_tot + self%be(phase_val) = fragments%be + self%pe(phase_val) = fragments%pe + self%te(phase_val) = self%ke_orbit(phase_val) + self%ke_spin(phase_val) + self%be(phase_val) + self%pe(phase_val) + end if end associate end select @@ -909,8 +833,8 @@ module subroutine collision_util_set_natural_scale_factors(self) associate(collider => self, fragments => self%fragments, impactors => self%impactors) ! Set primary scale factors (mass, length, and time) based on the impactor properties at the time of collision - collider%mscale = minval(fragments%mass(:)) - collider%dscale = minval(fragments%radius(:)) + collider%mscale = minval(impactors%mass(:)) + collider%dscale = minval(impactors%radius(:)) vesc = sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:))) collider%tscale = collider%dscale / vesc diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index b9b38fc4c..6ee9c36bd 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -124,7 +124,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur real(DP) :: dE real(DP), dimension(NDIM) :: dL character(len=STRMAX) :: message - real(DP), parameter :: fail_scale_initial = 1.001_DP + real(DP), parameter :: fail_scale_initial = 1.01_DP integer(I4B) :: nfrag_start ! The minimization and linear solvers can sometimes lead to floating point exceptions. Rather than halting the code entirely if this occurs, we @@ -284,6 +284,52 @@ module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t) end subroutine fraggle_generate_hitandrun + module subroutine fraggle_generate_merge(self, nbody_system, param, t) + !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Merge massive bodies in any collisional system. If the rotation is too high, switch to hit and run. + !! + !! Adapted from David E. Kaufmann's Swifter routines symba_merge_pl.f90 and symba_discard_merge_pl.f90 + !! + !! Adapted from Hal Levison's Swift routines symba5_merge.f and discard_mass_merge.f + implicit none + class(collision_fraggle), intent(inout) :: self !! Fraggle system object + class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! The time of the collision + + ! Internals + integer(I4B) :: i, j + real(DP), dimension(NDIM) :: L_spin_new, Ip, rot + real(DP) :: rotmag, mass, volume, radius + + select type(nbody_system) + class is (swiftest_nbody_system) + associate(impactors => self%impactors) + mass = sum(impactors%mass(:)) + volume = 4._DP / 3._DP * PI * sum(impactors%radius(:)**3) + radius = (3._DP * volume / (4._DP * PI))**(THIRD) + do concurrent(i = 1:NDIM) + Ip(i) = sum(impactors%mass(:) * impactors%Ip(i,:)) + L_spin_new(i) = sum(impactors%L_orbit(i,:) + impactors%L_spin(i,:)) + end do + Ip(:) = Ip(:) / mass + rot(:) = L_spin_new(:) / (Ip(3) * mass * radius**2) + rotmag = .mag.rot(:) + if (rotmag < self%max_rot) then + call self%collision_basic%merge(nbody_system, param, t) + else + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Merger would break the spin barrier. Converting to pure hit and run" ) + impactors%mass_dist(1:2) = impactors%mass(1:2) + call self%hitandrun(nbody_system, param, t) + end if + + end associate + end select + return + end subroutine fraggle_generate_merge + + module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailure) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! @@ -309,10 +355,8 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu integer(I4B), parameter :: MAXLOOP = 10000 real(DP), parameter :: cloud_size_scale_factor = 3.0_DP ! Scale factor to apply to the size of the cloud relative to the distance from the impact point. ! A larger value puts more space between fragments initially - real(DP) :: rbuffer ! Body radii are inflated by this scale factor to prevent secondary collisions - real(DP), parameter :: rbuffer_max = 1.2_DP + real(DP), parameter :: rbuffer = 1.01_DP ! Body radii are inflated by this scale factor to prevent secondary collisions real(DP), parameter :: pack_density = 0.5236_DP ! packing density of loose spheres - rbuffer = 1.01_DP associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody, & pl => nbody_system%pl, tp => nbody_system%tp, npl => nbody_system%pl%nbody, ntp => nbody_system%tp%nbody) @@ -328,20 +372,17 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu rdistance = 0.5_DP * sum(impactors%radius(:)) istart = 3 else - rdistance = 2 * impactors%radius(2) + rdistance = impactors%radius(2) istart = 3 end if - if (lhitandrun) then - mass_rscale(:) = 1.0_DP - else - ! Give the fragment positions a random value that is scaled with fragment mass so that the more massive bodies tend to be closer to the impact point - ! Later, velocities will be scaled such that the farther away a fragment is placed from the impact point, the higher will its velocity be. - call random_number(mass_rscale) - mass_rscale(:) = (mass_rscale(:) + 1.0_DP) / 2 - mass_rscale(:) = mass_rscale(:) * (fragments%mtot / fragments%mass(1:nfrag))**(0.125_DP) ! The power is arbitrary. It just gives the velocity a small mass dependence - mass_rscale(:) = mass_rscale(:) / maxval(mass_rscale(:)) - end if + mass_rscale(1:istart-1) = 1.0_DP + ! Give the fragment positions a random value that is scaled with fragment mass so that the more massive bodies tend to be closer to the impact point + ! Later, velocities will be scaled such that the farther away a fragment is placed from the impact point, the higher will its velocity be. + call random_number(mass_rscale(istart:nfrag)) + mass_rscale(istart:nfrag) = (mass_rscale(istart:nfrag) + 1.0_DP) / 2 + mass_rscale(istart:nfrag) = mass_rscale(istart:nfrag) * (sum(fragments%mass(istart:nfrag)) / fragments%mass(istart:nfrag))**(0.125_DP) ! The power is arbitrary. It just gives the velocity a small mass dependence + mass_rscale(istart:nfrag) = mass_rscale(istart:nfrag) / minval(mass_rscale(istart:nfrag)) loverlap(:) = .true. do loop = 1, MAXLOOP @@ -355,9 +396,8 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu else ! Keep the first and second bodies at approximately their original location, but so as not to be overlapping fragment_cloud_center(:,1) = impactors%rcimp(:) - rbuffer * max(fragments%radius(1),impactors%radius(1)) * impactors%y_unit(:) fragment_cloud_center(:,2) = impactors%rcimp(:) + rbuffer * max(fragments%radius(2),impactors%radius(2)) * impactors%y_unit(:) - fragment_cloud_radius(:) = cloud_size_scale_factor * rdistance - fragments%rc(:,1) = fragment_cloud_center(:,1) - fragments%rc(:,2) = fragment_cloud_center(:,2) + fragment_cloud_radius(:) = rdistance / pack_density + fragments%rc(:,1:2) = fragment_cloud_center(:,1:2) end if do i = 1, nfrag @@ -365,6 +405,8 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu call random_number(phi(i)) call random_number(theta(i)) call random_number(u(i)) + phi(i) = TWOPI * phi(i) + theta(i) = asin(2 * theta(i) - 1.0_DP) end if end do @@ -372,12 +414,12 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu do concurrent(i = istart:nfrag, loverlap(i)) j = fragments%origin_body(i) - ! Make a random cloud - phi(i) = TWOPI * phi(i) - theta(i) = acos(2 * theta(i) - 1.0_DP) - ! Scale the cloud size - fragments%rmag(i) = fragment_cloud_radius(j) * mass_rscale(i) * u(i)**(THIRD) + if (lhitandrun) then + fragments%rmag(i) = fragment_cloud_radius(j) * u(i)**(THIRD) + else + fragments%rmag(i) = fragment_cloud_radius(j) * mass_rscale(i) * u(i)**(THIRD) + end if ! Position the fragment in a random point within the cloud fragments%rc(1,i) = fragments%rmag(i) * sin(theta(i)) * cos(phi(i)) @@ -385,11 +427,16 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu fragments%rc(3,i) = fragments%rmag(i) * cos(theta(i)) ! Shift to the cloud center coordinates - fragments%rc(:,i) = fragments%rc(:,i) + fragment_cloud_center(:,j) ! Stretch out the hit and run cloud along the flight trajectory if (lhitandrun) then - fragments%rc(:,i) = fragments%rc(:,i) + cloud_size_scale_factor * rdistance * u(i) * impactors%bounce_unit(:) + fragments%rc(:,i) = fragments%rc(:,i) * (1.0_DP + 2 * fragment_cloud_radius(j) * mass_rscale(i) * impactors%bounce_unit(:)) + end if + + fragments%rc(:,i) = fragments%rc(:,i) + fragment_cloud_center(:,j) + + if (lhitandrun) then + fragments%rc(:,i) = fragments%rc(:,i) + 2 * fragment_cloud_radius(j) * mass_rscale(i) * impactors%bounce_unit(:) ! Shift the stretched cloud downrange else ! Make sure that the fragments are positioned away from the impact point direction = dot_product(fragments%rc(:,i) - impactors%rcimp(:), fragment_cloud_center(:,j) - impactors%rcimp(:)) @@ -415,14 +462,17 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu fragments%rmag(:) = .mag. fragments%rc(:,:) ! Check for any overlapping bodies. - loverlap(:) = .false. - do j = 1, nfrag + do j = nfrag, 1, -1 + if (.not.loverlap(j)) cycle + loverlap(j) = .false. ! Check for overlaps between fragments - do i = j + 1, nfrag + do i = 1,nfrag + if (i == j) cycle dis = .mag.(fragments%rc(:,j) - fragments%rc(:,i)) - loverlap(i) = loverlap(i) .or. (dis <= rbuffer * (fragments%radius(i) + fragments%radius(j))) + loverlap(j) = (dis <= rbuffer * (fragments%radius(i) + fragments%radius(j))) + if (loverlap(j)) exit end do - ! Check for overlaps with existing bodies that are not involved in the collision + if (loverlap(j)) cycle do i = 1, npl if (any(impactors%id(:) == i)) cycle dis = .mag. (fragments%rc(:,j) - (pl%rb(:,i) / collider%dscale - impactors%rbcom(:))) @@ -430,7 +480,6 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu end do end do rdistance = rdistance * collider%fail_scale - rbuffer = min(rbuffer * collider%fail_scale, rbuffer_max) end do lfailure = any(loverlap(:)) @@ -533,7 +582,8 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters logical, intent(out) :: lfailure !! Did the velocity computation fail? ! Internals - real(DP), parameter :: ENERGY_SUCCESS_METRIC = 1.0e-3_DP !! Relative energy error to accept as a success (success also must be energy-losing in addition to being within the metric amount) + real(DP), parameter :: ENERGY_SUCCESS_METRIC = 0.1_DP !! Relative energy error to accept as a success (success also must be energy-losing in addition to being within the metric amount) + real(DP), parameter :: ENERGY_CONVERGENCE_TOL = 1e-3_DP !! Relative change in error before giving up on energy convergence real(DP) :: MOMENTUM_SUCCESS_METRIC = 10*epsilon(1.0_DP) !! Relative angular momentum error to accept as a success (should be *much* stricter than energy) integer(I4B) :: i, j, loop, try, istart, nfrag, nsteps, nsteps_best, posloop logical :: lhitandrun, lsupercat @@ -546,8 +596,8 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu real(DP), parameter :: hitandrun_vscale = 0.25_DP real(DP) :: vmin_guess real(DP) :: vmax_guess - integer(I4B), parameter :: MAXLOOP = 100 - integer(I4B), parameter :: MAXTRY = 100 + integer(I4B), parameter :: MAXLOOP = 25 + integer(I4B), parameter :: MAXTRY = 10 integer(I4B), parameter :: MAXANGMTM = 1000 class(collision_fraggle), allocatable :: collider_local character(len=STRMAX) :: message @@ -681,11 +731,11 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu if (all(dL_metric(:) <= 1.0_DP)) exit angmtm do i = istart, fragments%nbody - dL(:) = -L_residual(:) / (fragments%nbody - istart + 1) + dL(:) = -L_residual(:) * fragments%mass(i) / sum(fragments%mass(istart:fragments%nbody)) call collision_util_velocity_torque(dL, fragments%mass(i), fragments%rc(:,i), fragments%vc(:,i)) - call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) - fragments%vmag(i) = .mag.fragments%vc(:,i) end do + call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) + fragments%vmag(:) = .mag.fragments%vc(:,:) end do angmtm call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) @@ -713,7 +763,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu dE_metric = abs(E_residual) / impactors%Qloss end if if ((dE_best < 0.0_DP) .and. (dE_metric <= ENERGY_SUCCESS_METRIC)) exit outer - if (dE_conv < ENERGY_SUCCESS_METRIC) exit inner + if (dE_conv < ENERGY_CONVERGENCE_TOL) exit inner end if ! Remove a constant amount of velocity from the bodies so we don't shift the center of mass and screw up the momentum diff --git a/src/fraggle/fraggle_module.f90 b/src/fraggle/fraggle_module.f90 index 7cdfbd925..8d118fbaf 100644 --- a/src/fraggle/fraggle_module.f90 +++ b/src/fraggle/fraggle_module.f90 @@ -18,9 +18,10 @@ module fraggle type, extends(collision_basic) :: collision_fraggle real(DP) :: fail_scale !! Scale factor to apply to distance values in the position model when overlaps occur. contains - procedure :: disrupt => fraggle_generate_disrupt !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. procedure :: generate => fraggle_generate !! A simple disruption models that does not constrain energy loss in collisions + procedure :: disrupt => fraggle_generate_disrupt !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. procedure :: hitandrun => fraggle_generate_hitandrun !! Generates either a pure hit and run, or one in which the runner is disrupted + procedure :: merge => fraggle_generate_merge !! Merges bodies unless the rotation would be too high, then it switches to pure hit and run. procedure :: set_mass_dist => fraggle_util_set_mass_dist !! Sets the distribution of mass among the fragments depending on the regime type procedure :: restructure => fraggle_util_restructure !! Restructures the fragment distribution after a failure to converge on a solution end type collision_fraggle @@ -28,7 +29,7 @@ module fraggle interface module subroutine fraggle_generate(self, nbody_system, param, t) implicit none - class(collision_fraggle), intent(inout) :: self !! Fraggle fragment system object + class(collision_fraggle), intent(inout) :: self !! Fraggle fragment system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! The time of the collision @@ -36,7 +37,7 @@ end subroutine fraggle_generate module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailure) implicit none - class(collision_fraggle), intent(inout) :: self !! Fraggle system object the outputs will be the fragmentation + class(collision_fraggle), intent(inout) :: self !! Fraggle system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Time of collision @@ -45,12 +46,20 @@ end subroutine fraggle_generate_disrupt module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t) implicit none - class(collision_fraggle), intent(inout) :: self !! Collision system object + class(collision_fraggle), intent(inout) :: self !! Fraggle system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions real(DP), intent(in) :: t !! Time of collision end subroutine fraggle_generate_hitandrun + module subroutine fraggle_generate_merge(self, nbody_system, param, t) + implicit none + class(collision_fraggle), intent(inout) :: self !! Fraggle system object + class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! The time of the collision + end subroutine fraggle_generate_merge + module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailure) implicit none class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 1f1e51b72..4daef8d40 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -86,20 +86,20 @@ module subroutine fraggle_util_set_mass_dist(self, param) class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters ! Internals - integer(I4B) :: i, j, jproj, jtarg, nfrag, istart, nfragmax, nrem + integer(I4B) :: i, j, k, jproj, jtarg, nfrag, istart, nfragmax, nrem real(DP), dimension(2) :: volume real(DP), dimension(NDIM) :: Ip_avg - real(DP) :: mfrag, mremaining, mtot, mcumul, G, mass_noise, Mslr, mscale, y, yp, Mrat, m0 + real(DP) :: mremaining, mtot, mcumul, G, mass_noise, Mslr, x0, x1, ymid, y0, y1, y, yp, eps, Mrat real(DP), dimension(:), allocatable :: mass real(DP) :: beta integer(I4B), parameter :: MASS_NOISE_FACTOR = 5 !! The number of digits of random noise that get added to the minimum mass value to prevent identical masses from being generated in a single run - integer(I4B), parameter :: NFRAGMAX_UNSCALED = 3000 !! Maximum number of fragments that can be generated + integer(I4B), parameter :: NFRAGMAX_UNSCALED = 100 !! Maximum number of fragments that can be generated integer(I4B), parameter :: iMlr = 1 integer(I4B), parameter :: iMslr = 2 integer(I4B), parameter :: iMrem = 3 integer(I4B), parameter :: NFRAGMIN = iMrem + 2 !! Minimum number of fragments that can be generated integer(I4B), dimension(:), allocatable :: ind - integer(I4B), parameter :: MAXLOOP = 10000 + integer(I4B), parameter :: MAXLOOP = 20 logical :: flipper associate(impactors => self%impactors, min_mfrag => self%min_mfrag) @@ -120,29 +120,19 @@ module subroutine fraggle_util_set_mass_dist(self, param) select case(impactors%regime) case(COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC, COLLRESOLVE_REGIME_HIT_AND_RUN) ! The first two bins of the mass_dist are the largest and second-largest fragments that came out of collision_regime. - ! The remainder from the third bin will be distributed among nfrag-2 bodies. The following code will determine nfrag based on - ! the limits bracketed above and the model size distribution of fragments. - ! Check to see if our size distribution would give us a smaller number of fragments than the maximum number + ! The remainder from the third bin will be distributed among nfrag-2 bodies. - ! ! Add a small amount of noise to the last digits of the minimum mass value so that multiple fragments don't get generated with identical mass values + !Add a small amount of noise to the last digits of the minimum mass value so that multiple fragments don't get generated with identical mass values call random_number(mass_noise) mass_noise = 1.0_DP + mass_noise * epsilon(1.0_DP) * 10**(MASS_NOISE_FACTOR) min_mfrag = (param%min_GMfrag / param%GU) * mass_noise mremaining = impactors%mass_dist(iMrem) nfrag = iMrem - 1 - beta = 2.85_DP ! From Leinhardt & Stewart (2012) Mslr = impactors%mass_dist(iMslr) - - mfrag = Mslr - do while (mremaining > 0.0_DP) - mfrag = (nfrag)**(-3._DP / BETA) * impactors%mass_dist(iMslr) - mfrag = max(mfrag, min_mfrag) - mremaining = mremaining - mfrag - nfrag = nfrag + 1 - end do nfragmax = ceiling(NFRAGMAX_UNSCALED / param%nfrag_reduction) - nfrag = max(min(ceiling(nfrag / param%nfrag_reduction), nfragmax), NFRAGMIN) + Mrat = max((mremaining + Mslr) / min_mfrag, 1.0_DP) + nfrag = max(ceiling(nfragmax * (log(Mrat) + 1.0_DP)), NFRAGMIN) call self%setup_fragments(nfrag) @@ -173,26 +163,58 @@ module subroutine fraggle_util_set_mass_dist(self, param) if (Mslr == min_mfrag) Mslr = Mslr + impactors%mass_dist(iMrem) / nfrag mremaining = impactors%mass_dist(iMrem) - ! The mass will be distributed evenly in logspace between the second-largest remnant and the minimum mass + ! The mass will be distributed in a power law where N>M=(M/Mslr)**(-beta/3) ! Use Newton's method solver to get the logspace slope of the mass function - Mrat = mremaining / min_mfrag - nrem = nfrag - 2 - mscale = 1.0_DP - 1.0_DP/Mrat + Mrat = (mremaining + Mslr) / Mslr + x0 = -3.0_DP + x1 = 3.0_DP + do j = 1, MAXLOOP + y0 = Mrat - 1.0_DP + y1 = Mrat - 1.0_DP + do i = iMrem, nfrag + y0 = y0 - (i-1)**(-x0/3.0_DP) + y1 = y1 - (i-1)**(-x1/3.0_DP) + end do + if (y0*y1 < 0.0_DP) exit + x1 = x1 * 1.6_DP + x0 = x0 / 1.6_DP + end do + + ! Find the mass scaling factor with bisection do j = 1, MAXLOOP - y = Mrat - (1.0_DP - mscale**nrem)/(1.0_DP - mscale) - yp = (mscale + mscale**nrem * ((nrem - 1) * mscale - nrem)) / (mscale * (mscale - 1.0_DP)**2) - m0 = mscale - mscale = m0 + y/yp - if (abs((mscale - m0)/mscale) < epsilon(1.0_DP)) exit + beta = 0.5_DP * (x0 + x1) + ymid = Mrat - 1.0_DP + do i = iMrem, nfrag + ymid = ymid - (i-1)**(-beta/3.0_DP) + end do + if (abs((y0 - ymid)/y0) < 1e-12_DP) exit + if (y0 * ymid < 0.0_DP) then + x1 = beta + y1 = ymid + else if (y1 * ymid < 0.0_DP) then + x0 = beta + y0 = ymid + else + exit + end if + ! Finish with Newton if we still haven't made it + if (j == MAXLOOP) then + do k = 1, MAXLOOP + y = Mrat - 1.0_DP + yp = 0.0_DP + do i = iMrem, nfrag + y = y - (i-1)**(-beta/3.0_DP) + if (i > 3) yp = yp - (i-1)**(-beta/3.0_DP) * log(real(i-1,kind=DP)) + end do + eps = y/yp + if (abs(eps) < epsilon(1.0_DP) * beta) exit + beta = beta + eps + end do + end if end do - Mslr = impactors%mass_dist(iMslr) - mass(iMslr) = Mslr - mfrag = min_mfrag do i = iMrem,nfrag - mass(i) = mfrag - mremaining = mremaining - mfrag - mfrag = max(mfrag * mscale, min_mfrag) + mass(i) = Mslr * (i-1)**(-beta/3.0_DP) end do ! There may still be some small residual due to round-off error. If so, simply add it to the last bin of the mass distribution. diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index 744c204ed..40c357fc2 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -909,6 +909,7 @@ module subroutine swiftest_io_netcdf_open(self, param, readonly) ! end if ! Optional variables The User Doesn't Need to Know About + status = nf90_inq_varid(nc%id, nc%mass_varname, nc%mass_varid) status = nf90_inq_varid(nc%id, nc%rhill_varname, nc%rhill_varid) param%lrhill_present = (status == NF90_NOERR) status = nf90_inq_varid(nc%id, nc%npl_varname, nc%npl_varid) @@ -1649,23 +1650,23 @@ module subroutine swiftest_io_netcdf_write_frame_cb(self, nc, param) associate(tslot => nc%tslot) call self%write_info(nc, param) - call netcdf_io_check( nf90_set_fill(nc%id, NF90_NOFILL, old_mode), "netcdf_io_write_frame_cb nf90_set_fill" ) + call netcdf_io_check( nf90_set_fill(nc%id, NF90_NOFILL, old_mode), "swiftest_io_netcdf_write_frame_cb nf90_set_fill" ) call nc%find_idslot(self%id, idslot) - call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, self%id, start=[idslot]), "netcdf_io_write_frame_cb nf90_put_var cb id_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%status_varid, ACTIVE, start=[idslot, tslot]), "netcdf_io_write_frame_cb nf90_put_var cb id_varid" ) - - call netcdf_io_check( nf90_put_var(nc%id, nc%Gmass_varid, self%Gmass, start=[idslot, tslot]), "netcdf_io_write_frame_cb nf90_put_var cb Gmass_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%mass_varid, self%mass, start=[idslot, tslot]), "netcdf_io_write_frame_cb nf90_put_var cb mass_varid" ) - if (param%lclose) call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, self%radius, start=[idslot, tslot]), "netcdf_io_write_frame_cb nf90_put_var cb radius_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%j2rp2_varid, self%j2rp2, start=[tslot]), "netcdf_io_write_frame_cb nf90_put_var cb j2rp2_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%j4rp4_varid, self%j4rp4, start=[tslot]), "netcdf_io_write_frame_cb nf90_put_var cb j4rp4_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, self%id, start=[idslot]), "swiftest_io_netcdf_write_frame_cb nf90_put_var cb id_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%status_varid, ACTIVE, start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_cb nf90_put_var cb id_varid" ) + + call netcdf_io_check( nf90_put_var(nc%id, nc%Gmass_varid, self%Gmass, start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_cb nf90_put_var cb Gmass_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%mass_varid, self%mass, start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_cb nf90_put_var cb mass_varid" ) + if (param%lclose) call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, self%radius, start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_cb nf90_put_var cb radius_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%j2rp2_varid, self%j2rp2, start=[tslot]), "swiftest_io_netcdf_write_frame_cb nf90_put_var cb j2rp2_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%j4rp4_varid, self%j4rp4, start=[tslot]), "swiftest_io_netcdf_write_frame_cb nf90_put_var cb j4rp4_varid" ) if (param%lrotation) then - call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, self%Ip(:), start=[1, idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_cb nf90_put_var cb Ip_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, self%rot(:) * RAD2DEG, start=[1, idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_cby nf90_put_var cb rot_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, self%Ip(:), start=[1, idslot, tslot], count=[NDIM,1,1]), "swiftest_io_netcdf_write_frame_cb nf90_put_var cb Ip_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, self%rot(:) * RAD2DEG, start=[1, idslot, tslot], count=[NDIM,1,1]), "swiftest_io_netcdf_write_frame_cby nf90_put_var cb rot_varid" ) end if - call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode), "netcdf_io_write_frame_cb nf90_set_fill old_mode" ) + call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode), "swiftest_io_netcdf_write_frame_cb nf90_set_fill old_mode" ) end associate return