diff --git a/src/io/io.f90 b/src/io/io.f90 index 2fdc0efa9..2d8334b60 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -172,6 +172,7 @@ module subroutine io_dump_system(self, param) real(DP) :: deltawall, wallperstep, tfrac integer(I8B) :: clock_count, count_rate, count_max character(*), parameter :: statusfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, "; Number of active pl, tp = ", I5, ", ", I5)' + character(*), parameter :: symbastatfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, "; Number of active plm, pl, tp = ", I5, ", ", I5, ", ", I5)' character(len=*), parameter :: walltimefmt = '(" Wall time (s): ", es12.5, "; Wall time/step in this interval (s): ", es12.5)' logical, save :: lfirst = .true. real(DP), save :: start, finish @@ -209,7 +210,12 @@ module subroutine io_dump_system(self, param) deltawall = clock_count / (count_rate * 1.0_DP) - finish wallperstep = deltawall / param%istep_dump finish = clock_count / (count_rate * 1.0_DP) - write(*, statusfmt) param%t, tfrac, self%pl%nbody, self%tp%nbody + select type(pl => self%pl) + class is (symba_pl) + write(*, symbastatfmt) param%t, tfrac, pl%nplm, pl%nbody, self%tp%nbody + class default + write(*, statusfmt) param%t, tfrac, pl%nbody, self%tp%nbody + end select write(*, walltimefmt) finish - start, wallperstep if (param%lenergy) call self%conservation_report(param, lterminal=.true.) end if diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index af8e09055..b259dfba9 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -220,7 +220,7 @@ module function symba_collision_casemerge(system, param, family, x, v, mass, rad ! Result integer(I4B) :: status !! Status flag assigned to this outcome ! Internals - integer(I4B) :: i, j, ibiggest, nfamily + integer(I4B) :: i, j, k, ibiggest, nfamily real(DP) :: volume_new, pe real(DP), dimension(NDIM) :: xc, vc, xcrossv real(DP), dimension(2) :: vol @@ -283,8 +283,27 @@ module function symba_collision_casemerge(system, param, family, x, v, mass, rad param%Ecollisions = param%Ecollisions + pe param%Euntracked = param%Euntracked - pe + ! Update any encounter lists that have the removedbodies in them so that they instead point to the new + do k = 1, system%plplenc_list%nenc + do j = 1, nfamily + i = family(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, family, id_frag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, status) + end select diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index c5796ccc8..c50e0d373 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -468,6 +468,7 @@ module subroutine symba_util_rearray_pl(self, system, param) logical, dimension(:), allocatable :: lmask class(symba_plplenc), allocatable :: plplenc_old logical :: lencounter + integer(I4B), dimension(:), allocatable :: levelg_orig_pl, levelm_orig_pl, levelg_orig_tp, levelm_orig_tp, nplenc_orig_pl, nplenc_orig_tp, ntpenc_orig_pl associate(pl => self, pl_adds => system%pl_adds, nadd => system%pl_adds%nbody) @@ -522,12 +523,27 @@ module subroutine symba_util_rearray_pl(self, system, param) pl%kin(1:npl)%nchild = 0 pl%kin(1:npl)%parent = [(i, i=1, npl)] - ! Re-build the encounter list + ! Re-build the zero-level encounter list, being sure to save the original level information for all bodies + allocate(levelg_orig_pl, source=pl%levelg) + allocate(levelm_orig_pl, source=pl%levelm) + allocate(nplenc_orig_pl, source=pl%nplenc) + allocate(ntpenc_orig_pl, source=pl%ntpenc) lencounter = pl%encounter_check(system, param%dt, 0) - select type(tp => system%tp) - class is (symba_tp) - lencounter = tp%encounter_check(system, param%dt, 0) - end select + if (system%tp%nbody > 0) then + select type(tp => system%tp) + class is (symba_tp) + allocate(levelg_orig_tp, source=tp%levelg) + allocate(levelm_orig_tp, source=tp%levelm) + allocate(nplenc_orig_tp, source=tp%nplenc) + lencounter = tp%encounter_check(system, param%dt, 0) + call move_alloc(levelg_orig_tp, tp%levelg) + call move_alloc(levelm_orig_tp, tp%levelm) + call move_alloc(nplenc_orig_tp, tp%nplenc) + end select + end if + call move_alloc(levelg_orig_pl, pl%levelg) + call move_alloc(levelm_orig_pl, pl%levelm) + call move_alloc(nplenc_orig_pl, pl%nplenc) associate(idnew1 => system%plplenc_list%id1, idnew2 => system%plplenc_list%id2, idold1 => plplenc_old%id1, idold2 => plplenc_old%id2) do k = 1, system%plplenc_list%nenc