From 85bcfd7fb92a146f8cecb325ddbda88653d585ba Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Mon, 13 Feb 2023 12:00:38 -0500 Subject: [PATCH 01/38] combined descriptions of ASCII parameter options and simulation kwargs --- README.md | 2 +- README_tables/param_options.md | 57 ---------------- README_tables/simulation_kwargs.md | 103 +++++++++++++++++------------ 3 files changed, 60 insertions(+), 102 deletions(-) delete mode 100644 README_tables/param_options.md diff --git a/README.md b/README.md index 75b22f0bc..2dd7969f7 100644 --- a/README.md +++ b/README.md @@ -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: 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 From b6fe17158a61860539c7e5e7946036e798454912 Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Mon, 13 Feb 2023 12:05:51 -0500 Subject: [PATCH 02/38] de-linked a few non-links --- README.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 2dd7969f7..109789768 100644 --- a/README.md +++ b/README.md @@ -195,7 +195,7 @@ 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). @@ -296,8 +296,8 @@ The number and type of output files generated by Swiftest depends on the input p - **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```). +- **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```. From 4e2c47cf2c58f725aa9121520539357a3e5d279c Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Mon, 13 Feb 2023 12:11:02 -0500 Subject: [PATCH 03/38] fixed some language --- README.md | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/README.md b/README.md index 109789768..d3ffcc335 100644 --- a/README.md +++ b/README.md @@ -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. @@ -199,7 +199,7 @@ After all desired parameters have been set, the parameters can be saved to the * 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: @@ -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,7 +295,7 @@ 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. +- **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```). @@ -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. --- From 167fddd0878ff46a3d5712919fe5b7e501a2f116 Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Mon, 13 Feb 2023 12:11:50 -0500 Subject: [PATCH 04/38] re-titled document --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index d3ffcc335..2a1863fcf 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 From 34582285c55ea61d221c8160f12c2f886a18ae1f Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Mon, 13 Feb 2023 13:04:27 -0500 Subject: [PATCH 05/38] adjusted some FAQs --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 2a1863fcf..bee2a181e 100644 --- a/README.md +++ b/README.md @@ -522,11 +522,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. From c0ac7d7e09b4cf9097aebd5dda3f53c9aa4fadbd Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Mon, 13 Feb 2023 13:08:31 -0500 Subject: [PATCH 06/38] updated fragmentation example --- examples/Fragmentation/swiftest_fragmentation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/Fragmentation/swiftest_fragmentation.py b/examples/Fragmentation/swiftest_fragmentation.py index 4c2368111..cf962b12e 100644 --- a/examples/Fragmentation/swiftest_fragmentation.py +++ b/examples/Fragmentation/swiftest_fragmentation.py @@ -42,7 +42,7 @@ 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, +run_arguments = {"tstart": 0.0, "tstop":1e-5, "dt": 1e-5, "istep_out": 1, "collision_model": "fraggle", "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. From f60a382918dd8775050ae74770100e35b201ffc4 Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Mon, 13 Feb 2023 13:14:45 -0500 Subject: [PATCH 07/38] deleted duplicate fragmentation example --- examples/Fragmentation/.gitignore | 1 - examples/Fragmentation/README.txt | 3 +- .../Fragmentation/swiftest_fragmentation.py | 79 ------------------- 3 files changed, 1 insertion(+), 82 deletions(-) delete mode 100644 examples/Fragmentation/swiftest_fragmentation.py 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/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 cf962b12e..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, "collision_model": "fraggle", "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() From a726b0c5a382af33c2dbb2ccb7aba4181b784b84 Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Mon, 13 Feb 2023 13:15:02 -0500 Subject: [PATCH 08/38] updated docstring header --- examples/Fragmentation/Fragmentation_Movie.py | 28 ++++++++++++++----- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 5afb9484d..b0b82d519 100755 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -11,16 +11,30 @@ """ """ -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 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. Inputs _______ -param.in : ASCII Swiftest parameter input file. -data.nc : A NetCDF file containing the simulation output. - -Returns -------- -fragmentation.mp4 : A .mp4 file of a fragmentation event. +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 From 8592f9fc92f52ea59449d3a34f88d72457a11dc3 Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Mon, 13 Feb 2023 13:22:09 -0500 Subject: [PATCH 09/38] updated docstring header again --- examples/Fragmentation/Fragmentation_Movie.py | 28 ++++++++----------- 1 file changed, 12 insertions(+), 16 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index b0b82d519..7d96c231e 100755 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -12,9 +12,7 @@ """ 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. +outputs are stored in the subdirectory named after their collisional regime. Inputs _______ @@ -22,19 +20,17 @@ 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. +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. + """ import swiftest From c317e14516d20057de41e618efe9f6cc4b126938 Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Mon, 13 Feb 2023 13:35:59 -0500 Subject: [PATCH 10/38] updated the description of the fragmentation example in the user manual --- README.md | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index bee2a181e..ce8a1f98a 100644 --- a/README.md +++ b/README.md @@ -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 2 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** From 4021d3fb26fab91f3bb3f23ec97530401dfb9abe Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Mon, 13 Feb 2023 13:55:13 -0500 Subject: [PATCH 11/38] updated the number of nodes needed for the fragmentation example --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index ce8a1f98a..c9f52f948 100644 --- a/README.md +++ b/README.md @@ -486,7 +486,7 @@ This example highlights the functionality of the Fraggle algorithm. It can be fo - A Merger Crossing the Spin Barrier - All of the Above -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 2 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. +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** From 0c9e63c3ac2cc18dd1007a9e162cab01a52bd8ee Mon Sep 17 00:00:00 2001 From: Kaustub Anand Date: Thu, 16 Feb 2023 10:32:46 -0500 Subject: [PATCH 12/38] kaustub testing 3 --- README.md | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index c9f52f948..592304422 100644 --- a/README.md +++ b/README.md @@ -570,4 +570,7 @@ 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 . + + +# kaustub testing \ No newline at end of file From b400f6b80d773041ea93d417a5cc17072281544f Mon Sep 17 00:00:00 2001 From: Kaustub Anand Date: Thu, 16 Feb 2023 11:03:27 -0500 Subject: [PATCH 13/38] kaustub untesting --- README.md | 3 --- 1 file changed, 3 deletions(-) diff --git a/README.md b/README.md index 592304422..5db472fec 100644 --- a/README.md +++ b/README.md @@ -571,6 +571,3 @@ 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 . - - -# kaustub testing \ No newline at end of file From 7f8da8ae2efd60dbec487f387a37ac0158b4809b Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 18 Feb 2023 18:19:31 -0500 Subject: [PATCH 14/38] Adjusted the number of tries to take to keep from taking forever when there is no hope --- src/fraggle/fraggle_generate.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index b9b38fc4c..667cf00be 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -533,7 +533,7 @@ 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 = 1.0e-2_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) :: 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 +546,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 From 33d38b3c31ad33226036dd974e88b43be3f7d6c5 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 18 Feb 2023 19:23:53 -0500 Subject: [PATCH 15/38] Relaxed the energy constraint so that we just need to be within 10% of the target energy loss --- src/fraggle/fraggle_generate.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 667cf00be..2f8c95e87 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -533,7 +533,7 @@ 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-2_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) :: 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 From a74e9f5eab8b2eb1061c75ca3d4041b5cd11ac4e Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 18 Feb 2023 19:42:29 -0500 Subject: [PATCH 16/38] Made a separate criteria for checking if energy solution is still converging --- src/fraggle/fraggle_generate.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 2f8c95e87..099417734 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -533,7 +533,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 = 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_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 @@ -713,7 +714,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 From eba83764923994a5109f93e46a598af2aebe11f1 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 19 Feb 2023 02:03:37 -0500 Subject: [PATCH 17/38] Improved overlap adjustment --- src/collision/collision_resolve.f90 | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) 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 From 175d83cdf68234b6274af14adb1aced3407b4e2a Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 19 Feb 2023 02:33:03 -0500 Subject: [PATCH 18/38] Got rid of the overly-complex constraint system that slowed down Fraggle when nfrag gets big --- src/collision/collision_module.f90 | 12 +-- src/collision/collision_util.f90 | 152 ++++++++--------------------- 2 files changed, 40 insertions(+), 124 deletions(-) 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_util.f90 b/src/collision/collision_util.f90 index 2b468c6f7..bdbf35fb9 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 From 91d2ad86f68be00584ef5512dd55da03edebdf89 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 19 Feb 2023 11:08:00 -0500 Subject: [PATCH 19/38] Added an argument that allows output files to be read in by Dask (using open_mfdataset), which is needed for very large data files. --- python/swiftest/swiftest/io.py | 4 +++- python/swiftest/swiftest/tool.py | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 98e6c154a..d070bbb91 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -842,7 +842,9 @@ 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) + #ds = xr.open_dataset(param['BIN_OUT'], mask_and_scale=False) + ds = xr.open_mfdataset(param['BIN_OUT'], parallel=True, engine='h5netcdf', 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/tool.py b/python/swiftest/swiftest/tool.py index 7fe2c8884..225619204 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], input_core_dims=[[dim]], kwargs={"ord": ord, "axis": -1}, dask="parallelized" ) def wrap_angle(angle): From 2b5ac9f46cc92de805d2c5b8a432a309de342bd0 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 19 Feb 2023 19:31:33 -0500 Subject: [PATCH 20/38] Added mass to variables that are inquired when opening an old file --- src/swiftest/swiftest_io.f90 | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) 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 From d374379f8fa99f0ce004616c24012a23c6f16e68 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Feb 2023 10:56:54 -0500 Subject: [PATCH 21/38] Fixed bug that was causing the number of fragments to be computed incorrectly when it grew to be very large --- src/fraggle/fraggle_util.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 1f1e51b72..a0dc08af5 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -134,14 +134,14 @@ module subroutine fraggle_util_set_mass_dist(self, param) 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) + nfragmax = ceiling(NFRAGMAX_UNSCALED / param%nfrag_reduction) + do i = 1, nfragmax + mfrag = (nfrag)**(-3._DP / BETA) * Mslr mfrag = max(mfrag, min_mfrag) mremaining = mremaining - mfrag nfrag = nfrag + 1 + if (mremaining < 0.0_DP) exit end do - nfragmax = ceiling(NFRAGMAX_UNSCALED / param%nfrag_reduction) nfrag = max(min(ceiling(nfrag / param%nfrag_reduction), nfragmax), NFRAGMIN) call self%setup_fragments(nfrag) From ff898874a530930415a7120b7bb9776b7536baea Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Feb 2023 13:23:46 -0500 Subject: [PATCH 22/38] Fixed issue of convergence on mass multiplier for SFD --- src/fraggle/fraggle_util.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index a0dc08af5..846122e73 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -89,7 +89,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) integer(I4B) :: i, j, 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) :: mfrag, mremaining, mtot, mcumul, G, mass_noise, Mslr, mscale, mscale0, y, yp, 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 @@ -177,13 +177,13 @@ module subroutine fraggle_util_set_mass_dist(self, param) ! 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 + mscale = Mrat**(1.0_DP/nrem) 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 + mscale0 = mscale + mscale = mscale0 + y/yp + if (abs((mscale - mscale0)/mscale) < epsilon(1.0_DP)) exit end do Mslr = impactors%mass_dist(iMslr) From 7073413086d11bd0c583c5c01a8f5009766ea72b Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Feb 2023 17:09:33 -0500 Subject: [PATCH 23/38] Improved convergence of mass factor in SFD calc using combination of bisection and Newton methods --- src/fraggle/fraggle_util.f90 | 50 +++++++++++++++++++++++++++--------- 1 file changed, 38 insertions(+), 12 deletions(-) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 846122e73..d1327d723 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -89,7 +89,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) integer(I4B) :: i, j, 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, mscale0, y, yp, Mrat + real(DP) :: mfrag, mremaining, mtot, mcumul, G, mass_noise, Mslr, mscale, x0, x1, xmid, 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 @@ -99,7 +99,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) 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) @@ -134,14 +134,14 @@ module subroutine fraggle_util_set_mass_dist(self, param) beta = 2.85_DP ! From Leinhardt & Stewart (2012) Mslr = impactors%mass_dist(iMslr) - nfragmax = ceiling(NFRAGMAX_UNSCALED / param%nfrag_reduction) - do i = 1, nfragmax + do i = 1, NFRAGMAX_UNSCALED mfrag = (nfrag)**(-3._DP / BETA) * Mslr mfrag = max(mfrag, min_mfrag) mremaining = mremaining - mfrag nfrag = nfrag + 1 if (mremaining < 0.0_DP) exit end do + nfragmax = ceiling(NFRAGMAX_UNSCALED / param%nfrag_reduction) nfrag = max(min(ceiling(nfrag / param%nfrag_reduction), nfragmax), NFRAGMIN) call self%setup_fragments(nfrag) @@ -177,16 +177,42 @@ module subroutine fraggle_util_set_mass_dist(self, param) ! Use Newton's method solver to get the logspace slope of the mass function Mrat = mremaining / min_mfrag nrem = nfrag - 2 - mscale = Mrat**(1.0_DP/nrem) - 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) - mscale0 = mscale - mscale = mscale0 + y/yp - if (abs((mscale - mscale0)/mscale) < epsilon(1.0_DP)) exit + x0 = 1.0_DP + 100*epsilon(1.0_DP) + x1 = Mrat**(1.0/nrem) + do j = 1, MAXLOOP + y0 = Mrat - (1.0_DP - x0**nrem)/(1.0_DP - x0) + y1 = Mrat - (1.0_DP - x1**nrem)/(1.0_DP - x1) + if (y0*y1 < 0.0_DP) exit + x1 = x1 * 1.6_DP end do - Mslr = impactors%mass_dist(iMslr) + ! Find the mass scaling factor with bisection + do i = 1, MAXLOOP + mscale= 0.5_DP * (x0 + x1) + ymid = Mrat - (1.0_DP - mscale**nrem)/(1.0_DP - mscale) + if (abs((y0 - ymid)/y0) < 1e-12_DP) exit + if (y0 * ymid < 0.0_DP) then + x1 = mscale + y1 = ymid + else if (y1 * ymid < 0.0_DP) then + x0 = mscale + y0 = ymid + else + exit + end if + ! Finish with Newton if we still haven't made it + if (i == MAXLOOP) then + 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) + eps = y/yp + if (abs(eps) < epsilon(1.0_DP) * mscale) exit + mscale = mscale + eps + end do + end if + end do + + Mslr = impactors%mass_dist(iMslr) mass(iMslr) = Mslr mfrag = min_mfrag do i = iMrem,nfrag From 8df230d05fa4dd2693e3352b749050bcfc6f5d67 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Feb 2023 18:18:05 -0500 Subject: [PATCH 24/38] Made position finding more efficient --- src/fraggle/fraggle_generate.f90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 099417734..911065308 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 @@ -310,7 +310,7 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu 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_max = 1.02_DP real(DP), parameter :: pack_density = 0.5236_DP ! packing density of loose spheres rbuffer = 1.01_DP @@ -415,12 +415,14 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu fragments%rmag(:) = .mag. fragments%rc(:,:) ! Check for any overlapping bodies. - loverlap(:) = .false. + !loverlap(:) = .false. do j = 1, nfrag + if (.not.loverlap(j)) cycle + loverlap(j) = .false. ! Check for overlaps between fragments do i = j + 1, nfrag dis = .mag.(fragments%rc(:,j) - fragments%rc(:,i)) - loverlap(i) = loverlap(i) .or. (dis <= rbuffer * (fragments%radius(i) + fragments%radius(j))) + loverlap(j) = loverlap(j) .or. (dis <= rbuffer * (fragments%radius(i) + fragments%radius(j))) end do ! Check for overlaps with existing bodies that are not involved in the collision do i = 1, npl From 13f6a8a62a4dbe8155f3bf55112d24f11de88fc9 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Feb 2023 18:28:13 -0500 Subject: [PATCH 25/38] Fixed problem that was causing bodies to overlap --- src/fraggle/fraggle_generate.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 911065308..5ef782c0a 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -420,7 +420,7 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu if (.not.loverlap(j)) cycle loverlap(j) = .false. ! Check for overlaps between fragments - do i = j + 1, nfrag + do concurrent(i = 1:nfrag, i/=j) dis = .mag.(fragments%rc(:,j) - fragments%rc(:,i)) loverlap(j) = loverlap(j) .or. (dis <= rbuffer * (fragments%radius(i) + fragments%radius(j))) end do From 2e7f0becf108579bbf7b6f58262272c38fb1d491 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Feb 2023 18:28:58 -0500 Subject: [PATCH 26/38] Added the chunk method in order to use the interpolation on dask-compatible ufuncs --- examples/Fragmentation/Fragmentation_Movie.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index b60dbc1db..adb37d371 100755 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -166,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").chunk(dict(time=-1)).interpolate_na(dim="time", method="akima") # 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] @@ -359,7 +359,7 @@ def vec_props(self, c): minimum_fragment_gmass = 0.01 * body_Gmass[style][1] gmtiny = 0.10 * body_Gmass[style][1] sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, nfrag_reduction=nfrag_reduction[style]) - sim.run(dt=5e-4, tstop=tstop[style], istep_out=1, dump_cadence=0) + sim.run(dt=5e-4, tstop=tstop[style], istep_out=1, dump_cadence=0, dask=False) print("Generating animation") anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1) From 927e13773edf0f03689bb81d0d8dfe48b97aec9a Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Feb 2023 19:56:16 -0500 Subject: [PATCH 27/38] Added dask support to Python side --- examples/Chambers2013/init_cond.py | 2 +- examples/Chambers2013/scattermovie.py | 2 +- python/swiftest/swiftest/io.py | 10 ++- python/swiftest/swiftest/simulation_class.py | 71 ++++++++++++-------- python/swiftest/swiftest/tool.py | 2 +- 5 files changed, 54 insertions(+), 33 deletions(-) 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/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index d070bbb91..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,8 +844,10 @@ 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) - ds = xr.open_mfdataset(param['BIN_OUT'], parallel=True, engine='h5netcdf', 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: 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 225619204..285c74829 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}, dask="parallelized" + np.linalg.norm, ds[x].where(~np.isnan(ds[x])), input_core_dims=[[dim]], kwargs={"ord": ord, "axis": -1}, dask="parallelized" ) def wrap_angle(angle): From 110d3343282abd59911b2170ffde5b9293dc7c78 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Feb 2023 20:03:33 -0500 Subject: [PATCH 28/38] Updated the SFD code to pin the top of the SFD rather than the bottom. This means that bodies smaller than the minimum can be made, but the minimum still sets the total number --- examples/Fragmentation/Fragmentation_Movie.py | 4 ++-- src/fraggle/fraggle_util.f90 | 9 +++++---- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index adb37d371..b60dbc1db 100755 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -166,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").chunk(dict(time=-1)).interpolate_na(dim="time", method="akima") + ds = xr.combine_nested([data,enc],concat_dim='time').sortby("time").interpolate_na(dim="time", method="akima") # 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] @@ -359,7 +359,7 @@ def vec_props(self, c): minimum_fragment_gmass = 0.01 * body_Gmass[style][1] gmtiny = 0.10 * body_Gmass[style][1] sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, nfrag_reduction=nfrag_reduction[style]) - sim.run(dt=5e-4, tstop=tstop[style], istep_out=1, dump_cadence=0, dask=False) + sim.run(dt=5e-4, tstop=tstop[style], istep_out=1, dump_cadence=0) print("Generating animation") anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index d1327d723..92bbb4aba 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -175,15 +175,16 @@ module subroutine fraggle_util_set_mass_dist(self, param) ! The mass will be distributed evenly in logspace between the second-largest remnant and the minimum mass ! Use Newton's method solver to get the logspace slope of the mass function - Mrat = mremaining / min_mfrag + Mrat = (mremaining + Mslr) / Mslr nrem = nfrag - 2 - x0 = 1.0_DP + 100*epsilon(1.0_DP) + x0 = 1.0_DP - 100*epsilon(1.0_DP) x1 = Mrat**(1.0/nrem) do j = 1, MAXLOOP y0 = Mrat - (1.0_DP - x0**nrem)/(1.0_DP - x0) y1 = Mrat - (1.0_DP - x1**nrem)/(1.0_DP - x1) 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 @@ -214,11 +215,11 @@ module subroutine fraggle_util_set_mass_dist(self, param) Mslr = impactors%mass_dist(iMslr) mass(iMslr) = Mslr - mfrag = min_mfrag + mfrag = Mslr do i = iMrem,nfrag + mfrag = mfrag * mscale mass(i) = mfrag mremaining = mremaining - mfrag - mfrag = max(mfrag * mscale, min_mfrag) 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. From d3b9a93d094e35660d1ab70492c6adc7d854161a Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Feb 2023 20:44:10 -0500 Subject: [PATCH 29/38] Fixed issues converging on position and fragment number --- src/fraggle/fraggle_generate.f90 | 16 +++++++--------- src/fraggle/fraggle_util.f90 | 2 +- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 5ef782c0a..554da1d31 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -309,10 +309,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.02_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) @@ -415,16 +413,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 concurrent(i = 1:nfrag, i/=j) + do i = 1,nfrag + if (i == j) cycle dis = .mag.(fragments%rc(:,j) - fragments%rc(:,i)) - loverlap(j) = loverlap(j) .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(:))) @@ -432,7 +431,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(:)) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 92bbb4aba..5863a3073 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -177,7 +177,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) ! Use Newton's method solver to get the logspace slope of the mass function Mrat = (mremaining + Mslr) / Mslr nrem = nfrag - 2 - x0 = 1.0_DP - 100*epsilon(1.0_DP) + x0 = 1.0_DP - 1.0_DP/Mrat x1 = Mrat**(1.0/nrem) do j = 1, MAXLOOP y0 = Mrat - (1.0_DP - x0**nrem)/(1.0_DP - x0) From 5264393e4d9276ec951ce78e91ce5b396b2300f2 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Feb 2023 21:02:01 -0500 Subject: [PATCH 30/38] Fixed issues arising when the mass contrast in fragments gets too big. --- src/collision/collision_util.f90 | 4 ++-- src/fraggle/fraggle_generate.f90 | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index bdbf35fb9..19a688c8a 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -833,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 554da1d31..5dac7024f 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -682,7 +682,7 @@ 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) From 0e00217ffb3b6c241c45b04b540614ce3346cf94 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 21 Feb 2023 13:07:26 -0500 Subject: [PATCH 31/38] Changed the fragmentation model to use a power law SFD for fragments. min_mfrag is now just a reference mass to set the number of fragments rather than a hard lower limit. --- src/fraggle/fraggle_util.f90 | 75 +++++++++++++++++------------------- 1 file changed, 36 insertions(+), 39 deletions(-) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 5863a3073..0d1a0604d 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -86,14 +86,14 @@ 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, x0, x1, xmid, ymid, y0, y1, y, yp, eps, Mrat + 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 @@ -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) - - do i = 1, NFRAGMAX_UNSCALED - mfrag = (nfrag)**(-3._DP / BETA) * Mslr - mfrag = max(mfrag, min_mfrag) - mremaining = mremaining - mfrag - nfrag = nfrag + 1 - if (mremaining < 0.0_DP) exit - 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,53 +163,60 @@ 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 + Mslr) / Mslr - nrem = nfrag - 2 - x0 = 1.0_DP - 1.0_DP/Mrat - x1 = Mrat**(1.0/nrem) + nrem = nfrag - iMslr + 1 + x0 = 0.1_DP + x1 = 5.0_DP do j = 1, MAXLOOP - y0 = Mrat - (1.0_DP - x0**nrem)/(1.0_DP - x0) - y1 = Mrat - (1.0_DP - x1**nrem)/(1.0_DP - x1) + y0 = Mrat - 1.0_DP + y1 = Mrat - 1.0_DP + do i = 2, nrem + y0 = y0 - (i)**(-x0/3.0_DP) + y1 = y1 - (i)**(-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 i = 1, MAXLOOP - mscale= 0.5_DP * (x0 + x1) - ymid = Mrat - (1.0_DP - mscale**nrem)/(1.0_DP - mscale) + do j = 1, MAXLOOP + beta = 0.5_DP * (x0 + x1) + ymid = Mrat - 1.0_DP + do i = 2, nrem + ymid = ymid - (i)**(-beta/3.0_DP) + end do if (abs((y0 - ymid)/y0) < 1e-12_DP) exit if (y0 * ymid < 0.0_DP) then - x1 = mscale + x1 = beta y1 = ymid else if (y1 * ymid < 0.0_DP) then - x0 = mscale + x0 = beta y0 = ymid else exit end if ! Finish with Newton if we still haven't made it - if (i == MAXLOOP) then - 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) + if (j == MAXLOOP) then + do k = 1, MAXLOOP + y = Mrat - 1.0_DP + yp = 0.0_DP + do i = 2, nrem + y = y - (i)**(-beta/3.0_DP) + if (i > 3) yp = yp - (i)**(-beta/3.0_DP) * log(real(i,kind=DP)) + end do eps = y/yp - if (abs(eps) < epsilon(1.0_DP) * mscale) exit - mscale = mscale + eps + 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 = Mslr do i = iMrem,nfrag - mfrag = mfrag * mscale - mass(i) = mfrag - mremaining = mremaining - 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. From 5cd7c93b3ecf0a950e187a69c36f270c8b62fbd5 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 21 Feb 2023 15:21:12 -0500 Subject: [PATCH 32/38] Fixed bugs in the fragment SFD calculation --- src/fraggle/fraggle_util.f90 | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 0d1a0604d..8a7a795f2 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -166,15 +166,14 @@ module subroutine fraggle_util_set_mass_dist(self, param) ! 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 + Mslr) / Mslr - nrem = nfrag - iMslr + 1 x0 = 0.1_DP x1 = 5.0_DP do j = 1, MAXLOOP y0 = Mrat - 1.0_DP y1 = Mrat - 1.0_DP - do i = 2, nrem - y0 = y0 - (i)**(-x0/3.0_DP) - y1 = y1 - (i)**(-x1/3.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 @@ -185,8 +184,8 @@ module subroutine fraggle_util_set_mass_dist(self, param) do j = 1, MAXLOOP beta = 0.5_DP * (x0 + x1) ymid = Mrat - 1.0_DP - do i = 2, nrem - ymid = ymid - (i)**(-beta/3.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 @@ -203,9 +202,9 @@ module subroutine fraggle_util_set_mass_dist(self, param) do k = 1, MAXLOOP y = Mrat - 1.0_DP yp = 0.0_DP - do i = 2, nrem - y = y - (i)**(-beta/3.0_DP) - if (i > 3) yp = yp - (i)**(-beta/3.0_DP) * log(real(i,kind=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 @@ -214,7 +213,6 @@ module subroutine fraggle_util_set_mass_dist(self, param) end if end do - Mslr = impactors%mass_dist(iMslr) do i = iMrem,nfrag mass(i) = Mslr * (i-1)**(-beta/3.0_DP) end do From abc674adcd730e1503a9bc825ac9ba7e9a816892 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 21 Feb 2023 17:39:47 -0500 Subject: [PATCH 33/38] Changes to the positining that mainly affects hit and runs to increase the success rate of fragmentations after hit and run. --- src/fraggle/fraggle_generate.f90 | 53 +++++++++++++++++--------------- 1 file changed, 28 insertions(+), 25 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 5dac7024f..8dec5ae39 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -326,20 +326,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 @@ -353,9 +350,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 @@ -363,6 +359,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 @@ -370,12 +368,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)) @@ -383,11 +381,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 + 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) + (fragment_cloud_radius(j) + maxval(mass_rscale)) * 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(:)) @@ -549,7 +552,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu real(DP) :: vmax_guess integer(I4B), parameter :: MAXLOOP = 25 integer(I4B), parameter :: MAXTRY = 10 - integer(I4B), parameter :: MAXANGMTM = 1000 + integer(I4B), parameter :: MAXANGMTM = 30000 class(collision_fraggle), allocatable :: collider_local character(len=STRMAX) :: message @@ -682,11 +685,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%mass(i) / sum(fragments%mass(istart:fragments%nbody)) + 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) From 49c9a720bc73be24ada5278117475a38bd0083af Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 21 Feb 2023 17:41:47 -0500 Subject: [PATCH 34/38] Increased stretching factor for hit and run fragment cloud --- src/fraggle/fraggle_generate.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 8dec5ae39..0b31b98a8 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -384,7 +384,7 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu ! Stretch out the hit and run cloud along the flight trajectory if (lhitandrun) then - fragments%rc(:,i) = fragments%rc(:,i) * (1.0_DP + fragment_cloud_radius(j) * mass_rscale(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) From 7845d71b9fbc98c7855a28afce1f8b324169e545 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 21 Feb 2023 17:50:51 -0500 Subject: [PATCH 35/38] More position tweaks to hit and run --- src/fraggle/fraggle_generate.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 0b31b98a8..6dcf6ab0e 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -390,7 +390,7 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu fragments%rc(:,i) = fragments%rc(:,i) + fragment_cloud_center(:,j) if (lhitandrun) then - fragments%rc(:,i) = fragments%rc(:,i) + (fragment_cloud_radius(j) + maxval(mass_rscale)) * impactors%bounce_unit(:) ! Shift the stretched cloud downrange + 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(:)) From 03b790f5af4078900a4d3dbc1ce157229c2de7e8 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 21 Feb 2023 17:56:23 -0500 Subject: [PATCH 36/38] Adjusted the size of the angular momentum loop --- src/fraggle/fraggle_generate.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 6dcf6ab0e..a8d85b442 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -552,7 +552,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu real(DP) :: vmax_guess integer(I4B), parameter :: MAXLOOP = 25 integer(I4B), parameter :: MAXTRY = 10 - integer(I4B), parameter :: MAXANGMTM = 30000 + integer(I4B), parameter :: MAXANGMTM = 1000 class(collision_fraggle), allocatable :: collider_local character(len=STRMAX) :: message From aac25ad2acbca4f0cc8a27246cc1f136c1220b3f Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 22 Feb 2023 14:15:19 -0500 Subject: [PATCH 37/38] Updated the SFD to allow for negative slopes (when fragment numbers get small) --- examples/Fragmentation/Fragmentation_Movie.py | 2 +- python/swiftest/swiftest/tool.py | 2 +- src/fraggle/fraggle_generate.f90 | 46 +++++++++++++++++++ src/fraggle/fraggle_module.f90 | 17 +++++-- src/fraggle/fraggle_util.f90 | 4 +- 5 files changed, 63 insertions(+), 8 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index b60dbc1db..c44fe0b56 100755 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -166,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/python/swiftest/swiftest/tool.py b/python/swiftest/swiftest/tool.py index 285c74829..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].where(~np.isnan(ds[x])), input_core_dims=[[dim]], kwargs={"ord": ord, "axis": -1}, dask="parallelized" + 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/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index a8d85b442..bbca213f6 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -160,6 +160,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur call fraggle_generate_rot_vec(self, nbody_system, param) call fraggle_generate_vel_vec(self, nbody_system, param, lfailure) end if + lfailure = .true. if (.not.lfailure) then if (self%fragments%nbody /= nfrag_start) then @@ -284,6 +285,51 @@ 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 + 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 !! 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 8a7a795f2..4daef8d40 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -166,8 +166,8 @@ module subroutine fraggle_util_set_mass_dist(self, param) ! 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 + Mslr) / Mslr - x0 = 0.1_DP - x1 = 5.0_DP + x0 = -3.0_DP + x1 = 3.0_DP do j = 1, MAXLOOP y0 = Mrat - 1.0_DP y1 = Mrat - 1.0_DP From 2213b911bbe6bae5a6d6eddd33c5d7b7010aec66 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 22 Feb 2023 14:22:42 -0500 Subject: [PATCH 38/38] Added ability to convert over to hit and run if merger body spins too fast. --- src/fraggle/fraggle_generate.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index bbca213f6..6ee9c36bd 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -160,7 +160,6 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur call fraggle_generate_rot_vec(self, nbody_system, param) call fraggle_generate_vel_vec(self, nbody_system, param, lfailure) end if - lfailure = .true. if (.not.lfailure) then if (self%fragments%nbody /= nfrag_start) then @@ -320,6 +319,7 @@ module subroutine fraggle_generate_merge(self, nbody_system, param, t) 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