diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 807f22419..6932223a8 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -127,7 +127,7 @@ def encounter_combiner(sim): ds = xr.combine_nested([data,enc],concat_dim='time').sortby("time").interpolate_na(dim="time") # Interpolate in time to make a smooth, constant time step dataset - smooth_time = np.linspace(start=tgood.isel(time=0), stop=tgood.isel(time=-1), num=2400) + smooth_time = np.linspace(start=tgood.isel(time=0), stop=ds.time[-1], num=2400) ds = ds.interp(time=smooth_time) return ds diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index e66848183..e4fa80049 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -2775,6 +2775,9 @@ def _preprocess(ds, param): # Remove any overlapping time values tgood,tid = np.unique(self.encounters.time,return_index=True) self.encounters = self.encounters.isel(time=tid) + # Remove any NaN values + tgood=self.encounters.time.where(~np.isnan(self.encounters.time),drop=True) + self.encounters = self.encounters.sel(time=tgood) return diff --git a/src/collision/collision_check.f90 b/src/collision/collision_check.f90 index eb5c26f8d..1f88c98aa 100644 --- a/src/collision/collision_check.f90 +++ b/src/collision/collision_check.f90 @@ -79,9 +79,7 @@ module subroutine collision_check_plpl(self, nbody_system, param, t, dt, irec, l real(DP), dimension(NDIM) :: xr, vr integer(I4B) :: i, j, k, nenc real(DP) :: rlim, Gmtot - logical :: isplpl, lany_closest - character(len=STRMAX) :: timestr, idstri, idstrj, message - class(collision_list_plpl), allocatable :: tmp + logical :: lany_closest lany_collision = .false. if (self%nenc == 0) return diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index ce49e0cbd..cba31cfb8 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -17,8 +17,8 @@ module subroutine encounter_io_netcdf_dump(self, param) !! Dumps the time history of an encounter to file. implicit none ! Arguments - class(encounter_storage(*)), intent(inout) :: self !! Encounter storage object - class(base_parameters), intent(inout) :: param !! Current run configuration parameters + class(encounter_storage(*)), intent(inout) :: self !! Encounter storage object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i @@ -33,7 +33,7 @@ module subroutine encounter_io_netcdf_dump(self, param) write(nc%file_name, '("encounter_",I0.6,".nc")') nc%file_number call nc%initialize(param) - do i = 1, self%nframes + do i = 1, self%iframe if (allocated(self%frame(i)%item)) then select type(snapshot => self%frame(i)%item) class is (encounter_snapshot) diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index 45e6f48fb..0377c290c 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -332,8 +332,6 @@ module subroutine encounter_util_setup_list(self, n) ! Arguments class(encounter_list), intent(inout) :: self !! Swiftest encounter structure integer(I8B), intent(in) :: n !! Number of encounters to allocate space for - ! Internals - integer(I8B) :: i if (n < 0) return call self%dealloc() @@ -450,7 +448,6 @@ subroutine encounter_util_save_snapshot(encounter_history, snapshot) end do deallocate(encounter_history) call move_alloc(tmp,encounter_history) - nnew = nbig end if ! Find out which time slot this belongs in by searching for an existing slot diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index aa14c8bd2..f00080f6b 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -131,7 +131,6 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur lfailure_local = .false. - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle try " // trim(adjustl(message))) ! Use the disruption collision model to generate initial conditions ! Compute the "before" energy/momentum and the budgets @@ -140,6 +139,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur do try = 1, MAXTRY self%fail_scale = (fail_scale_initial)**try write(message,*) try + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle try " // trim(adjustl(message))) call fraggle_generate_pos_vec(self) call fraggle_generate_rot_vec(self) call fraggle_generate_vel_vec(self,lfailure_local) @@ -187,7 +187,6 @@ module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t) logical :: lpure character(len=STRMAX) :: message - select type(nbody_system) class is (swiftest_nbody_system) select type(pl => nbody_system%pl) @@ -250,12 +249,12 @@ module subroutine fraggle_generate_pos_vec(collider) ! Internals real(DP) :: dis real(DP), dimension(NDIM,2) :: fragment_cloud_center - real(DP), dimension(2) :: fragment_cloud_radius + real(DP), dimension(2) :: fragment_cloud_radius, rdistance logical, dimension(collider%fragments%nbody) :: loverlap integer(I4B) :: i, j, loop logical :: lcat, lhitandrun integer(I4B), parameter :: MAXLOOP = 10000 - real(DP) :: rdistance + real(DP), parameter :: rdistance_scale_factor = 0.01_DP ! Scale factor to apply to distance scaling in the event of a fail associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) lcat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) @@ -263,16 +262,13 @@ module subroutine fraggle_generate_pos_vec(collider) ! We will treat the first two fragments of the list as special cases. ! Place the first two bodies at the centers of the two fragment clouds, but be sure they are sufficiently far apart to avoid overlap - rdistance = .mag. (impactors%rc(:,2) - impactors%rc(:,1)) - sum(fragments%radius(1:2)) - rdistance = 0.5_DP*rdistance - fragment_cloud_radius(:) = impactors%radius(:) - loverlap(:) = .true. + rdistance(:) = rdistance_scale_factor * .mag.impactors%vc(:,:) do loop = 1, MAXLOOP if (.not.any(loverlap(:))) exit - fragment_cloud_center(:,1) = impactors%rc(:,1) !+ rdistance * impactors%bounce_unit(:) - fragment_cloud_center(:,2) = impactors%rc(:,2) !- rdistance * impactors%bounce_unit(:) + fragment_cloud_center(:,1) = impactors%rc(:,1) - rdistance(1) * impactors%bounce_unit(:) + fragment_cloud_center(:,2) = impactors%rc(:,2) + rdistance(2) * impactors%bounce_unit(:) do concurrent(i = 1:nfrag, loverlap(i)) if (i < 3) then fragments%rc(:,i) = fragment_cloud_center(:,i) @@ -301,7 +297,7 @@ module subroutine fraggle_generate_pos_vec(collider) loverlap(i) = loverlap(i) .or. (dis <= (fragments%radius(i) + fragments%radius(j))) end do end do - rdistance = rdistance * collider%fail_scale + rdistance(:) = rdistance(:) * collider%fail_scale fragment_cloud_radius(:) = fragment_cloud_radius(:) * collider%fail_scale end do diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 8d641485e..fcbb734fd 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -34,19 +34,17 @@ module subroutine symba_step_system(self, param, t, dt) class is (symba_tp) select type(param) class is (swiftest_parameters) - associate(encounter_history => self%encounter_history) - call self%reset(param) - lencounter = pl%encounter_check(param, self, dt, 0) .or. tp%encounter_check(param, self, dt, 0) - if (lencounter) then - if (param%lenc_save_trajectory) call encounter_history%take_snapshot(param, self, t, "trajectory") - call self%interp(param, t, dt) - if (param%lenc_save_trajectory) call encounter_history%take_snapshot(param, self, t+dt, "trajectory") - else - self%irec = -1 - call helio_step_system(self, param, t, dt) - end if - param%lfirstkick = pl%lfirst - end associate + call self%reset(param) + lencounter = pl%encounter_check(param, self, dt, 0) .or. tp%encounter_check(param, self, dt, 0) + if (lencounter) then + if (param%lenc_save_trajectory) call self%encounter_history%take_snapshot(param, self, t, "trajectory") + call self%interp(param, t, dt) + if (param%lenc_save_trajectory) call self%encounter_history%take_snapshot(param, self, t+dt, "trajectory") + else + self%irec = -1 + call helio_step_system(self, param, t, dt) + end if + param%lfirstkick = pl%lfirst end select end select end select @@ -192,7 +190,7 @@ recursive module subroutine symba_step_recur_system(self, param, t, ireci) class is (symba_list_plpl) select type(pltp_encounter => self%pltp_encounter) class is (symba_list_pltp) - associate(nbody_system => self, lplpl_collision => plpl_encounter%lcollision, lpltp_collision => pltp_encounter%lcollision, encounter_history => self%encounter_history) + associate(nbody_system => self, lplpl_collision => plpl_encounter%lcollision, lpltp_collision => pltp_encounter%lcollision) nbody_system%irec = ireci dtl = param%dt / (NTENC**ireci) dth = 0.5_DP * dtl @@ -249,7 +247,7 @@ recursive module subroutine symba_step_recur_system(self, param, t, ireci) if (lplpl_collision) call plpl_encounter%resolve_collision(nbody_system, param, t+j*dtl, dtl, ireci) if (lpltp_collision) call pltp_encounter%resolve_collision(nbody_system, param, t+j*dtl, dtl, ireci) end if - if (param%lenc_save_trajectory) call encounter_history%take_snapshot(param, self, t+j*dtl, "trajectory") + if (param%lenc_save_trajectory) call self%encounter_history%take_snapshot(param, self, t+j*dtl, "trajectory") call self%set_recur_levels(ireci)