diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 0fb803d8d..ab666fb00 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -43,6 +43,7 @@ movie_title_list = ["Head-on Disrutption", "Off-axis Supercatastrophic", "Hit and Run"] movie_titles = dict(zip(available_movie_styles, movie_title_list)) +# These initial conditions were generated by trial and error pos_vectors = {"disruption_headon" : [np.array([1.0, -1.807993e-05, 0.0]), np.array([1.0, 1.807993e-05 ,0.0])], "supercatastrophic_off_axis": [np.array([1.0, -4.2e-05, 0.0]), @@ -67,7 +68,7 @@ np.array([0.0, 0.0, 1.0e5])] } -body_Gmass = {"disruption_headon" : [1e-7, 7e-10], +body_Gmass = {"disruption_headon" : [1e-7, 1e-10], "supercatastrophic_off_axis": [1e-7, 1e-8], "hitandrun" : [1e-7, 7e-10] } @@ -131,8 +132,10 @@ def animate(i,ds,movie_title): sim.add_body(Gmass=body_Gmass[style], radius=body_radius[style], xh=pos_vectors[style], vh=vel_vectors[style], rot=rot_vectors[style]) # Set fragmentation parameters - sim.set_parameter(fragmentation = True, gmtiny=1e-11, minimum_fragment_gmass=1e-11) - sim.run(dt=1e-8, tstop=1e-5) + minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body + gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades + sim.set_parameter(fragmentation = True, gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass) + sim.run(dt=1e-8, tstop=1.e-5) # Calculate the number of frames in the dataset. nframes = int(sim.data['time'].size) @@ -144,5 +147,5 @@ def animate(i,ds,movie_title): # Set up the figure and the animation. fig, ax = plt.subplots(figsize=(4,4)) # Generate the movie. - ani = animation.FuncAnimation(fig, animate, interval=1, frames=nframes, repeat=False) + ani = animation.FuncAnimation(fig, animate, fargs=(sim.data, movie_titles[style]), interval=1, frames=nframes, repeat=False) ani.save(movie_filename, fps=60, dpi=300, extra_args=['-vcodec', 'libx264']) \ No newline at end of file diff --git a/examples/Fragmentation/cb.in b/examples/Fragmentation/cb.in deleted file mode 100644 index a1275bf77..000000000 --- a/examples/Fragmentation/cb.in +++ /dev/null @@ -1,7 +0,0 @@ -Sun -39.47841760435743 ! G*Mass -0.005 ! Radius -0.0 ! J2 -0.0 ! J4 -0.4 0.4 0.4 !Ip -0.0 0.0 0.0 !rot !11.2093063 -38.75937204 82.25088158 ! rot (radian / year) \ No newline at end of file diff --git a/examples/Fragmentation/disruption_headon.in b/examples/Fragmentation/disruption_headon.in deleted file mode 100644 index 72c14edd1..000000000 --- a/examples/Fragmentation/disruption_headon.in +++ /dev/null @@ -1,13 +0,0 @@ -2 -Body1 1e-07 0.0009 -7e-06 -1.0 -1.807993e-05 0.0 --2.562596e-04 6.280005 0.0 -0.4 0.4 0.4 !Ip -0.0 0.0 0.0 !rot -Body2 7e-10 0.0004 -3.25e-06 -1.0 1.807993e-05 0.0 --2.562596e-04 -6.280005 0.0 -0.4 0.4 0.4 !Ip -0.0 0.0 0.0 !rot diff --git a/examples/Fragmentation/hitandrun.in b/examples/Fragmentation/hitandrun.in deleted file mode 100644 index 554d953ee..000000000 --- a/examples/Fragmentation/hitandrun.in +++ /dev/null @@ -1,13 +0,0 @@ -2 -Body1 1e-07 0.0009 -7e-06 -1.0 -4.20E-05 0.0 -0.00 6.28 0.0 -0.4 0.4 0.4 !Ip -0.0 0.0 6.0e4 !rot -Body2 7e-10 0.0004 -3.25e-06 -1.0 4.20E-05 0.0 --1.50 -6.28 0.0 -0.4 0.4 0.4 !Ip -0.0 0.0 1.0e5 !rot diff --git a/examples/Fragmentation/supercatastrophic_off_axis.in b/examples/Fragmentation/supercatastrophic_off_axis.in deleted file mode 100644 index cd09c9070..000000000 --- a/examples/Fragmentation/supercatastrophic_off_axis.in +++ /dev/null @@ -1,13 +0,0 @@ -2 -Body1 1e-07 0.0009 -7e-06 -1.0 -4.2e-05 0.0 -0.00 6.28 0.0 -0.4 0.4 0.4 !Ip -0.0 0.0 -6.0e4 !rot -Body2 1e-08 0.0004 -3.25e-06 -1.0 4.2e-05 0.0 -1.00 -6.28 0.0 -0.4 0.4 0.4 !Ip -0.0 0.0 1.0e5 !rot diff --git a/examples/Fragmentation/tp.in b/examples/Fragmentation/tp.in deleted file mode 100644 index 3c6f40630..000000000 --- a/examples/Fragmentation/tp.in +++ /dev/null @@ -1,10 +0,0 @@ -!! Copyright 2022 - 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. - -0 diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index fbd0a3360..4a2687af7 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -31,6 +31,7 @@ Literal, Dict, List, + Tuple, Any ) @@ -387,7 +388,7 @@ def _type_scrub(output_data): if "nplm" in self.data: post_message += f" nplm: {self.data['nplm'].values[0]}" if self.param['ENERGY']: - post_message += f" dL/L0: {0.0:.5e} dE/|E0|: {0.0:.5e}" + post_message += f" dL/L0: {0.0:.5e} dE/|E0|: {0.0:+.5e}" pbar = tqdm(total=noutput, desc=pre_message, postfix=post_message, bar_format='{l_bar}{bar}{postfix}') try: with subprocess.Popen(shlex.split(cmd), @@ -410,7 +411,7 @@ def _type_scrub(output_data): if "LTOTERR" in output_data: post_message += f" dL/L0: {output_data['LTOTERR']:.5e}" if "ETOTERR" in output_data: - post_message += f" dE/|E0|: {output_data['ETOTERR']:.5e}" + post_message += f" dE/|E0|: {output_data['ETOTERR']:+.5e}" interval = output_data['ILOOP'] - iloop if interval > 0: pbar.update(interval) @@ -2278,8 +2279,8 @@ def add_body(self, v4: float | List[float] | npt.NDArray[np.float_] | None = None, v5: float | List[float] | npt.NDArray[np.float_] | None = None, v6: float | List[float] | npt.NDArray[np.float_] | None = None, - xh: List[float] | npt.NDArray[np.float_] | None = None, - vh: List[float] | npt.NDArray[np.float_] | None = None, + xh: List[float] | List[npt.NDArray[np.float_]] | npt.NDArray[np.float_] | None = None, + vh: List[float] | List[npt.NDArray[np.float_]] | npt.NDArray[np.float_] | None = None, mass: float | List[float] | npt.NDArray[np.float_] | None=None, Gmass: float | List[float] | npt.NDArray[np.float_] | None=None, radius: float | List[float] | npt.NDArray[np.float_] | None=None, @@ -2291,7 +2292,7 @@ def add_body(self, rotx: float | List[float] | npt.NDArray[np.float_] | None=None, roty: float | List[float] | npt.NDArray[np.float_] | None=None, rotz: float | List[float] | npt.NDArray[np.float_] | None=None, - rot: List[float] | npt.NDArray[np.float_] | None=None, + rot: List[float] | List[npt.NDArray[np.float_]] | npt.NDArray[np.float_] | None=None, J2: float | List[float] | npt.NDArray[np.float_] | None=None, J4: float | List[float] | npt.NDArray[np.float_] | None=None): """ diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index d59e2a9b7..3ec23ef99 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -57,8 +57,12 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa end if f_spin = F_SPIN_FIRST - lk_plpl = allocated(pl%k_plpl) - if (lk_plpl) deallocate(pl%k_plpl) + if (param%lflatten_interactions) then + lk_plpl = allocated(pl%k_plpl) + if (lk_plpl) deallocate(pl%k_plpl) + else + lk_plpl = .false. + end if call frag%set_natural_scale(colliders) @@ -250,19 +254,21 @@ subroutine fraggle_generate_spins(frag, f_spin, lfailure) frag%rot(:,:) = 0.0_DP frag%ke_spin = 0.0_DP - do i = 1, nfrag - ! Convert a fraction (f_spin) of either the remaining angular momentum or kinetic energy budget into spin, whichever gives the smaller rotation so as not to blow any budgets - rot_ke(:) = sqrt(2 * f_spin * frag%ke_budget / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i))) & - * L_remainder(:) / norm2(L_remainder(:)) - rot_L(:) = f_spin * L_remainder(:) / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i)) - if (norm2(rot_ke) < norm2(rot_L)) then - frag%rot(:,i) = rot_ke(:) - else - frag%rot(:, i) = rot_L(:) - end if - frag%ke_spin = frag%ke_spin + frag%mass(i) * frag%Ip(3, i) * frag%radius(i)**2 & - * dot_product(frag%rot(:, i), frag%rot(:, i)) - end do + if (norm2(L_remainder(:)) > FRAGGLE_LTOL) then + do i = 1, nfrag + ! Convert a fraction (f_spin) of either the remaining angular momentum or kinetic energy budget into spin, whichever gives the smaller rotation so as not to blow any budgets + rot_ke(:) = sqrt(2 * f_spin * frag%ke_budget / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i))) & + * L_remainder(:) / norm2(L_remainder(:)) + rot_L(:) = f_spin * L_remainder(:) / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i)) + if (norm2(rot_ke) < norm2(rot_L)) then + frag%rot(:,i) = rot_ke(:) + else + frag%rot(:, i) = rot_L(:) + end if + frag%ke_spin = frag%ke_spin + frag%mass(i) * frag%Ip(3, i) * frag%radius(i)**2 & + * dot_product(frag%rot(:, i), frag%rot(:, i)) + end do + end if frag%ke_spin = 0.5_DP * frag%ke_spin lfailure = ((frag%ke_budget - frag%ke_spin - frag%ke_orbit) < 0.0_DP) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 8d9594974..e03e30eb5 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -87,7 +87,6 @@ module subroutine fraggle_util_construct_temporary_system(frag, system, param, t !! Author: David A. Minton !! !! Constructs a temporary internal system consisting of active bodies and additional fragments. This internal temporary system is used to calculate system energy with and without fragments - !! and optionally including fragments. implicit none ! Arguments class(fraggle_fragments), intent(in) :: frag !! Fraggle fragment system object @@ -147,7 +146,6 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the system, with colliders included and fragments excluded or vice versa ! Internals - logical, dimension(:), allocatable :: lexclude class(swiftest_nbody_system), allocatable, save :: tmpsys class(swiftest_parameters), allocatable, save :: tmpparam integer(I4B) :: npl_before, npl_after @@ -162,24 +160,24 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para npl_before = pl%nbody npl_after = npl_before + nfrag - ! Build the exluded body logical mask - allocate(lexclude(npl_after)) if (lbefore) then - lexclude(1:npl_before) = .false. - lexclude(npl_before+1:npl_after) = .true. call fraggle_util_construct_temporary_system(frag, system, param, tmpsys, tmpparam) + ! Build the exluded body logical mask for the *before* case: Only the original bodies are used to compute energy and momentum + tmpsys%pl%status(colliders%idx(1:colliders%ncoll)) = ACTIVE + tmpsys%pl%status(npl_before+1:npl_after) = INACTIVE else - lexclude(1:npl_after) = .false. - lexclude(colliders%idx(1:colliders%ncoll)) = .true. if (.not.allocated(tmpsys)) then write(*,*) "Error in fraggle_util_get_energy_momentum. " // & " This must be called with lbefore=.true. at least once before calling it with lbefore=.false." call util_exit(FAILURE) end if + ! Build the exluded body logical mask for the *after* case: Only the new bodies are used to compute energy and momentum call fraggle_util_add_fragments_to_system(frag, colliders, tmpsys, tmpparam) + tmpsys%pl%status(colliders%idx(1:colliders%ncoll)) = INACTIVE + tmpsys%pl%status(npl_before+1:npl_after) = ACTIVE end if - call tmpsys%pl%flatten(param) + if (param%lflatten_interactions) call tmpsys%pl%flatten(param) call tmpsys%get_energy_and_momentum(param) diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index 621ea80a6..ed7119d8b 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -23,7 +23,6 @@ module subroutine util_get_energy_momentum_system(self, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i - integer(I8B) :: nplpl real(DP) :: kecb, kespincb real(DP), dimension(self%pl%nbody) :: kepl, kespinpl real(DP), dimension(self%pl%nbody) :: Lplorbitx, Lplorbity, Lplorbitz @@ -32,7 +31,6 @@ module subroutine util_get_energy_momentum_system(self, param) real(DP) :: hx, hy, hz associate(system => self, pl => self%pl, npl => self%pl%nbody, cb => self%cb) - nplpl = pl%nplpl system%Lorbit(:) = 0.0_DP system%Lspin(:) = 0.0_DP system%Ltot(:) = 0.0_DP @@ -89,7 +87,7 @@ module subroutine util_get_energy_momentum_system(self, param) end if if (param%lflatten_interactions) then - call util_get_energy_potential_flat(npl, nplpl, pl%k_plpl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%xb, system%pe) + call util_get_energy_potential_flat(npl, pl%nplpl, pl%k_plpl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%xb, system%pe) else call util_get_energy_potential_triangular(npl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%xb, system%pe) end if