diff --git a/examples/Basic_Simulation/cb.in b/examples/Basic_Simulation/cb.in deleted file mode 100644 index b2cb85c35..000000000 --- a/examples/Basic_Simulation/cb.in +++ /dev/null @@ -1,7 +0,0 @@ -Sun -39.476926408897626 -0.004650467260962157 -4.7535806948127355e-12 --2.2473967953572827e-18 -0.0 0.0 0.07 -11.209306302144773 -38.759372036774764 82.25088158389266 diff --git a/examples/Basic_Simulation/initial_conditions.py b/examples/Basic_Simulation/initial_conditions.py index 7e0218a03..861556958 100644 --- a/examples/Basic_Simulation/initial_conditions.py +++ b/examples/Basic_Simulation/initial_conditions.py @@ -11,24 +11,37 @@ #!/usr/bin/env python3 """ -Generates a set of Swiftest input files from initial conditions. +Generates and runs a set of Swiftest input files from initial conditions with the SyMBA integrator. All simulation +outputs are stored in the /simdata subdirectory. -Returns -------- -Updates sim.data with the simulation data +Input +------ +None. + +Output +------ +bin.nc : A NetCDF file containing the simulation output. +dump_bin1.nc : A NetCDF file containing the necessary inputs to restart a simulation from t!=0. +dump_bin2.nc : A NetCDF file containing the necessary inputs to restart a simulation from t!=0. +dump_param1.in : An ASCII file containing the necessary parameters to restart a simulation. +dump_param2.in : An ASCII file containing the necessary parameters to restart a simulation. +fraggle.log : An ASCII file containing the information of any collisional events that occured. +init_cond.nc : A NetCDF file containing the initial conditions for the simulation. +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 -# Initialize the simulation object as a variable +# Initialize the simulation object as a variable. Arguments may be defined here or through the sim.run() method. sim = swiftest.Simulation(fragmentation=True, minimum_fragment_mass = 2.5e-11, mtiny=2.5e-8) -# Add the modern planets and the Sun using the JPL Horizons Database +# Add the modern planets and the Sun using the JPL Horizons Database. sim.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune","Pluto"]) -# Add 5 user-defined massive bodies +# Add 5 user-defined massive bodies. npl = 5 density_pl = 3000.0 / (sim.param['MU2KG'] / sim.param['DU2M'] ** 3) @@ -47,7 +60,7 @@ sim.add_body(name=name_pl, a=a_pl, e=e_pl, inc=inc_pl, capom=capom_pl, omega=omega_pl, capm=capm_pl, Gmass=GM_pl, radius=R_pl, rhill=Rh_pl, Ip=Ip_pl, rot=rot_pl) -# Add 10 user-defined test particles +# Add 10 user-defined test particles. ntp = 10 name_tp = ["TestParticle_01", "TestParticle_02", "TestParticle_03", "TestParticle_04", "TestParticle_05", "TestParticle_06", "TestParticle_07", "TestParticle_08", "TestParticle_09", "TestParticle_10"] @@ -59,8 +72,8 @@ capm_tp = default_rng().uniform(0.0, 360.0, ntp) sim.add_body(name=name_tp, a=a_tp, e=e_tp, inc=inc_tp, capom=capom_tp, omega=omega_tp, capm=capm_tp) -# Display the run configuration parameters +# Display the run configuration parameters. sim.get_parameter() -# Run the simulation +# Run the simulation. Arguments may be defined here or thorugh the swiftest.Simulation() method. sim.run(tstart=0.0, tstop=1.0e3, dt=0.01, tstep_out=1.0e0, dump_cadence=0) diff --git a/examples/Basic_Simulation/output_reader.py b/examples/Basic_Simulation/output_reader.py index b171eb402..a41103ccd 100644 --- a/examples/Basic_Simulation/output_reader.py +++ b/examples/Basic_Simulation/output_reader.py @@ -14,33 +14,32 @@ Reads and processes a Swiftest output file. Input -------- -out.nc : NetCDF file - Swiftest output file. +------ +bin.nc : A NetCDF file containing the simulation output. -Returns -------- -output.eps : Encapsulated PostScript file. - A figure containing the eccentricity vs. semi-major axis for all bodies at the start of the simulation. +Output +------ +output.eps : Encapsulated PostScript file depicting the eccentricity vs. semi-major axis for all bodies at the start + of the simulation. """ import swiftest import xarray as xr import matplotlib.pyplot as plt -# Read in the simulation output and store it as an Xarray dataset -ds = swiftest.Simulation(param_file="simdata/param.in", read_old_output_file=True).data +# Read in the simulation output and store it as an Xarray dataset. +sim = swiftest.Simulation(read_old_output_file=True) -# Plot of the data and save the output plot -colors = ['white' if x == 'Massive Body' else 'black' for x in ds['particle_type']] -sizes = [100 if x == 'Massive Body' else 10 for x in ds['particle_type']] +# Plot of the data and save the output plot. +colors = ['white' if x == 'Massive Body' else 'black' for x in sim.data['particle_type']] +sizes = [100 if x == 'Massive Body' else 10 for x in sim.data['particle_type']] fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5,3)) -ax.set(xlabel="Semimajor Axis (AU)", ylabel="Eccentricity", title="Simulation Start") -ax.scatter(ds['a'].isel(time=0), ds['e'].isel(time=0), c=colors, s=sizes, edgecolor='black') +ax.set(xlabel="Semimajor Axis (AU)", ylabel="Eccentricity", title="Simulation Initial Conditions (t=0)") +ax.scatter(sim.data['a'].isel(time=0), sim.data['e'].isel(time=0), c=colors, s=sizes, edgecolor='black') ax.set_xlim(0, 2.0) ax.set_ylim(0, 0.4) -ax.text(1.5, 0.35, f"t = {ds['time'].isel(time=0).values} years", size=10, ha="left") +ax.text(1.5, 0.35, f"t = {sim.data['time'].isel(time=0).values} years", size=10, ha="left") plt.tight_layout() plt.show() -fig.savefig("output.eps", dpi=300, facecolor='white', transparent=False, bbox_inches="tight") \ No newline at end of file +fig.savefig("output.eps", dpi=300, facecolor='white', transparent=False, bbox_inches="tight") diff --git a/examples/Basic_Simulation/pl.in b/examples/Basic_Simulation/pl.in deleted file mode 100644 index 215b91f1d..000000000 --- a/examples/Basic_Simulation/pl.in +++ /dev/null @@ -1,85 +0,0 @@ -14 -Mercury 6.553709809565314e-06 -1.6306381826061646e-05 -0.3870985843095394 0.2056234010897001 7.003302508001384 -48.29611837378607 29.20442403952454 338.3394874682879 -0.0 0.0 0.346 -3.573018077015318 -18.380317085416586 34.36143850429876 -Venus 9.663313399581537e-05 -4.0453784346544176e-05 -0.7233297579736101 0.006717605698865438 3.394439273342282 -76.60235891771119 54.96037946082961 200.4789339550648 -0.0 0.0 0.4 -0.17650282045605922 -3.6612475825356214 8.702866268072764 -Earth 0.00012002693582795245 -4.25875607065041e-05 -0.9999904874543223 0.01671400376545858 0.003637862608863003 -175.025172600231 287.9619628812575 114.3482934042427 -0.0 0.0 0.3307 -6.157239621449141 0.057224133705898524 2301.2082528350797 -Mars 1.2739802010675942e-05 -2.2657408050928896e-05 -1.523711925589535 0.09344151133508208 1.847441673557901 -49.4728572124747 286.7379771285891 209.3396773477138 -0.0 0.0 0.3644 -997.9224351226384 -909.5549030011778 1783.3823046046184 -Jupiter 0.037692251088985676 -0.0004673261703049093 -5.2027278008516 0.04824497711637968 1.303631134570075 -100.5192588433081 273.5898402882514 129.5536700659942 -0.0 0.0 0.2756 --80.9864396731672 -2388.0246092955053 5008.722318533006 -Saturn 0.011285899820091273 -0.00038925687730393614 -9.532011952667288 0.05486329870433341 2.487906363280301 -113.6305781676206 339.5467356402391 290.8995806568904 -0.0 0.0 0.22 -441.95954822014636 378.52638822638795 5135.909115928892 -Uranus 0.001723658947826773 -0.00016953449859497232 -19.24498838290236 0.04796174942301296 0.7730102596086205 -74.0125809801658 93.59554912280227 262.8658637277515 -0.0 0.0 0.23 --677.3000258209181 -3008.1099071905787 -836.3013266185699 -Neptune 0.0020336100526728304 -0.00016458790412449367 -30.03895991152209 0.008955570138096731 1.771119354296142 -131.8221159748827 284.4748429674216 308.4513720536233 -0.0 0.0 0.23 -1232.224106980634 -2177.3040821077648 2329.8227878119233 -Pluto 2.924216771029454e-07 -7.943294877391593e-06 -39.36791814672583 0.2487178537481577 17.1705505990969 -110.3314332962701 113.0826635900664 55.11416408345664 -0.0 0.0 0.4 --243.59404988903637 261.28663002814227 -38.57352022187049 -MassiveBody_01 1.1912109366578089e-05 -2.425055692051245e-05 -0.7452368716298337 0.0633011418780484 0.11151363780595558 -203.01823417718037 284.9353898127118 266.79344592519305 -0.4 0.4 0.4 -0.0 0.0 0.0 -MassiveBody_02 1.5882812488770783e-05 -2.6691191565570074e-05 -0.7671203280602826 0.10149029388964753 53.46706656938751 -61.74738152068808 68.20565593722856 271.3352706902475 -0.4 0.4 0.4 -0.0 0.0 0.0 -MassiveBody_03 1.985351561096348e-05 -2.8752214513575297e-05 -1.4698824276418418 0.13621250684495437 23.635498327264845 -38.071905339231236 283.134612455057 250.67457601578352 -0.4 0.4 0.4 -0.0 0.0 0.0 -MassiveBody_04 5.9560546832890444e-05 -4.14678690275904e-05 -0.9741731590760117 0.1519713326893784 20.51335588582416 -350.53780805624825 304.05941264938997 142.62713592644738 -0.4 0.4 0.4 -0.0 0.0 0.0 -MassiveBody_05 9.926757805481742e-05 -4.916559523190238e-05 -0.6633500122731075 0.13550930562584215 10.680323453316653 -328.45970148163457 49.93948991697533 316.2109831817007 -0.4 0.4 0.4 -0.0 0.0 0.0 diff --git a/examples/Basic_Simulation/read_old_run.ipynb b/examples/Basic_Simulation/read_old_run.ipynb deleted file mode 100644 index 0a7459084..000000000 --- a/examples/Basic_Simulation/read_old_run.ipynb +++ /dev/null @@ -1,853 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "id": "86b309f3-d400-4164-88fc-cf56b6cfcf5f", - "metadata": {}, - "outputs": [], - "source": [ - "import swiftest" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "353b2559-83a2-496f-8648-60f6121333a5", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Reading Swiftest file /home/daminton/git_debug/swiftest/examples/Basic_Simulation/simdata/param.in\n", - "\n", - "Creating Dataset from NetCDF file\n", - "Successfully converted 78 output frames.\n", - "Swiftest simulation data stored as xarray DataSet .data\n" - ] - } - ], - "source": [ - "sim = swiftest.Simulation(read_old_output_file=True)" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "id": "c7bceca2-678d-469d-a3a7-06cfa9f29fc0", - "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "
<xarray.Dataset>\n",
-       "Dimensions:          (time: 78, name: 25)\n",
-       "Coordinates:\n",
-       "  * time             (time) float64 0.0 2.0 nan 4.0 nan ... nan 76.0 nan 78.0\n",
-       "  * name             (name) <U32 'Sun' 'Mercury' ... 'TestParticle_10'\n",
-       "Data variables: (12/45)\n",
-       "    id               (name) int64 0 1 2 3 4 5 6 7 8 ... 17 18 19 20 21 22 23 24\n",
-       "    npl              (time) int64 14 14 -2147483647 14 ... 14 -2147483647 14\n",
-       "    ntp              (time) int64 10 10 -2147483647 10 ... 10 -2147483647 10\n",
-       "    nplm             (time) int64 13 13 -2147483647 13 ... 13 -2147483647 13\n",
-       "    particle_type    (name) <U32 'Central Body' ... 'Test Particle'\n",
-       "    status           (name) <U32 'ACTIVE' 'ACTIVE' ... 'ACTIVE' 'ACTIVE'\n",
-       "    ...               ...\n",
-       "    Ip3              (time, name) float64 0.07 0.346 0.4 0.3307 ... nan nan nan\n",
-       "    rotx             (time, name) float64 11.21 3.573 0.1765 ... nan nan nan\n",
-       "    roty             (time, name) float64 -38.76 -18.38 -3.661 ... nan nan nan\n",
-       "    rotz             (time, name) float64 82.25 34.36 8.703 ... nan nan nan\n",
-       "    j2rp2            (time) float64 4.754e-12 4.754e-12 nan ... nan 4.754e-12\n",
-       "    j4rp4            (time) float64 -2.247e-18 -2.247e-18 nan ... nan -2.247e-18
" - ], - "text/plain": [ - "\n", - "Dimensions: (time: 78, name: 25)\n", - "Coordinates:\n", - " * time (time) float64 0.0 2.0 nan 4.0 nan ... nan 76.0 nan 78.0\n", - " * name (name) system%pl) - class is (symba_pl) - pl%status(colliders%idx(:)) = status - pl%ldiscard(colliders%idx(:)) = .false. - pl%lcollision(colliders%idx(:)) = .false. - end select - else - ! Populate the list of new bodies - nfrag = frag%nbody - write(message, *) nfrag - call io_log_one_message(FRAGGLE_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") - select case(frag%regime) - case(COLLRESOLVE_REGIME_DISRUPTION) - status = DISRUPTION - case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - status = SUPERCATASTROPHIC - end select - frag%id(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)] - param%maxid = frag%id(nfrag) - call symba_collision_mergeaddsub(system, param, colliders, frag, status) - end if - - return - end function symba_collision_casedisruption - - - module function symba_collision_casehitandrun(system, param, colliders, frag) result(status) - !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Create the fragments resulting from a non-catastrophic hit-and-run collision - !! - implicit none - ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object - class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object - ! Result - integer(I4B) :: status !! Status flag assigned to this outcom - ! Internals - integer(I4B) :: i, ibiggest, nfrag, jtarg, jproj - logical :: lpure - character(len=STRMAX) :: message - - message = "Hit and run between" - call symba_collision_collider_message(system%pl, colliders%idx, message) - call io_log_one_message(FRAGGLE_LOG_OUT, trim(adjustl(message))) - - if (colliders%mass(1) > colliders%mass(2)) then - jtarg = 1 - jproj = 2 - else - jtarg = 2 - jproj = 1 - end if - - if (frag%mass_dist(2) > 0.9_DP * colliders%mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched - call io_log_one_message(FRAGGLE_LOG_OUT, "Pure hit and run. No new fragments generated.") - nfrag = 0 - lpure = .true. - else ! Imperfect hit and run, so we'll keep the largest body and destroy the other - lpure = .false. - call frag%set_mass_dist(colliders, param) - - ! Generate the position and velocity distributions of the fragments - call frag%generate_fragments(colliders, system, param, lpure) - - if (lpure) then - call io_log_one_message(FRAGGLE_LOG_OUT, "Should have been a pure hit and run instead") - nfrag = 0 - else - nfrag = frag%nbody - write(message, *) nfrag - call io_log_one_message(FRAGGLE_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") - end if - end if - if (lpure) then ! Reset these bodies back to being active so that nothing further is done to them - status = HIT_AND_RUN_PURE - select type(pl => system%pl) - class is (symba_pl) - pl%status(colliders%idx(:)) = ACTIVE - pl%ldiscard(colliders%idx(:)) = .false. - pl%lcollision(colliders%idx(:)) = .false. - end select - else - ibiggest = colliders%idx(maxloc(system%pl%Gmass(colliders%idx(:)), dim=1)) - frag%id(1) = system%pl%id(ibiggest) - frag%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] - param%maxid = frag%id(nfrag) - status = HIT_AND_RUN_DISRUPT - call symba_collision_mergeaddsub(system, param, colliders, frag, status) - end if - - return - end function symba_collision_casehitandrun - - - module function symba_collision_casemerge(system, param, colliders, frag) result(status) - !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Merge massive bodies. - !! - !! 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 - ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object - class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object - ! Result - integer(I4B) :: status !! Status flag assigned to this outcome - ! Internals - integer(I4B) :: i, j, k, ibiggest - real(DP) :: pe - real(DP), dimension(NDIM) :: L_spin_new - character(len=STRMAX) :: message - - message = "Merging" - call symba_collision_collider_message(system%pl, colliders%idx, message) - call io_log_one_message(FRAGGLE_LOG_OUT, message) - - select type(pl => system%pl) - class is (symba_pl) - - call frag%set_mass_dist(colliders, param) - ibiggest = colliders%idx(maxloc(pl%Gmass(colliders%idx(:)), dim=1)) - frag%id(1) = pl%id(ibiggest) - frag%xb(:,1) = frag%xbcom(:) - frag%vb(:,1) = frag%vbcom(:) - - if (param%lrotation) then - ! Conserve angular momentum by putting pre-impact orbital momentum into spin of the new body - L_spin_new(:) = colliders%L_orbit(:,1) + colliders%L_orbit(:,2) + colliders%L_spin(:,1) + colliders%L_spin(:,2) - - ! Assume prinicpal axis rotation on 3rd Ip axis - frag%rot(:,1) = L_spin_new(:) / (frag%Ip(3,1) * frag%mass(1) * frag%radius(1)**2) - else ! If spin is not enabled, we will consider the lost pre-collision angular momentum as "escaped" and add it to our bookkeeping variable - param%Lescape(:) = param%Lescape(:) + colliders%L_orbit(:,1) + colliders%L_orbit(:,2) - end if - - ! Keep track of the component of potential energy due to the pre-impact colliders%idx for book-keeping - pe = 0.0_DP - do j = 1, colliders%ncoll - do i = j + 1, colliders%ncoll - pe = pe - pl%Gmass(i) * pl%mass(j) / norm2(pl%xb(:, i) - pl%xb(:, j)) - end do - end do - system%Ecollisions = system%Ecollisions + pe - system%Euntracked = system%Euntracked - pe - - ! Update any encounter lists that have the removed bodies in them so that they instead point to the new - do k = 1, system%plplenc_list%nenc - do j = 1, colliders%ncoll - i = colliders%idx(j) - if (i == ibiggest) cycle - if (system%plplenc_list%id1(k) == pl%id(i)) then - system%plplenc_list%id1(k) = pl%id(ibiggest) - system%plplenc_list%index1(k) = i - end if - if (system%plplenc_list%id2(k) == pl%id(i)) then - system%plplenc_list%id2(k) = pl%id(ibiggest) - system%plplenc_list%index2(k) = i - end if - if (system%plplenc_list%id1(k) == system%plplenc_list%id2(k)) system%plplenc_list%status(k) = INACTIVE - end do - end do - - status = MERGED - - call symba_collision_mergeaddsub(system, param, colliders, frag, status) - - end select - - return - end function symba_collision_casemerge - - - subroutine symba_collision_collider_message(pl, collidx, collider_message) - !! author: David A. Minton - !! - !! Prints a nicely formatted message about which bodies collided, including their names and ids. - !! This subroutine appends the body names and ids to an input message. - implicit none - ! Arguments - class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object - integer(I4B), dimension(:), intent(in) :: collidx !! Index of collisional colliders%idx members - character(*), intent(inout) :: collider_message !! The message to print to the screen. - ! Internals - integer(I4B) :: i, n - character(len=STRMAX) :: idstr - - n = size(collidx) - if (n == 0) return - - do i = 1, n - if (i > 1) collider_message = trim(adjustl(collider_message)) // " and " - collider_message = " " // trim(adjustl(collider_message)) // " " // trim(adjustl(pl%info(collidx(i))%name)) - write(idstr, '(I10)') pl%id(collidx(i)) - collider_message = trim(adjustl(collider_message)) // " (" // trim(adjustl(idstr)) // ") " - end do - - return - end subroutine symba_collision_collider_message - - - module function symba_collision_check_encounter(self, system, param, t, dt, irec) result(lany_collision) - !! author: David A. Minton - !! - !! Check for merger between massive bodies and test particles in SyMBA - !! - !! Adapted from David E. Kaufmann's Swifter routine symba_merge.f90 and symba_merge_tp.f90 - !! - !! Adapted from Hal Levison's Swift routine symba5_merge.f - implicit none - ! Arguments - class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter list object - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! current time - real(DP), intent(in) :: dt !! step size - integer(I4B), intent(in) :: irec !! Current recursion level - ! Result - logical :: lany_collision !! Returns true if cany pair of encounters resulted in a collision - ! Internals - logical, dimension(:), allocatable :: lcollision, lmask - real(DP), dimension(NDIM) :: xr, vr - integer(I4B) :: i, j, k, nenc - real(DP) :: rlim, Gmtot - logical :: isplpl - character(len=STRMAX) :: timestr, idstri, idstrj, message - class(symba_encounter), allocatable :: tmp - - lany_collision = .false. - if (self%nenc == 0) return - - select type(self) - class is (symba_plplenc) - isplpl = .true. - class default - isplpl = .false. - end select - - select type(pl => system%pl) - class is (symba_pl) - select type(tp => system%tp) - class is (symba_tp) - nenc = self%nenc - allocate(lmask(nenc)) - lmask(:) = ((self%status(1:nenc) == ACTIVE) .and. (pl%levelg(self%index1(1:nenc)) >= irec)) - if (isplpl) then - lmask(:) = lmask(:) .and. (pl%levelg(self%index2(1:nenc)) >= irec) - else - lmask(:) = lmask(:) .and. (tp%levelg(self%index2(1:nenc)) >= irec) - end if - if (.not.any(lmask(:))) return - - allocate(lcollision(nenc)) - lcollision(:) = .false. - - if (isplpl) then - do concurrent(k = 1:nenc, lmask(k)) - i = self%index1(k) - j = self%index2(k) - xr(:) = pl%rh(:, i) - pl%rh(:, j) - vr(:) = pl%vb(:, i) - pl%vb(:, j) - rlim = pl%radius(i) + pl%radius(j) - Gmtot = pl%Gmass(i) + pl%Gmass(j) - lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), & - Gmtot, rlim, dt, self%lvdotr(k)) - end do - else - do concurrent(k = 1:nenc, lmask(k)) - i = self%index1(k) - j = self%index2(k) - xr(:) = pl%rh(:, i) - tp%rh(:, j) - vr(:) = pl%vb(:, i) - tp%vb(:, j) - lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), & - pl%Gmass(i), pl%radius(i), dt, self%lvdotr(k)) - end do - end if - - lany_collision = any(lcollision(:)) - if (lany_collision) then - call pl%xh2xb(system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary - do k = 1, nenc - i = self%index1(k) - j = self%index2(k) - if (lcollision(k)) self%status(k) = COLLISION - self%t(k) = t - self%x1(:,k) = pl%rh(:,i) + system%cb%xb(:) - self%v1(:,k) = pl%vb(:,i) - if (isplpl) then - self%x2(:,k) = pl%rh(:,j) + system%cb%xb(:) - self%v2(:,k) = pl%vb(:,j) - if (lcollision(k)) then - ! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional colliders%idx - if (pl%lcollision(i) .or. pl%lcollision(j)) call pl%make_colliders([i,j]) - - ! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step - pl%lcollision([i, j]) = .true. - pl%status([i, j]) = COLLISION - call pl%info(i)%set_value(status="COLLISION", discard_time=t, discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i)) - call pl%info(j)%set_value(status="COLLISION", discard_time=t, discard_rh=pl%rh(:,j), discard_vh=pl%vh(:,j)) - end if - else - self%x2(:,k) = tp%rh(:,j) + system%cb%xb(:) - self%v2(:,k) = tp%vb(:,j) - if (lcollision(k)) then - tp%status(j) = DISCARDED_PLR - tp%ldiscard(j) = .true. - write(idstri, *) pl%id(i) - write(idstrj, *) tp%id(j) - write(timestr, *) t - call tp%info(j)%set_value(status="DISCARDED_PLR", discard_time=t, discard_rh=tp%rh(:,j), discard_vh=tp%vh(:,j)) - write(message, *) "Particle " // trim(adjustl(tp%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & - // " collided with massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & - // " at t = " // trim(adjustl(timestr)) - call io_log_one_message(FRAGGLE_LOG_OUT, message) - end if - end if - end do - end if - end select - end select - - ! Extract the pl-pl or pl-tp encounter list and return the pl-pl or pl-tp collision_list - if (lany_collision) then - select type(self) - class is (symba_plplenc) - call self%extract_collisions(system, param) - class default - allocate(tmp, mold=self) - call self%spill(tmp, lcollision, ldestructive=.true.) ! Remove this encounter pair from the encounter list - end select - end if - - return - end function symba_collision_check_encounter - - - pure elemental function symba_collision_check_one(xr, yr, zr, vxr, vyr, vzr, Gmtot, rlim, dt, lvdotr) result(lcollision) - !! author: David A. Minton - !! - !! Check for a merger between a single pair of particles - !! - !! Adapted from David E. Kaufmann's Swifter routines symba_merge_tp.f90 and symba_merge_pl.f90 - !! - !! Adapted from Hal Levison's Swift routine symba5_merge.f - implicit none - ! Arguments - real(DP), intent(in) :: xr, yr, zr !! Relative position vector components - real(DP), intent(in) :: vxr, vyr, vzr !! Relative velocity vector components - real(DP), intent(in) :: Gmtot !! Sum of G*mass of colliding bodies - real(DP), intent(in) :: rlim !! Collision limit - Typically the sum of the radii of colliding bodies - real(DP), intent(in) :: dt !! Step size - logical, intent(in) :: lvdotr !! Logical flag indicating that these two bodies are approaching in the current substep - ! Result - logical :: lcollision !! Logical flag indicating whether these two bodies will collide or not - ! Internals - real(DP) :: r2, rlim2, a, e, q, vdotr, tcr2, dt2 - - r2 = xr**2 + yr**2 + zr**2 - rlim2 = rlim**2 - - if (r2 <= rlim2) then ! checks if bodies are actively colliding in this time step - lcollision = .true. - else ! if they are not actively colliding in this time step, checks if they are going to collide next time step based on velocities and q - lcollision = .false. - vdotr = xr * vxr + yr * vyr + zr * vzr - if (lvdotr .and. (vdotr > 0.0_DP)) then - tcr2 = r2 / (vxr**2 + vyr**2 + vzr**2) - dt2 = dt**2 - if (tcr2 <= dt2) then - call orbel_xv2aeq(Gmtot, xr, yr, zr, vxr, vyr, vzr, a, e, q) - lcollision = (q < rlim) - end if - end if - end if - - return - end function symba_collision_check_one - - - function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, colliders) result(lflag) - !! author: David A. Minton - !! - !! Loops through the pl-pl collision list and groups families together by index. Outputs the indices of all colliders%idx members, - !! and pairs of quantities (x and v vectors, mass, radius, L_spin, and Ip) that can be used to resolve the collisional outcome. - implicit none - ! Arguments - class(symba_pl), intent(inout) :: pl !! SyMBA massive body object - class(symba_cb), intent(inout) :: cb !! SyMBA central body object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B), dimension(2), intent(inout) :: idx_parent !! Index of the two bodies considered the "parents" of the collision - class(fraggle_colliders), intent(out) :: colliders - ! Result - logical :: lflag !! Logical flag indicating whether a colliders%idx was successfully created or not - ! Internals - type collidx_array - integer(I4B), dimension(:), allocatable :: id - integer(I4B), dimension(:), allocatable :: idx - end type collidx_array - type(collidx_array), dimension(2) :: parent_child_index_array - integer(I4B), dimension(2) :: nchild - integer(I4B) :: i, j, ncolliders, idx_child - real(DP), dimension(2) :: volume, density - real(DP) :: mchild, volchild - real(DP), dimension(NDIM) :: xc, vc, xcom, vcom, xchild, vchild, xcrossv - real(DP), dimension(NDIM,2) :: mxc, vcc - - nchild(:) = pl%kin(idx_parent(:))%nchild - ! If all of these bodies share a parent, but this is still a unique collision, move the last child - ! out of the parent's position and make it the secondary body - if (idx_parent(1) == idx_parent(2)) then - if (nchild(1) == 0) then ! There is only one valid body recorded in this pair (this could happen due to restructuring of the kinship relationships, though it should be rare) - lflag = .false. - call pl%reset_kinship([idx_parent(1)]) - return - end if - idx_parent(2) = pl%kin(idx_parent(1))%child(nchild(1)) - nchild(1) = nchild(1) - 1 - nchild(2) = 0 - pl%kin(idx_parent(:))%nchild = nchild(:) - pl%kin(idx_parent(2))%parent = idx_parent(1) - end if - - colliders%mass(:) = pl%mass(idx_parent(:)) ! Note: This is meant to mass, not G*mass, as the collisional regime determination uses mass values that will be converted to Si - colliders%radius(:) = pl%radius(idx_parent(:)) - volume(:) = (4.0_DP / 3.0_DP) * PI * colliders%radius(:)**3 - - ! Group together the ids and indexes of each collisional parent and its children - do j = 1, 2 - allocate(parent_child_index_array(j)%idx(nchild(j)+ 1)) - allocate(parent_child_index_array(j)%id(nchild(j)+ 1)) - associate(idx_arr => parent_child_index_array(j)%idx, & - id_arr => parent_child_index_array(j)%id, & - ncj => nchild(j), & - pl => pl, & - plkinj => pl%kin(idx_parent(j))) - idx_arr(1) = idx_parent(j) - if (ncj > 0) idx_arr(2:ncj + 1) = plkinj%child(1:ncj) - id_arr(:) = pl%id(idx_arr(:)) - end associate - end do - - ! Consolidate the groups of collsional parents with any children they may have into a single "colliders%idx" index array - ncolliders = 2 + sum(nchild(:)) - allocate(colliders%idx(ncolliders)) - colliders%idx = [parent_child_index_array(1)%idx(:),parent_child_index_array(2)%idx(:)] - - colliders%ncoll = count(pl%lcollision(colliders%idx(:))) - colliders%idx = pack(colliders%idx(:), pl%lcollision(colliders%idx(:))) - colliders%L_spin(:,:) = 0.0_DP - colliders%Ip(:,:) = 0.0_DP - - ! Find the barycenter of each body along with its children, if it has any - do j = 1, 2 - colliders%xb(:, j) = pl%rh(:, idx_parent(j)) + cb%xb(:) - colliders%vb(:, j) = pl%vb(:, idx_parent(j)) - ! Assume principal axis rotation about axis corresponding to highest moment of inertia (3rd Ip) - if (param%lrotation) then - colliders%Ip(:, j) = colliders%mass(j) * pl%Ip(:, idx_parent(j)) - colliders%L_spin(:, j) = colliders%Ip(3, j) * colliders%radius(j)**2 * pl%rot(:, idx_parent(j)) - end if - - if (nchild(j) > 0) then - do i = 1, nchild(j) ! Loop over all children and take the mass weighted mean of the properties - idx_child = parent_child_index_array(j)%idx(i + 1) - if (.not. pl%lcollision(idx_child)) cycle - mchild = pl%mass(idx_child) - xchild(:) = pl%rh(:, idx_child) + cb%xb(:) - vchild(:) = pl%vb(:, idx_child) - volchild = (4.0_DP / 3.0_DP) * PI * pl%radius(idx_child)**3 - volume(j) = volume(j) + volchild - ! Get angular momentum of the child-parent pair and add that to the spin - ! Add the child's spin - if (param%lrotation) then - xcom(:) = (colliders%mass(j) * colliders%xb(:,j) + mchild * xchild(:)) / (colliders%mass(j) + mchild) - vcom(:) = (colliders%mass(j) * colliders%vb(:,j) + mchild * vchild(:)) / (colliders%mass(j) + mchild) - xc(:) = colliders%xb(:, j) - xcom(:) - vc(:) = colliders%vb(:, j) - vcom(:) - xcrossv(:) = xc(:) .cross. vc(:) - colliders%L_spin(:, j) = colliders%L_spin(:, j) + colliders%mass(j) * xcrossv(:) - - xc(:) = xchild(:) - xcom(:) - vc(:) = vchild(:) - vcom(:) - xcrossv(:) = xc(:) .cross. vc(:) - colliders%L_spin(:, j) = colliders%L_spin(:, j) + mchild * xcrossv(:) - - colliders%L_spin(:, j) = colliders%L_spin(:, j) + mchild * pl%Ip(3, idx_child) & - * pl%radius(idx_child)**2 & - * pl%rot(:, idx_child) - colliders%Ip(:, j) = colliders%Ip(:, j) + mchild * pl%Ip(:, idx_child) - end if - - ! Merge the child and parent - colliders%mass(j) = colliders%mass(j) + mchild - colliders%xb(:, j) = xcom(:) - colliders%vb(:, j) = vcom(:) - end do - end if - density(j) = colliders%mass(j) / volume(j) - colliders%radius(j) = (3 * volume(j) / (4 * PI))**(1.0_DP / 3.0_DP) - if (param%lrotation) colliders%Ip(:, j) = colliders%Ip(:, j) / colliders%mass(j) - end do - lflag = .true. - - xcom(:) = (colliders%mass(1) * colliders%xb(:, 1) + colliders%mass(2) * colliders%xb(:, 2)) / sum(colliders%mass(:)) - vcom(:) = (colliders%mass(1) * colliders%vb(:, 1) + colliders%mass(2) * colliders%vb(:, 2)) / sum(colliders%mass(:)) - mxc(:, 1) = colliders%mass(1) * (colliders%xb(:, 1) - xcom(:)) - mxc(:, 2) = colliders%mass(2) * (colliders%xb(:, 2) - xcom(:)) - vcc(:, 1) = colliders%vb(:, 1) - vcom(:) - vcc(:, 2) = colliders%vb(:, 2) - vcom(:) - colliders%L_orbit(:,:) = mxc(:,:) .cross. vcc(:,:) - - ! Destroy the kinship relationships for all members of this colliders%idx - call pl%reset_kinship(colliders%idx(:)) - - return - end function symba_collision_consolidate_colliders - - - module subroutine symba_collision_encounter_extract_collisions(self, system, param) - !! author: David A. Minton - !! - !! Processes the pl-pl encounter list remove only those encounters that led to a collision - !! - implicit none - ! Arguments - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - ! Internals - logical, dimension(:), allocatable :: lplpl_collision - logical, dimension(:), allocatable :: lplpl_unique_parent - integer(I4B), dimension(:), pointer :: plparent - integer(I4B), dimension(:), allocatable :: collision_idx, unique_parent_idx - integer(I4B) :: i, index_coll, ncollisions, nunique_parent, nplplenc - - select type (pl => system%pl) - class is (symba_pl) - associate(plplenc_list => self, idx1 => self%index1, idx2 => self%index2, plparent => pl%kin%parent) - nplplenc = plplenc_list%nenc - allocate(lplpl_collision(nplplenc)) - lplpl_collision(:) = plplenc_list%status(1:nplplenc) == COLLISION - if (.not.any(lplpl_collision)) return - ! Collisions have been detected in this step. So we need to determine which of them are between unique bodies. - - ! Get the subset of pl-pl encounters that lead to a collision - ncollisions = count(lplpl_collision(:)) - allocate(collision_idx(ncollisions)) - collision_idx = pack([(i, i=1, nplplenc)], lplpl_collision) - - ! Get the subset of collisions that involve a unique pair of parents - allocate(lplpl_unique_parent(ncollisions)) - - lplpl_unique_parent(:) = plparent(idx1(collision_idx(:))) /= plparent(idx2(collision_idx(:))) - nunique_parent = count(lplpl_unique_parent(:)) - allocate(unique_parent_idx(nunique_parent)) - unique_parent_idx = pack(collision_idx(:), lplpl_unique_parent(:)) - - ! Scrub all pl-pl collisions involving unique pairs of parents, which will remove all duplicates and leave behind - ! all pairs that have themselves as parents but are not part of the unique parent list. This can hapepn in rare cases - ! due to restructuring of parent/child relationships when there are large numbers of multi-body collisions in a single - ! step - lplpl_unique_parent(:) = .true. - do index_coll = 1, ncollisions - associate(ip1 => plparent(idx1(collision_idx(index_coll))), ip2 => plparent(idx2(collision_idx(index_coll)))) - lplpl_unique_parent(:) = .not. ( any(plparent(idx1(unique_parent_idx(:))) == ip1) .or. & - any(plparent(idx2(unique_parent_idx(:))) == ip1) .or. & - any(plparent(idx1(unique_parent_idx(:))) == ip2) .or. & - any(plparent(idx2(unique_parent_idx(:))) == ip2) ) - end associate - end do - - ! Reassemble collision index list to include only those containing the unique pairs of parents, plus all the non-unique pairs that don't - ! contain a parent body on the unique parent list. - ncollisions = nunique_parent + count(lplpl_unique_parent) - collision_idx = [unique_parent_idx(:), pack(collision_idx(:), lplpl_unique_parent(:))] - - ! Create a mask that contains only the pl-pl encounters that did not result in a collision, and then discard them - lplpl_collision(:) = .false. - lplpl_collision(collision_idx(:)) = .true. - call plplenc_list%spill(system%plplcollision_list, lplpl_collision, ldestructive=.true.) ! Extract any encounters that are not collisions from the list. - end associate - end select - - return - end subroutine symba_collision_encounter_extract_collisions - - - module subroutine symba_collision_make_colliders_pl(self, idx) - !! author: Jennifer L.L. Pouplin, Carlisle A. wishard, and David A. Minton - !! - !! When a single body is involved in more than one collision in a single step, it becomes part of a colliders%idx. - !! The largest body involved in a multi-body collision is the "parent" and all bodies that collide with it are its "children," - !! including those that collide with the children. - !! - !! Adapted from David E. Kaufmann's Swifter routine symba_merge_pl.f90 - !! - !! Adapted from Hal Levison's Swift routine symba5_merge.f - implicit none - ! Arguments - class(symba_pl), intent(inout) :: self !! SyMBA massive body object - integer(I4B), dimension(2), intent(in) :: idx !! Array holding the indices of the two bodies involved in the collision - ! Internals - integer(I4B) :: i, j, index_parent, index_child, p1, p2 - integer(I4B) :: nchild_inherit, nchild_orig, nchild_new - integer(I4B), dimension(:), allocatable :: temp - - associate(pl => self) - p1 = pl%kin(idx(1))%parent - p2 = pl%kin(idx(2))%parent - if (p1 == p2) return ! This is a collision between to children of a shared parent. We will ignore it. - - if (pl%mass(p1) > pl%mass(p2)) then - index_parent = p1 - index_child = p2 - else - index_parent = p2 - index_child = p1 - end if - - ! Expand the child array (or create it if necessary) and copy over the previous lists of children - nchild_orig = pl%kin(index_parent)%nchild - nchild_inherit = pl%kin(index_child)%nchild - nchild_new = nchild_orig + nchild_inherit + 1 - allocate(temp(nchild_new)) - - if (nchild_orig > 0) temp(1:nchild_orig) = pl%kin(index_parent)%child(1:nchild_orig) - ! Find out if the child body has any children of its own. The new parent wil inherit these children - if (nchild_inherit > 0) then - temp(nchild_orig+1:nchild_orig+nchild_inherit) = pl%kin(index_child)%child(1:nchild_inherit) - do i = 1, nchild_inherit - j = pl%kin(index_child)%child(i) - ! Set the childrens' parent to the new parent - pl%kin(j)%parent = index_parent - end do - end if - call pl%reset_kinship([index_child]) - ! Add the new child to its parent - pl%kin(index_child)%parent = index_parent - temp(nchild_new) = index_child - ! Save the new child array to the parent - pl%kin(index_parent)%nchild = nchild_new - call move_alloc(from=temp, to=pl%kin(index_parent)%child) - end associate - - return - end subroutine symba_collision_make_colliders_pl - - - subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) - !! author: David A. Minton - !! - !! Fills the pl_discards and pl_adds with removed and added bodies - !! - implicit none - ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object - class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object - integer(I4B), intent(in) :: status !! Status flag to assign to adds - ! Internals - integer(I4B) :: i, ibiggest, ismallest, iother, nstart, nend, ncolliders, nfrag - logical, dimension(system%pl%nbody) :: lmask - class(symba_pl), allocatable :: plnew, plsub - character(*), parameter :: FRAGFMT = '("Newbody",I0.7)' - character(len=NAMELEN) :: newname - - select type(pl => system%pl) - class is (symba_pl) - select type(pl_discards => system%pl_discards) - class is (symba_merger) - associate(info => pl%info, pl_adds => system%pl_adds, cb => system%cb, npl => pl%nbody) - ! Add the colliders%idx bodies to the subtraction list - ncolliders = colliders%ncoll - nfrag = frag%nbody - - param%maxid_collision = max(param%maxid_collision, maxval(system%pl%info(:)%collision_id)) - param%maxid_collision = param%maxid_collision + 1 - - ! Setup new bodies - allocate(plnew, mold=pl) - call plnew%setup(nfrag, param) - ibiggest = colliders%idx(maxloc(pl%Gmass(colliders%idx(:)), dim=1)) - ismallest = colliders%idx(minloc(pl%Gmass(colliders%idx(:)), dim=1)) - - ! Copy over identification, information, and physical properties of the new bodies from the fragment list - plnew%id(1:nfrag) = frag%id(1:nfrag) - plnew%xb(:, 1:nfrag) = frag%xb(:, 1:nfrag) - plnew%vb(:, 1:nfrag) = frag%vb(:, 1:nfrag) - call pl%vb2vh(cb) - call pl%xh2xb(cb) - do i = 1, nfrag - plnew%rh(:,i) = frag%xb(:, i) - cb%xb(:) - plnew%vh(:,i) = frag%vb(:, i) - cb%vb(:) - end do - plnew%mass(1:nfrag) = frag%mass(1:nfrag) - plnew%Gmass(1:nfrag) = param%GU * frag%mass(1:nfrag) - plnew%radius(1:nfrag) = frag%radius(1:nfrag) - plnew%density(1:nfrag) = frag%mass(1:nfrag) / frag%radius(1:nfrag) - call plnew%set_rhill(cb) - - select case(status) - case(DISRUPTION) - plnew%status(1:nfrag) = NEW_PARTICLE - do i = 1, nfrag - write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Disruption", origin_time=system%t, name=newname, & - origin_rh=plnew%rh(:,i), & - origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision) - end do - do i = 1, ncolliders - if (colliders%idx(i) == ibiggest) then - iother = ismallest - else - iother = ibiggest - end if - call pl%info(colliders%idx(i))%set_value(status="Disruption", discard_time=system%t, & - discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) - end do - case(SUPERCATASTROPHIC) - plnew%status(1:nfrag) = NEW_PARTICLE - do i = 1, nfrag - write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=system%t, name=newname, & - origin_rh=plnew%rh(:,i), origin_vh=plnew%vh(:,i), & - collision_id=param%maxid_collision) - end do - do i = 1, ncolliders - if (colliders%idx(i) == ibiggest) then - iother = ismallest - else - iother = ibiggest - end if - call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=system%t, & - discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), & - discard_body_id=iother) - end do - case(HIT_AND_RUN_DISRUPT) - call plnew%info(1)%copy(pl%info(ibiggest)) - plnew%status(1) = OLD_PARTICLE - do i = 2, nfrag - write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Hit and run fragment", origin_time=system%t, name=newname, & - origin_rh=plnew%rh(:,i), origin_vh=plnew%vh(:,i), & - collision_id=param%maxid_collision) - end do - do i = 1, ncolliders - if (colliders%idx(i) == ibiggest) cycle - iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status="Hit and run fragmention", discard_time=system%t, & - discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), & - discard_body_id=iother) - end do - case(MERGED) - call plnew%info(1)%copy(pl%info(ibiggest)) - plnew%status(1) = OLD_PARTICLE - do i = 1, ncolliders - if (colliders%idx(i) == ibiggest) cycle - - iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=system%t, discard_rh=pl%rh(:,i), & - discard_vh=pl%vh(:,i), discard_body_id=iother) - end do - end select - - if (param%lrotation) then - plnew%Ip(:, 1:nfrag) = frag%Ip(:, 1:nfrag) - plnew%rot(:, 1:nfrag) = frag%rot(:, 1:nfrag) - end if - - ! if (param%ltides) then - ! plnew%Q = pl%Q(ibiggest) - ! plnew%k2 = pl%k2(ibiggest) - ! plnew%tlag = pl%tlag(ibiggest) - ! end if - - !Copy over or set integration parameters for new bodies - plnew%lcollision(1:nfrag) = .false. - plnew%ldiscard(1:nfrag) = .false. - plnew%levelg(1:nfrag) = pl%levelg(ibiggest) - plnew%levelm(1:nfrag) = pl%levelm(ibiggest) - - ! Log the properties of the new bodies - call fraggle_io_log_pl(plnew, param) - - ! Append the new merged body to the list - nstart = pl_adds%nbody + 1 - nend = pl_adds%nbody + nfrag - call pl_adds%append(plnew, lsource_mask=[(.true., i=1, nfrag)]) - ! Record how many bodies were added in this event - pl_adds%ncomp(nstart:nend) = plnew%nbody - - ! Add the discarded bodies to the discard list - pl%status(colliders%idx(:)) = MERGED - pl%ldiscard(colliders%idx(:)) = .true. - pl%lcollision(colliders%idx(:)) = .true. - lmask(:) = .false. - lmask(colliders%idx(:)) = .true. - - call plnew%setup(0, param) - deallocate(plnew) - - allocate(plsub, mold=pl) - call pl%spill(plsub, lmask, ldestructive=.false.) - - nstart = pl_discards%nbody + 1 - nend = pl_discards%nbody + ncolliders - call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, ncolliders)]) - - ! Record how many bodies were subtracted in this event - pl_discards%ncomp(nstart:nend) = ncolliders - - call plsub%setup(0, param) - deallocate(plsub) - end associate - end select - end select - - return - end subroutine symba_collision_mergeaddsub - - - module subroutine symba_collision_resolve_fragmentations(self, system, param) - !! author: David A. Minton - !! - !! Process list of collisions, determine the collisional regime, and then create fragments. - !! - implicit none - ! Arguments - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - ! Internals - ! Internals - integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision - logical :: lgoodcollision - integer(I4B) :: i - type(fraggle_colliders) :: colliders !! Fraggle colliders object - type(fraggle_fragments) :: frag !! Fraggle fragmentation system object - - associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2) - select type(pl => system%pl) - class is (symba_pl) - select type (cb => system%cb) - class is (symba_cb) - do i = 1, ncollisions - idx_parent(1) = pl%kin(idx1(i))%parent - idx_parent(2) = pl%kin(idx2(i))%parent - lgoodcollision = symba_collision_consolidate_colliders(pl, cb, param, idx_parent, colliders) - if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLISION)) cycle - - call colliders%regime(frag, system, param) - - select case (frag%regime) - case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - plplcollision_list%status(i) = symba_collision_casedisruption(system, param, colliders, frag) - case (COLLRESOLVE_REGIME_HIT_AND_RUN) - plplcollision_list%status(i) = symba_collision_casehitandrun(system, param, colliders, frag) - case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - plplcollision_list%status(i) = symba_collision_casemerge(system, param, colliders, frag) - case default - write(*,*) "Error in symba_collision, unrecognized collision regime" - call util_exit(FAILURE) - end select - end do - end select - end select - end associate - - return - end subroutine symba_collision_resolve_fragmentations - - - module subroutine symba_collision_resolve_mergers(self, system, param) - !! author: David A. Minton - !! - !! Process list of collisions and merge colliding bodies together. - !! - implicit none - ! Arguments - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - ! Internals - integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision - logical :: lgoodcollision - integer(I4B) :: i - type(fraggle_colliders) :: colliders !! Fraggle colliders object - type(fraggle_fragments) :: frag !! Fraggle fragmentation system object - - associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2) - select type(pl => system%pl) - class is (symba_pl) - select type(cb => system%cb) - class is (symba_cb) - do i = 1, ncollisions - idx_parent(1) = pl%kin(idx1(i))%parent - idx_parent(2) = pl%kin(idx2(i))%parent - lgoodcollision = symba_collision_consolidate_colliders(pl, cb, param, idx_parent, colliders) - if (.not. lgoodcollision) cycle - if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved - - frag%regime = COLLRESOLVE_REGIME_MERGE - frag%mtot = sum(colliders%mass(:)) - frag%mass_dist(1) = frag%mtot - frag%mass_dist(2) = 0.0_DP - frag%mass_dist(3) = 0.0_DP - frag%xbcom(:) = (colliders%mass(1) * colliders%xb(:,1) + colliders%mass(2) * colliders%xb(:,2)) / frag%mtot - frag%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / frag%mtot - plplcollision_list%status(i) = symba_collision_casemerge(system, param, colliders, frag) - end do - end select - end select - end associate - - return - end subroutine symba_collision_resolve_mergers - - - module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, irec) - !! author: David A. Minton - !! - !! Process the pl-pl collision list, then modifiy the massive bodies based on the outcome of the collision - !! - implicit none - ! Arguments - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - real(DP), intent(in) :: t !! Current simulation time - real(DP), intent(in) :: dt !! Current simulation step size - integer(I4B), intent(in) :: irec !! Current recursion level - ! Internals - real(DP) :: Eorbit_before, Eorbit_after - logical :: lplpl_collision - character(len=STRMAX) :: timestr - class(symba_parameters), allocatable :: tmp_param - - associate(plplenc_list => self, plplcollision_list => system%plplcollision_list) - select type(pl => system%pl) - class is (symba_pl) - select type(param) - class is (symba_parameters) - if (plplcollision_list%nenc == 0) return ! No collisions to resolve - ! Make sure that the heliocentric and barycentric coordinates are consistent with each other - call pl%vb2vh(system%cb) - call pl%xh2xb(system%cb) - - ! Get the energy before the collision is resolved - if (param%lenergy) then - call system%get_energy_and_momentum(param) - Eorbit_before = system%te - end if - - do - write(timestr,*) t - call io_log_one_message(FRAGGLE_LOG_OUT, "") - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & - "***********************************************************") - call io_log_one_message(FRAGGLE_LOG_OUT, "Collision between massive bodies detected at time t = " // & - trim(adjustl(timestr))) - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & - "***********************************************************") - allocate(tmp_param, source=param) - if (param%lfragmentation) then - call plplcollision_list%resolve_fragmentations(system, param) - else - call plplcollision_list%resolve_mergers(system, param) - end if - - ! Destroy the collision list now that the collisions are resolved - call plplcollision_list%setup(0_I8B) - - if ((system%pl_adds%nbody == 0) .and. (system%pl_discards%nbody == 0)) exit - - ! Save the add/discard information to file - call system%write_discard(tmp_param) - - ! Rearrange the arrays: Remove discarded bodies, add any new bodies, resort, and recompute all indices and encounter lists - call pl%rearray(system, tmp_param) - - ! Destroy the add/discard list so that we don't append the same body multiple times if another collision is detected - call system%pl_discards%setup(0, param) - call system%pl_adds%setup(0, param) - deallocate(tmp_param) - - ! Check whether or not any of the particles that were just added are themselves in a collision state. This will generate a new plplcollision_list - lplpl_collision = plplenc_list%collision_check(system, param, t, dt, irec) - - if (.not.lplpl_collision) exit - end do - - if (param%lenergy) then - call system%get_energy_and_momentum(param) - Eorbit_after = system%te - system%Ecollisions = system%Ecollisions + (Eorbit_after - Eorbit_before) - end if - - end select - end select - end associate - - return - end subroutine symba_collision_resolve_plplenc - - - module subroutine symba_collision_resolve_pltpenc(self, system, param, t, dt, irec) - !! author: David A. Minton - !! - !! Process the pl-tp collision list, then modifiy the massive bodies based on the outcome of the collision - !! - implicit none - ! Arguments - class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - real(DP), intent(in) :: t !! Current simulation tim - real(DP), intent(in) :: dt !! Current simulation step size - integer(I4B), intent(in) :: irec !! Current recursion level - - ! Make sure coordinate systems are all synced up due to being inside the recursion at this point - call system%pl%vb2vh(system%cb) - call system%tp%vb2vh(system%cb%vb) - call system%pl%b2h(system%cb) - call system%tp%b2h(system%cb) - - ! Discard the collider - call system%tp%discard(system, param) - - return - end subroutine symba_collision_resolve_pltpenc - -end submodule s_symba_collision \ No newline at end of file diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 06c474078..eb891eb23 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -26,7 +26,7 @@ module function symba_collision_casedisruption(system, param, colliders, frag) ! Result integer(I4B) :: status !! Status flag assigned to this outcome ! Internals - integer(I4B) :: i, nfrag + integer(I4B) :: i, ibiggest, nfrag logical :: lfailure character(len=STRMAX) :: message @@ -63,11 +63,16 @@ module function symba_collision_casedisruption(system, param, colliders, frag) select case(frag%regime) case(COLLRESOLVE_REGIME_DISRUPTION) status = DISRUPTION + ibiggest = colliders%idx(maxloc(system%pl%Gmass(colliders%idx(:)), dim=1)) + frag%id(1) = system%pl%id(ibiggest) + frag%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] + param%maxid = frag%id(nfrag) case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) status = SUPERCATASTROPHIC + frag%id(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)] + param%maxid = frag%id(nfrag) end select - frag%id(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)] - param%maxid = frag%id(nfrag) + call symba_collision_mergeaddsub(system, param, colliders, frag, status) end if @@ -720,7 +725,7 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) logical, dimension(system%pl%nbody) :: lmask class(symba_pl), allocatable :: plnew, plsub character(*), parameter :: FRAGFMT = '("Newbody",I0.7)' - character(len=NAMELEN) :: newname + character(len=NAMELEN) :: newname, origin_type select type(pl => system%pl) class is (symba_pl) @@ -757,23 +762,6 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) call plnew%set_rhill(cb) select case(status) - case(DISRUPTION) - plnew%status(1:nfrag) = NEW_PARTICLE - do i = 1, nfrag - write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Disruption", origin_time=system%t, name=newname, & - origin_rh=plnew%rh(:,i), & - origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision) - end do - do i = 1, ncolliders - if (colliders%idx(i) == ibiggest) then - iother = ismallest - else - iother = ibiggest - end if - call pl%info(colliders%idx(i))%set_value(status="Disruption", discard_time=system%t, & - discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) - end do case(SUPERCATASTROPHIC) plnew%status(1:nfrag) = NEW_PARTICLE do i = 1, nfrag @@ -792,19 +780,24 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), & discard_body_id=iother) end do - case(HIT_AND_RUN_DISRUPT) + case(DISRUPTION,HIT_AND_RUN_DISRUPT) + if (status == DISRUPTION) then + write(origin_type,*) "Disruption" + else if (status == HIT_AND_RUN_DISRUPT) then + write(origin_type,*) "Hit and run fragmention" + end if call plnew%info(1)%copy(pl%info(ibiggest)) plnew%status(1) = OLD_PARTICLE do i = 2, nfrag write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Hit and run fragment", origin_time=system%t, name=newname, & + call plnew%info(i)%set_value(origin_type=origin_type, origin_time=system%t, name=newname, & origin_rh=plnew%rh(:,i), origin_vh=plnew%vh(:,i), & collision_id=param%maxid_collision) end do do i = 1, ncolliders if (colliders%idx(i) == ibiggest) cycle iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status="Hit and run fragmention", discard_time=system%t, & + call pl%info(colliders%idx(i))%set_value(status=origin_type, discard_time=system%t, & discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), & discard_body_id=iother) end do