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