From 155dc899f6ce98486e93fd53dbe299cc6a36e7d6 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 6 Dec 2022 10:49:42 -0500 Subject: [PATCH] Started implementing the encounter saving with user inputs --- src/modules/symba_classes.f90 | 2 +- src/setup/setup.f90 | 1 - src/symba/symba_io.f90 | 15 +++- src/symba/symba_step.f90 | 129 +++++++++++++++++----------------- 4 files changed, 81 insertions(+), 66 deletions(-) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 871637c3d..c715e70da 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -222,7 +222,7 @@ module symba_classes procedure :: dealloc => symba_util_dealloc_system !! Deallocates all allocatable arrays procedure :: resize_storage => symba_util_resize_storage !! Resizes the encounter history storage object so that it contains enough spaces for the number of snapshots needed procedure :: snapshot => symba_util_take_encounter_snapshot !! Take a minimal snapshot of the system through an encounter - procedure :: start_encounter => symba_io_start_encounter !! Initializes the new encounter and/or fragmentation save file(s) + procedure :: start_encounter => symba_io_start_encounter !! Initializes the new encounter history procedure :: stop_encounter => symba_io_stop_encounter !! Saves the encounter and/or fragmentation data to file(s) final :: symba_util_final_system !! Finalizes the SyMBA nbody system object - deallocates all allocatables end type symba_nbody_system diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index ea6d3db23..e95505b9b 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -68,7 +68,6 @@ module subroutine setup_construct_system(system, param) allocate(symba_pltpenc :: system%pltpenc_list) allocate(symba_plplenc :: system%plplenc_list) allocate(symba_plplenc :: system%plplcollision_list) - allocate(symba_encounter_storage :: system%encounter_history) end select case (RINGMOONS) write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled' diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index f531d992e..0300c6dc5 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -356,19 +356,27 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms write(*,*) "Error writing parameter file for SyMBA: " // trim(adjustl(iomsg)) end subroutine symba_io_param_writer + module subroutine symba_io_start_encounter(self, param, t) !! author: David A. Minton !! - !! Initializes the new encounter and/or fragmentation save file(s) + !! Initializes the new encounter and/or fragmentation history implicit none ! Arguments class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object class(symba_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time + if (allocated(self%encounter_history)) deallocate(self%encounter_history) + allocate(symba_encounter_storage :: self%encounter_history) + + ! Take the snapshot at the start of the encounter + call self%snapshot(param, t) + return end subroutine symba_io_start_encounter + module subroutine symba_io_stop_encounter(self, param, t) !! author: David A. Minton !! @@ -379,6 +387,11 @@ module subroutine symba_io_stop_encounter(self, param, t) class(symba_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time + ! Take the final snapshot + call self%snapshot(param, t) + call self%encounter_history%dump(param) + deallocate(self%encounter_history) + return end subroutine symba_io_stop_encounter diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 7b26f80af..bdf73309b 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -181,71 +181,74 @@ recursive module subroutine symba_step_recur_system(self, param, t, ireci) logical :: lencounter, lplpl_collision, lpltp_collision associate(system => self, plplenc_list => self%plplenc_list, pltpenc_list => self%pltpenc_list) - select type(pl => self%pl) - class is (symba_pl) - select type(tp => self%tp) - class is (symba_tp) - system%irec = ireci - dtl = param%dt / (NTENC**ireci) - dth = 0.5_DP * dtl - IF (dtl / param%dt < VSMALL) THEN - write(*, *) "SWIFTEST Warning:" - write(*, *) " In symba_step_recur_system, local time step is too small" - write(*, *) " Roundoff error will be important!" - call util_exit(FAILURE) - END IF - irecp = ireci + 1 - if (ireci == 0) then - nloops = 1 - else - nloops = NTENC - end if - do j = 1, nloops - lencounter = plplenc_list%encounter_check(param, system, dtl, irecp) & - .or. pltpenc_list%encounter_check(param, system, dtl, irecp) - - call plplenc_list%kick(system, dth, irecp, 1) - call pltpenc_list%kick(system, dth, irecp, 1) - if (ireci /= 0) then - call plplenc_list%kick(system, dth, irecp, -1) - call pltpenc_list%kick(system, dth, irecp, -1) - end if - - if (param%lgr) then - call pl%gr_pos_kick(system, param, dth) - call tp%gr_pos_kick(system, param, dth) - end if - - call pl%drift(system, param, dtl) - call tp%drift(system, param, dtl) - - if (lencounter) call system%recursive_step(param, t+dth,irecp) + select type(param) + class is (symba_parameters) + select type(pl => self%pl) + class is (symba_pl) + select type(tp => self%tp) + class is (symba_tp) system%irec = ireci - - if (param%lgr) then - call pl%gr_pos_kick(system, param, dth) - call tp%gr_pos_kick(system, param, dth) + dtl = param%dt / (NTENC**ireci) + dth = 0.5_DP * dtl + IF (dtl / param%dt < VSMALL) THEN + write(*, *) "SWIFTEST Warning:" + write(*, *) " In symba_step_recur_system, local time step is too small" + write(*, *) " Roundoff error will be important!" + call util_exit(FAILURE) + END IF + irecp = ireci + 1 + if (ireci == 0) then + nloops = 1 + else + nloops = NTENC end if - - call plplenc_list%kick(system, dth, irecp, 1) - call pltpenc_list%kick(system, dth, irecp, 1) - if (ireci /= 0) then - call plplenc_list%kick(system, dth, irecp, -1) - call pltpenc_list%kick(system, dth, irecp, -1) - end if - - if (param%lclose) then - lplpl_collision = plplenc_list%collision_check(system, param, t+dtl, dtl, ireci) - lpltp_collision = pltpenc_list%collision_check(system, param, t+dtl, dtl, ireci) - - if (lplpl_collision) call plplenc_list%resolve_collision(system, param, t+dtl, dtl, ireci) - if (lpltp_collision) call pltpenc_list%resolve_collision(system, param, t+dtl, dtl, ireci) - end if - call system%snapshot(param, t+dtl) - - call self%set_recur_levels(ireci) - - end do + do j = 1, nloops + lencounter = plplenc_list%encounter_check(param, system, dtl, irecp) & + .or. pltpenc_list%encounter_check(param, system, dtl, irecp) + + call plplenc_list%kick(system, dth, irecp, 1) + call pltpenc_list%kick(system, dth, irecp, 1) + if (ireci /= 0) then + call plplenc_list%kick(system, dth, irecp, -1) + call pltpenc_list%kick(system, dth, irecp, -1) + end if + + if (param%lgr) then + call pl%gr_pos_kick(system, param, dth) + call tp%gr_pos_kick(system, param, dth) + end if + + call pl%drift(system, param, dtl) + call tp%drift(system, param, dtl) + + if (lencounter) call system%recursive_step(param, t+dth,irecp) + system%irec = ireci + + if (param%lgr) then + call pl%gr_pos_kick(system, param, dth) + call tp%gr_pos_kick(system, param, dth) + end if + + call plplenc_list%kick(system, dth, irecp, 1) + call pltpenc_list%kick(system, dth, irecp, 1) + if (ireci /= 0) then + call plplenc_list%kick(system, dth, irecp, -1) + call pltpenc_list%kick(system, dth, irecp, -1) + end if + + if (param%lclose) then + lplpl_collision = plplenc_list%collision_check(system, param, t+dtl, dtl, ireci) + lpltp_collision = pltpenc_list%collision_check(system, param, t+dtl, dtl, ireci) + + if (lplpl_collision) call plplenc_list%resolve_collision(system, param, t+dtl, dtl, ireci) + if (lpltp_collision) call pltpenc_list%resolve_collision(system, param, t+dtl, dtl, ireci) + end if + if (param%lencounter_save) call system%snapshot(param, t+dtl) + + call self%set_recur_levels(ireci) + + end do + end select end select end select end associate