From ad8520d83aad1aadad01e0110fe959badef176ff Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 2 Aug 2021 17:58:37 -0400 Subject: [PATCH] Made somewhat functional SyMBA version of the discard method for pl --- src/io/io.f90 | 2 +- src/modules/symba_classes.f90 | 1 - src/setup/setup.f90 | 1 - src/symba/symba_collision.f90 | 6 +++--- src/symba/symba_drift.f90 | 8 ++++---- src/symba/symba_io.f90 | 27 ++++++++++++++++----------- src/symba/symba_kick.f90 | 6 +++--- src/util/util_append.f90 | 2 +- src/util/util_coord.f90 | 4 ++-- src/util/util_resize.f90 | 2 +- 10 files changed, 31 insertions(+), 28 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index e3556b11c..2900fc7f4 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -1074,7 +1074,7 @@ module subroutine io_write_discard(self, param) write(LUN, HDRFMT) param%t, nsp, param%lbig_discard do i = 1, nsp - write(LUN, NAMEFMT) sub, tp_discards%id(i), tp_discards%status(i) + write(LUN, NAMEFMT) SUB, tp_discards%id(i), tp_discards%status(i) write(LUN, VECFMT) tp_discards%xh(1, i), tp_discards%xh(2, i), tp_discards%xh(3, i) write(LUN, VECFMT) tp_discards%vh(1, i), tp_discards%vh(2, i), tp_discards%vh(3, i) end do diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 2094b533a..c4b399471 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -164,7 +164,6 @@ module symba_classes class(symba_pl), allocatable :: mergesub_list !! List of subtracted bodies in mergers or collisions class(symba_pltpenc), allocatable :: pltpenc_list !! List of massive body-test particle encounters in a single step class(symba_plplenc), allocatable :: plplenc_list !! List of massive body-massive body encounters in a single step - class(symba_pl), allocatable :: pl_discards !! Discarded test particle data structure integer(I4B) :: irec !! System recursion level contains procedure :: write_discard => symba_io_write_discard !! Write out information about discarded and merged planets and test particles in SyMBA diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index f1e44f19c..7de8ce5c3 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -53,7 +53,6 @@ module subroutine setup_construct_system(system, param) allocate(symba_cb :: system%cb) allocate(symba_pl :: system%pl) allocate(symba_tp :: system%tp) - allocate(symba_pl :: system%pl_discards) allocate(symba_tp :: system%tp_discards) allocate(symba_pl :: system%mergeadd_list) allocate(symba_pl :: system%mergesub_list) diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 7c0a05c91..ed2e1e0a1 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -67,11 +67,11 @@ module subroutine symba_collision_check_plplenc(self, system, param, t, dt, irec lmask(:) = .false. lmask(ind1(k)) = .not.pl%lcollision(ind1(k)) lmask(ind2(k)) = .not.pl%lcollision(ind2(k)) - call system%mergesub_list%append(pl, lmask) - ! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step pl%lcollision([ind1(k), ind2(k)]) = .true. - pl%ldiscard([ind1(k), ind2(k)]) = .true. + pl%lcollision([ind1(k), ind2(k)]) = .true. + pl%status([ind1(k), ind2(k)]) = COLLISION + call system%mergesub_list%append(pl, lmask) end do end if diff --git a/src/symba/symba_drift.f90 b/src/symba/symba_drift.f90 index ac06cbb6a..c4efee05f 100644 --- a/src/symba/symba_drift.f90 +++ b/src/symba/symba_drift.f90 @@ -17,9 +17,9 @@ module subroutine symba_drift_pl(self, system, param, dt) select type(system) class is (symba_nbody_system) - self%lmask(:) = self%status(:) == ACTIVE .and. self%levelg(:) == system%irec + self%lmask(:) = self%status(:) /= INACTIVE .and. self%levelg(:) == system%irec call helio_drift_body(self, system, param, dt) - self%lmask(:) = self%status(:) == ACTIVE + self%lmask(:) = self%status(:) /= INACTIVE end select return @@ -41,9 +41,9 @@ module subroutine symba_drift_tp(self, system, param, dt) select type(system) class is (symba_nbody_system) - self%lmask(:) = self%status(:) == ACTIVE .and. self%levelg(:) == system%irec + self%lmask(:) = self%status(:) /= INACTIVE .and. self%levelg(:) == system%irec call helio_drift_body(self, system, param, dt) - self%lmask(:) = self%status(:) == ACTIVE + self%lmask(:) = self%status(:) /= INACTIVE end select return diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 97e4ab3f6..f35e45408 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -225,10 +225,10 @@ module subroutine symba_io_write_discard(self, param) character(*), parameter :: PLNAMEFMT = '(I8, 2(1X, E23.16))' class(swiftest_body), allocatable :: pltemp - associate(pl_discards => self%pl_discards, nsppl => self%pl_discards%nbody, pl => self%pl, npl => self%pl%nbody) + associate(pl => self%pl, npl => self%pl%nbody, mergesub_list => self%mergesub_list, mergeadd_list => self%mergeadd_list) if (self%tp_discards%nbody > 0) call io_write_discard(self, param) - if (nsppl == 0) return + if (mergesub_list%nbody == 0) return select case(param%out_stat) case('APPEND') open(unit = LUN, file = param%discard_out, status = 'OLD', position = 'APPEND', form = 'FORMATTED', iostat = ierr) @@ -239,14 +239,19 @@ module subroutine symba_io_write_discard(self, param) call util_exit(FAILURE) end select lfirst = .false. - if (param%lgr) call pl_discards%pv2v(param) - - ! write(LUN, HDRFMT) param%t, nsp, param%lbig_discard - ! do i = 1, nsp - ! write(LUN, NAMEFMT) sub, tp_discards%id(i), tp_discards%status(i) - ! write(LUN, VECFMT) tp_discards%xh(1, i), tp_discards%xh(2, i), tp_discards%xh(3, i) - ! write(LUN, VECFMT) tp_discards%vh(1, i), tp_discards%vh(2, i), tp_discards%vh(3, i) - ! end do + if (param%lgr) then + call mergesub_list%pv2v(param) + call mergeadd_list%pv2v(param) + end if + + write(LUN, HDRFMT) param%t, mergesub_list%nbody, param%lbig_discard + do i = 1, mergesub_list%nbody + write(LUN, NAMEFMT) SUB, mergesub_list%id(i), mergesub_list%status(i) + write(LUN, VECFMT) mergesub_list%xh(1, i), mergesub_list%xh(2, i), mergesub_list%xh(3, i) + write(LUN, VECFMT) mergesub_list%vh(1, i), mergesub_list%vh(2, i), mergesub_list%vh(3, i) + end do + + ! This is incomplete until the mergeadd_list methods are completed ! if (param%lbig_discard) then ! if (param%lgr) then ! allocate(pltemp, source = pl) @@ -265,7 +270,7 @@ module subroutine symba_io_write_discard(self, param) ! end do ! deallocate(vh) ! end if - ! close(LUN) + close(LUN) end associate return diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index aebb6bb2b..8625b3d81 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -124,8 +124,8 @@ module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn) select type(tp => system%tp) class is (symba_tp) associate(ind1 => self%index1, ind2 => self%index2) - if (pl%nbody > 0) pl%lmask(:) = pl%status(:) == ACTIVE - if (tp%nbody > 0) tp%lmask(:) = tp%status(:) == ACTIVE + if (pl%nbody > 0) pl%lmask(:) = pl%status(:) /= INACTIVE + if (tp%nbody > 0) tp%lmask(:) = tp%status(:) /= INACTIVE irm1 = irec - 1 if (sgn < 0) then @@ -145,7 +145,7 @@ module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn) else lgoodlevel = (pl%levelg(ind1(k)) >= irm1) .and. (tp%levelg(ind2(k)) >= irm1) end if - if ((self%status(k) == ACTIVE) .and. lgoodlevel) then + if ((self%status(k) /= INACTIVE) .and. lgoodlevel) then if (isplpl) then ri = ((pl%rhill(ind1(k)) + pl%rhill(ind2(k)))**2) * (RHSCALE**2) * (RSHELL**(2*irecl)) rim1 = ri * (RSHELL**2) diff --git a/src/util/util_append.f90 b/src/util/util_append.f90 index a13103cfa..0f7ac0bde 100644 --- a/src/util/util_append.f90 +++ b/src/util/util_append.f90 @@ -232,7 +232,7 @@ module subroutine util_append_body(self, source, lsource_mask) call util_append(self%omega, source%omega, lsource_mask) call util_append(self%capm, source%capm, lsource_mask) - self%nbody = count(self%status(:) == ACTIVE) + self%nbody = count(self%status(:) /= INACTIVE) return end subroutine util_append_body diff --git a/src/util/util_coord.f90 b/src/util/util_coord.f90 index bdc772d21..938d04951 100644 --- a/src/util/util_coord.f90 +++ b/src/util/util_coord.f90 @@ -54,7 +54,7 @@ module subroutine util_coord_h2b_tp(self, cb) associate(ntp => self%nbody, xbcb => cb%xb, vbcb => cb%vb, status => self%status, & xb => self%xb, xh => self%xh, vb => self%vb, vh => self%vh) - where(status(1:ntp) == ACTIVE) + where(status(1:ntp) /= INACTIVE) xb(1, 1:ntp) = xh(1, 1:ntp) + xbcb(1) xb(2, 1:ntp) = xh(2, 1:ntp) + xbcb(2) xb(3, 1:ntp) = xh(3, 1:ntp) + xbcb(3) @@ -109,7 +109,7 @@ module subroutine util_coord_b2h_tp(self, cb) associate(ntp => self%nbody, xbcb => cb%xb, vbcb => cb%vb, xb => self%xb, xh => self%xh, & vb => self%vb, vh => self%vh, status => self%status) - where(status(1:ntp) == ACTIVE) + where(status(1:ntp) /= INACTIVE) xh(1, 1:ntp) = xb(1, 1:ntp) - xbcb(1) xh(2, 1:ntp) = xb(2, 1:ntp) - xbcb(2) xh(3, 1:ntp) = xb(3, 1:ntp) - xbcb(3) diff --git a/src/util/util_resize.f90 b/src/util/util_resize.f90 index 4a84a003b..9abd66189 100644 --- a/src/util/util_resize.f90 +++ b/src/util/util_resize.f90 @@ -201,7 +201,7 @@ module subroutine util_resize_body(self, nnew) call util_resize(self%capom, nnew) call util_resize(self%omega, nnew) call util_resize(self%capm, nnew) - self%nbody = count(self%status(1:nnew) == ACTIVE) + self%nbody = count(self%status(1:nnew) /= INACTIVE) return end subroutine util_resize_body