From a237141cf85df5af5078bf4259c8b7335cbd4f03 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 8 Jan 2023 17:12:04 -0500 Subject: [PATCH] Major update aimed at improving the stability and memory usage of the code. Trying to beat down out of memory errors. --- src/base/base_module.f90 | 245 ++++++++--- src/collision/collision_generate.f90 | 2 +- src/collision/collision_io.f90 | 6 +- src/collision/collision_module.f90 | 109 ++--- src/collision/collision_resolve.f90 | 4 +- src/collision/collision_util.f90 | 181 +++++---- src/encounter/encounter_io.f90 | 6 +- src/encounter/encounter_module.f90 | 90 ++-- src/encounter/encounter_util.f90 | 125 +++--- src/fraggle/fraggle_generate.f90 | 11 +- src/fraggle/fraggle_module.f90 | 67 +-- src/fraggle/fraggle_util.f90 | 54 +-- src/globals/globals_module.f90 | 2 +- src/helio/helio_module.f90 | 46 --- src/misc/solver_module.f90 | 2 +- src/netcdf_io/netcdf_io_implementations.f90 | 2 +- src/operator/operator_cross.f90 | 9 + src/operator/operator_mag.f90 | 8 +- src/operator/operator_unit.f90 | 6 +- src/rmvs/rmvs_module.f90 | 24 +- src/rmvs/rmvs_step.f90 | 8 +- src/rmvs/rmvs_util.f90 | 28 +- src/swiftest/swiftest_driver.f90 | 2 +- src/swiftest/swiftest_io.f90 | 52 +-- src/swiftest/swiftest_module.f90 | 200 +++++---- src/swiftest/swiftest_util.f90 | 429 +++++++++++++++----- src/symba/symba_encounter_check.f90 | 8 +- src/symba/symba_kick.f90 | 2 +- src/symba/symba_module.f90 | 98 +---- src/symba/symba_step.f90 | 2 +- src/symba/symba_util.f90 | 35 +- src/whm/whm_drift.f90 | 2 +- src/whm/whm_module.f90 | 4 +- src/whm/whm_util.f90 | 6 +- 34 files changed, 1034 insertions(+), 841 deletions(-) diff --git a/src/base/base_module.f90 b/src/base/base_module.f90 index e3f818d7e..85fe74169 100644 --- a/src/base/base_module.f90 +++ b/src/base/base_module.f90 @@ -107,10 +107,11 @@ module base logical :: lyarkovsky = .false. !! Turn on Yarkovsky effect logical :: lyorp = .false. !! Turn on YORP effect contains - procedure(abstract_io_dump_param), deferred :: dump - procedure(abstract_io_param_reader), deferred :: reader - procedure(abstract_io_param_writer), deferred :: writer - procedure(abstract_io_read_in_param), deferred :: read_in + procedure :: dealloc => base_util_dealloc_param + procedure(abstract_io_dump_param), deferred :: dump + procedure(abstract_io_param_reader), deferred :: reader + procedure(abstract_io_param_writer), deferred :: writer + procedure(abstract_io_read_in_param), deferred :: read_in end type base_parameters abstract interface @@ -157,27 +158,31 @@ end subroutine abstract_io_read_in_param type :: base_storage_frame class(*), allocatable :: item contains - procedure :: store => copy_store !! Stores a snapshot of the nbody system so that later it can be retrieved for saving to file. + procedure :: store => base_util_copy_store !! Stores a snapshot of the nbody system so that later it can be retrieved for saving to file. generic :: assignment(=) => store - final :: final_storage_frame + final :: base_final_storage_frame end type - type, abstract :: base_storage(nframes) + type, abstract :: base_storage !! An class that establishes the pattern for various storage objects - integer(I4B), len :: nframes = 4096 !! Total number of frames that can be stored + integer(I4B) :: nframes !! Total number of frames that can be stored !! An class that establishes the pattern for various storage objects - type(base_storage_frame), dimension(nframes) :: frame !! Array of stored frames - integer(I4B) :: iframe = 0 !! Index of the last frame stored in the system - integer(I4B) :: nid !! Number of unique id values in all saved snapshots - integer(I4B), dimension(:), allocatable :: idvals !! The set of unique id values contained in the snapshots - integer(I4B), dimension(:), allocatable :: idmap !! The id value -> index map - integer(I4B) :: nt !! Number of unique time values in all saved snapshots - real(DP), dimension(:), allocatable :: tvals !! The set of unique time values contained in the snapshots - integer(I4B), dimension(:), allocatable :: tmap !! The t value -> index map + type(base_storage_frame), dimension(:), allocatable :: frame !! Array of stored frames + integer(I4B) :: iframe = 0 !! Index of the last frame stored in the system + integer(I4B) :: nid !! Number of unique id values in all saved snapshots + integer(I4B), dimension(:), allocatable :: idvals !! The set of unique id values contained in the snapshots + integer(I4B), dimension(:), allocatable :: idmap !! The id value -> index map + integer(I4B) :: nt !! Number of unique time values in all saved snapshots + real(DP), dimension(:), allocatable :: tvals !! The set of unique time values contained in the snapshots + integer(I4B), dimension(:), allocatable :: tmap !! The t value -> index map contains - procedure :: reset => reset_storage !! Resets a storage object by deallocating all items and resetting the frame counter to 0 + procedure :: dealloc => base_util_dealloc_storage !! Deallocates all allocatables + procedure :: reset => base_util_reset_storage !! Resets the storage object back to its original state by removing all of the saved items from the storage frames + procedure :: resize => base_util_resize_storage !! Resizes storage if it is too small + procedure :: setup => base_util_setup_storage !! Sets up a storage system with a set number of frames + procedure :: save => base_util_snapshot_save !! Takes a snapshot of the current system end type base_storage @@ -189,8 +194,18 @@ end subroutine abstract_io_read_in_param !> An abstract class for a generic collection of Swiftest bodies type, abstract :: base_object + contains + procedure(abstract_util_dealloc_object), deferred :: dealloc end type base_object + abstract interface + subroutine abstract_util_dealloc_object(self) + import base_object + implicit none + class(base_object), intent(inout) :: self !! Generic Swiftest object type + end subroutine abstract_util_dealloc_object + end interface + type, abstract :: base_multibody(nbody) integer(I4B), len :: nbody @@ -209,7 +224,7 @@ end subroutine abstract_io_read_in_param contains - subroutine copy_store(self, source) + subroutine base_util_copy_store(self, source) !! author: David A. Minton !! !! Stores a snapshot of the nbody system so that later it can be retrieved for saving to file. @@ -221,70 +236,49 @@ subroutine copy_store(self, source) allocate(self%item, source=source) return - end subroutine copy_store + end subroutine base_util_copy_store - subroutine final_storage_frame(self) + subroutine base_util_dealloc_param(self) !! author: David A. Minton !! - !! Finalizer for the storage frame data type - implicit none - type(base_storage_frame) :: self - - if (allocated(self%item)) deallocate(self%item) - - return - end subroutine final_storage_frame - - - subroutine base_final_storage(self) - !! author: David A. Minton - !! - !! Finalizer for the storage object + !! Deallocates all allocatables implicit none ! Arguments - class(base_storage(*)), intent(inout) :: self - ! Internals - integer(I4B) :: i + class(base_parameters),intent(inout) :: self !! Collection of parameters + + if (allocated(self%integrator)) deallocate(self%integrator) + if (allocated(self%param_file_name)) deallocate(self%param_file_name) + if (allocated(self%display_style)) deallocate(self%display_style) + if (allocated(self%seed)) deallocate(self%seed) - do i = 1, self%nframes - call final_storage_frame(self%frame(i)) - end do return - end subroutine base_final_storage + end subroutine base_util_dealloc_param - subroutine reset_storage(self) + subroutine base_util_dealloc_storage(self) !! author: David A. Minton !! !! Resets a storage object by deallocating all items and resetting the frame counter to 0 implicit none ! Arguments - class(base_storage(*)), intent(inout) :: self !! Swiftest storage object - ! Internals - integer(I4B) :: i - - do i = 1, self%nframes - if (allocated(self%frame(i)%item)) deallocate(self%frame(i)%item) - end do - - if (allocated(self%idmap)) deallocate(self%idmap) - if (allocated(self%tmap)) deallocate(self%tmap) - self%nid = 0 - self%nt = 0 - self%iframe = 0 + class(base_storage), intent(inout) :: self !! Swiftest storage object + + call self%reset() + if (allocated(self%frame)) deallocate(self%frame) + self%nframes = 0 return - end subroutine reset_storage + end subroutine base_util_dealloc_storage - subroutine util_exit(code) + subroutine base_util_exit(code) !! author: David A. Minton !! !! Print termination message and exit program !! - !! Adapted from David E. Kaufmann's Swifter routine: util_exit.f90 - !! Adapted from Hal Levison's Swift routine util_exit.f + !! Adapted from David E. Kaufmann's Swifter routine: base_util_exit.f90 + !! Adapted from Hal Levison's Swift routine base_util_exit.f implicit none ! Arguments integer(I4B), intent(in) :: code @@ -311,6 +305,137 @@ subroutine util_exit(code) stop - end subroutine util_exit + end subroutine base_util_exit + + + subroutine base_util_reset_storage(self) + !! author: David A. Minton + !! + !! Resets the storage object back to its original state by removing all of the saved items from the storage frames, but does not deallocate the frames + implicit none + ! Arguments + class(base_storage), intent(inout) :: self + ! Internals + integer(I4B) :: i + + if (allocated(self%frame)) then + do i = 1, self%nframes + if (allocated(self%frame(i)%item)) deallocate(self%frame(i)%item) + end do + end if + + if (allocated(self%idmap)) deallocate(self%idmap) + if (allocated(self%idvals)) deallocate(self%idvals) + if (allocated(self%tmap)) deallocate(self%tmap) + if (allocated(self%tvals)) deallocate(self%tvals) + self%nid = 0 + self%nt = 0 + self%iframe = 0 + + return + end subroutine base_util_reset_storage + + + subroutine base_util_resize_storage(self, nnew) + !! author: David A. Minton + !! + !! Checks the current size of a Swiftest against the requested size and resizes it if it is too small. + implicit none + ! Arguments + class(base_storage), intent(inout) :: self !! Storage object + integer(I4B), intent(in) :: nnew !! New size + ! Internals + class(base_storage_frame), dimension(:), allocatable :: tmp + integer(I4B) :: i, n, nold, nbig + + nold = self%nframes + if (nnew <= nold) return + + nbig = nold + do while (nbig < nnew) + nbig = nbig * 2 + end do + call move_alloc(self%frame, tmp) + allocate(self%frame(nbig)) + self%nframes = nbig + do i = 1, nold + if (allocated(tmp(i)%item)) call move_alloc(tmp(i)%item, self%frame(i)%item) + end do + + return + end subroutine base_util_resize_storage + + + subroutine base_util_setup_storage(self, n) + !! author: David A. Minton + !! + !! Checks the current size of a Swiftest against the requested size and resizes it if it is too small. + implicit none + ! Arguments + class(base_storage), intent(inout) :: self !! Storage object + integer(I4B), intent(in) :: n !! New size + + if (allocated(self%frame)) deallocate(self%frame) + allocate(self%frame(n)) + self%nframes = n + + return + end subroutine base_util_setup_storage + + + subroutine base_util_snapshot_save(self, snapshot) + !! author: David A. Minton + !! + !! Checks the current size of the storage object against the required size and extends it by a factor of 2 more than requested if it is too small. + !! Note: The reason to extend it by a factor of 2 is for performance. When there are many enounters per step, resizing every time you want to add an + !! encounter takes significant computational effort. Resizing by a factor of 2 is a tradeoff between performance (fewer resize calls) and memory managment + !! Memory usage grows by a factor of 2 each time it fills up, but no more. + implicit none + ! Arguments + class(base_storage), intent(inout) :: self !! Storage ncounter storage object + class(*), intent(in) :: snapshot !! Object to snapshot + ! Internals + integer(I4B) :: nnew, nold + + ! Advance the snapshot frame counter + self%iframe = self%iframe + 1 + + nnew = self%iframe + nold = self%nframes + + ! Check to make sure the current storage object is big enough. If not, grow it by a factor of 2 + if (nnew > nold) call self%resize(nnew) + + self%frame(nnew) = snapshot + + return + end subroutine base_util_snapshot_save + + + subroutine base_final_storage(self) + !! author: David A. Minton + !! + !! Finalizer for the storage object + implicit none + ! Arguments + class(base_storage), intent(inout) :: self + + call self%dealloc() + return + end subroutine base_final_storage + + + subroutine base_final_storage_frame(self) + !! author: David A. Minton + !! + !! Finalizer for the storage frame data type + implicit none + type(base_storage_frame) :: self + + if (allocated(self%item)) deallocate(self%item) + + return + end subroutine base_final_storage_frame + end module base diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index cc0363b84..288c00ebd 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -87,7 +87,7 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t) call self%merge(nbody_system, param, t) ! Use the default collision model, which is merge case default call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Error in swiftest_collision, unrecognized collision regime") - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select end associate end select diff --git a/src/collision/collision_io.f90 b/src/collision/collision_io.f90 index 1fc3e39c1..476700a52 100644 --- a/src/collision/collision_io.f90 +++ b/src/collision/collision_io.f90 @@ -90,7 +90,7 @@ module subroutine collision_io_netcdf_dump(self, param) !! Dumps the time history of an encounter to file. implicit none ! Arguments - class(collision_storage(*)), intent(inout) :: self !! Encounter storage object + class(collision_storage), intent(inout) :: self !! Encounter storage object class(base_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i @@ -266,7 +266,7 @@ module subroutine collision_io_netcdf_initialize_output(self, param) 667 continue write(*,*) "Error creating fragmentation output file. " // trim(adjustl(errmsg)) - call util_exit(FAILURE) + call base_util_exit(FAILURE) end subroutine collision_io_netcdf_initialize_output @@ -359,7 +359,7 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) implicit none ! Arguments class(collision_snapshot), intent(in) :: self !! Swiftest encounter structure - class(encounter_storage(*)), intent(inout) :: history !! Collision history object + class(encounter_storage), intent(inout) :: history !! Collision history object class(base_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i, idslot, old_mode, npl, stage diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 5b6e041f5..3ce48a518 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -82,7 +82,7 @@ module collision contains procedure :: consolidate => collision_resolve_consolidate_impactors !! Consolidates a multi-body collision into an equivalent 2-body collision procedure :: get_regime => collision_regime_impactors !! Determine which fragmentation regime the set of impactors will be - procedure :: reset => collision_util_reset_impactors !! Resets the collider object variables to 0 and deallocates the index and mass distributions + procedure :: dealloc => collision_util_dealloc_impactors !! Resets the collider object variables to 0 and deallocates the index and mass distributions procedure :: set_coordinate_system => collision_util_set_coordinate_impactors !! Sets the coordinate system of the impactors final :: collision_final_impactors !! Finalizer will deallocate all allocatables end type collision_impactors @@ -124,9 +124,10 @@ module collision real(DP), dimension(nbody) :: ke_orbit !! Orbital kinetic energy of each individual fragment real(DP), dimension(nbody) :: ke_spin !! Spin kinetic energy of each individual fragment contains - procedure :: reset => collision_util_reset_fragments !! Deallocates all allocatable arrays and sets everything else to 0 + procedure :: dealloc => collision_util_dealloc_fragments !! Deallocates all allocatable arrays and sets everything else to 0 + procedure :: reset => collision_util_reset_fragments !! Resets all position and velocity-dependent fragment quantities in order to do a fresh calculation (does not reset mass, radius, or other values that get set prior to the call to fraggle_generate) procedure :: set_coordinate_system => collision_util_set_coordinate_fragments !! Sets the coordinate system of the fragments - final :: collision_final_fragments !! Finalizer deallocates all allocatables + final :: collision_final_fragments !! Finalizer deallocates all allocatables end type collision_fragments @@ -166,20 +167,20 @@ module collision procedure :: merge => collision_generate_merge !! Merges the impactors to make a single final body procedure :: add_fragments => collision_util_add_fragments_to_collider !! Add fragments to nbody_system procedure :: get_energy_and_momentum => collision_util_get_energy_and_momentum !! Calculates total nbody_system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) - procedure :: reset => collision_util_reset_system !! Deallocates all allocatables + procedure :: dealloc => collision_util_dealloc_basic !! Deallocates all allocatables procedure :: setup => collision_util_setup_collider !! Initializer for the encounter collision system and the before/after snapshots procedure :: setup_impactors => collision_util_setup_impactors_collider !! Initializer for the impactors for the encounter collision system. Deallocates old impactors before creating new ones procedure :: setup_fragments => collision_util_setup_fragments_collider !! Initializer for the fragments of the collision system. procedure :: set_coordinate_system => collision_util_set_coordinate_collider !! Sets the coordinate system of the collisional system procedure :: set_natural_scale => collision_util_set_natural_scale_factors !! Scales dimenional quantities to ~O(1) with respect to the collisional system. procedure :: set_original_scale => collision_util_set_original_scale_factors !! Restores dimenional quantities back to the original system units + final :: collision_final_basic end type collision_basic type, extends(collision_basic) :: collision_bounce contains procedure :: generate => collision_generate_bounce !! If a collision would result in a disruption, "bounce" the bodies instead. - final :: collision_final_bounce !! Finalizer will deallocate all allocatables end type collision_bounce @@ -199,17 +200,16 @@ module collision contains procedure :: initialize => collision_io_netcdf_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object procedure :: open => collision_io_netcdf_open !! Opens an old file - final :: collision_final_netcdf_parameters !! Finalizer closes the NetCDF file end type collision_netcdf_parameters type, extends(encounter_snapshot) :: collision_snapshot - logical :: lcollision !! Indicates that this snapshot contains at least one collision + logical :: lcollision !! Indicates that this snapshot contains at least one collision class(collision_basic), allocatable :: collider !! Collider object at this snapshot contains procedure :: write_frame => collision_io_netcdf_write_frame_snapshot !! Writes a frame of encounter data to file + procedure :: dealloc => collision_util_dealloc_snapshot !! Deallocates all allocatables procedure :: get_idvals => collision_util_get_idvalues_snapshot !! Gets an array of all id values saved in this snapshot - final :: collision_final_snapshot !! Finalizer deallocates all allocatables end type collision_snapshot @@ -219,7 +219,6 @@ module collision procedure :: dump => collision_io_netcdf_dump !! Dumps contents of encounter history to file procedure :: take_snapshot => collision_util_snapshot !! Take a minimal snapshot of the nbody_system through an encounter procedure :: make_index_map => collision_util_index_map !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id - final :: collision_final_storage !! Finalizer deallocates all allocatables end type collision_storage @@ -270,8 +269,8 @@ end subroutine collision_io_log_regime module subroutine collision_io_netcdf_dump(self, param) implicit none - class(collision_storage(*)), intent(inout) :: self !! Collision storage object - class(base_parameters), intent(inout) :: param !! Current run configuration parameters + class(collision_storage), intent(inout) :: self !! Collision storage object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine collision_io_netcdf_dump module subroutine collision_io_netcdf_initialize_output(self, param) @@ -290,7 +289,7 @@ end subroutine collision_io_netcdf_open module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) implicit none class(collision_snapshot), intent(in) :: self !! Swiftest encounter structure - class(encounter_storage(*)), intent(inout) :: history !! Collision history object + class(encounter_storage), intent(inout) :: history !! Collision history object class(base_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine collision_io_netcdf_write_frame_snapshot @@ -383,7 +382,7 @@ module subroutine collision_util_add_fragments_to_collider(self, nbody_system, p implicit none class(collision_basic), intent(in) :: self !! Collision system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters end subroutine collision_util_add_fragments_to_collider module subroutine collision_util_construct_after_system(collider, nbody_system, param, after_system) @@ -393,11 +392,21 @@ module subroutine collision_util_construct_after_system(collider, nbody_system, implicit none ! Arguments class(collision_basic), intent(inout) :: collider !! Collision system object - class(base_nbody_system), intent(in) :: nbody_system !! Original swiftest nbody system object - class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters - class(base_nbody_system), allocatable, intent(out) :: after_system !! Output temporary swiftest nbody system object + class(base_nbody_system), intent(in) :: nbody_system !! Original Swiftest nbody system object + class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters + class(base_nbody_system), allocatable, intent(out) :: after_system !! Output temporary Swiftest nbody system object end subroutine collision_util_construct_after_system + module subroutine collision_util_dealloc_fragments(self) + implicit none + class(collision_fragments(*)), intent(inout) :: self + end subroutine collision_util_dealloc_fragments + + module subroutine collision_util_dealloc_snapshot(self) + implicit none + class(collision_snapshot), intent(inout) :: self !! Collsion snapshot object + end subroutine collision_util_dealloc_snapshot + module subroutine collision_util_reset_fragments(self) implicit none class(collision_fragments(*)), intent(inout) :: self @@ -452,28 +461,28 @@ module subroutine collision_util_get_energy_and_momentum(self, nbody_system, par implicit none class(collision_basic), intent(inout) :: self !! Encounter collision system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(base_parameters), intent(inout) :: param !! Current swiftest run configuration parameters + class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters character(len=*), intent(in) :: phase !! One of "before" or "after", indicating which phase of the calculation this needs to be done end subroutine collision_util_get_energy_and_momentum module subroutine collision_util_index_map(self) implicit none - class(collision_storage(*)), intent(inout) :: self !! Collision storage object + class(collision_storage), intent(inout) :: self !! Collision storage object end subroutine collision_util_index_map - module subroutine collision_util_reset_impactors(self) + module subroutine collision_util_dealloc_impactors(self) implicit none class(collision_impactors), intent(inout) :: self !! Collision system object - end subroutine collision_util_reset_impactors + end subroutine collision_util_dealloc_impactors - module subroutine collision_util_reset_system(self) + module subroutine collision_util_dealloc_basic(self) implicit none class(collision_basic), intent(inout) :: self !! Collision system object - end subroutine collision_util_reset_system + end subroutine collision_util_dealloc_basic module subroutine collision_util_snapshot(self, param, nbody_system, t, arg) implicit none - class(collision_storage(*)), intent(inout) :: self !! Swiftest storage object + class(collision_storage), intent(inout) :: self !! Swiftest storage object class(base_parameters), intent(inout) :: param !! Current run configuration parameters class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object to store real(DP), intent(in), optional :: t !! Time of snapshot if different from nbody_system time @@ -514,24 +523,11 @@ subroutine collision_final_impactors(self) ! Arguments type(collision_impactors), intent(inout) :: self !! Collision impactors storage object - call self%reset() + call self%dealloc() return end subroutine collision_final_impactors - subroutine collision_final_netcdf_parameters(self) - !! author: David A. Minton - !! - !! Finalize the NetCDF by closing the file - implicit none - ! Arguments - type(collision_netcdf_parameters), intent(inout) :: self - - call self%close() - - return - end subroutine collision_final_netcdf_parameters - subroutine collision_final_plpl(self) !! author: David A. Minton !! @@ -558,50 +554,19 @@ subroutine collision_final_pltp(self) return end subroutine collision_final_pltp - subroutine collision_final_snapshot(self) + subroutine collision_final_basic(self) !! author: David A. Minton !! !! Finalizer will deallocate all allocatables implicit none ! Arguments - type(collision_snapshot), intent(inout) :: self !! Collsion snapshot object + type(collision_basic), intent(inout) :: self !! Collision system object - call encounter_final_snapshot(self%encounter_snapshot) - - return - end subroutine collision_final_snapshot - - subroutine collision_final_storage(self) - !! author: David A. Minton - !! - !! Finalizer will deallocate all allocatables - implicit none - ! Arguments - type(collision_storage(*)), intent(inout) :: self !! Collision storage object - ! Internals - integer(I4B) :: i - - do i = 1, self%nframes - if (allocated(self%frame(i)%item)) deallocate(self%frame(i)%item) - end do + call self%dealloc() return - end subroutine collision_final_storage + end subroutine collision_final_basic - subroutine collision_final_bounce(self) - !! author: David A. Minton - !! - !! Finalizer will deallocate all allocatables - implicit none - ! Arguments - type(collision_bounce), intent(inout) :: self !! Collision system object - - call self%reset() - if (allocated(self%impactors)) deallocate(self%impactors) - if (allocated(self%fragments)) deallocate(self%fragments) - - return - end subroutine collision_final_bounce end module collision diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index 834f15de0..4b88971cb 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -556,7 +556,7 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) call collision_history%take_snapshot(param,nbody_system, t, "after") plpl_collision%status(i) = collider%status - call impactors%reset() + call impactors%dealloc() end do ! Destroy the collision list now that the collisions are resolved @@ -581,7 +581,7 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) if (loop == MAXCASCADE) then call swiftest_io_log_one_message(COLLISION_LOG_OUT,"An runaway collisional cascade has been detected in collision_resolve_plpl.") call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Consider reducing the step size or changing the parameters in the collisional model to reduce the number of fragments.") - call util_exit(FAILURE) + call base_util_exit(FAILURE) end if end associate end do diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 5848f35ff..d1598c62f 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -19,7 +19,7 @@ module subroutine collision_util_add_fragments_to_collider(self, nbody_system, p ! Arguments class(collision_basic), intent(in) :: self !! Collision system system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters ! Internals integer(I4B) :: i, npl_before, npl_after logical, dimension(:), allocatable :: lexclude @@ -69,13 +69,13 @@ module subroutine collision_util_construct_after_system(collider, nbody_system, implicit none ! Arguments class(collision_basic), intent(inout) :: collider !! Collision system object - class(base_nbody_system), intent(in) :: nbody_system !! Original swiftest nbody system object - class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters - class(base_nbody_system), allocatable, intent(out) :: after_system !! Output temporary swiftest nbody system object + class(base_nbody_system), intent(in) :: nbody_system !! Original Swiftest nbody system object + class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters + class(base_nbody_system), allocatable, intent(out) :: after_system !! Output temporary Swiftest nbody system object ! Internals integer(I4B) :: i, status class(swiftest_nbody_system), allocatable :: tmpsys - class(swiftest_parameters), allocatable :: tmpparam + !class(swiftest_parameters), allocatable :: tmpparam select type(nbody_system) class is (swiftest_nbody_system) @@ -90,16 +90,11 @@ module subroutine collision_util_construct_after_system(collider, nbody_system, end do ! Set up a new temporary system based on the original - allocate(tmpparam, source=param) - call swiftest_util_setup_construct_system(tmpsys, tmpparam) - - if (allocated(tmpsys%cb)) deallocate(tmpsys%cb) + allocate(tmpsys, mold=nbody_system) allocate(tmpsys%cb, source=cb) - - if (allocated(tmpsys%pl)) deallocate(tmpsys%pl) allocate(tmpsys%pl, source=pl) - - if (allocated(tmpsys%collider)) deallocate(tmpsys%collider) + allocate(tmpsys%pl_adds, mold=pl) + allocate(tmpsys%pl_discards, mold=pl) allocate(tmpsys%collider, source=collider) call tmpsys%collider%set_original_scale() @@ -112,8 +107,8 @@ module subroutine collision_util_construct_after_system(collider, nbody_system, status = HIT_AND_RUN_DISRUPT end select - call collision_resolve_mergeaddsub(tmpsys, tmpparam, nbody_system%t, status) - call tmpsys%pl%rearray(tmpsys, tmpparam) + call collision_resolve_mergeaddsub(tmpsys, param, nbody_system%t, status) + call tmpsys%pl%rearray(tmpsys, param) call move_alloc(tmpsys, after_system) @@ -189,7 +184,7 @@ module subroutine collision_util_get_energy_and_momentum(self, nbody_system, par ! Arguments class(collision_basic), intent(inout) :: self !! Encounter collision system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(base_parameters), intent(inout) :: param !! Current swiftest run configuration parameters + class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters character(len=*), intent(in) :: phase !! One of "before" or "after", indicating which phase of the calculation this needs to be done ! Internals class(base_nbody_system), allocatable :: after_system @@ -253,7 +248,7 @@ module subroutine collision_util_index_map(self) !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id implicit none ! Arguments - class(collision_storage(*)), intent(inout) :: self !! Swiftest storage object + class(collision_storage), intent(inout) :: self !! Swiftest storage object ! Internals integer(I4B), dimension(:), allocatable :: idvals real(DP), dimension(:), allocatable :: tvals @@ -271,7 +266,25 @@ module subroutine collision_util_index_map(self) end subroutine collision_util_index_map - module subroutine collision_util_reset_impactors(self) + module subroutine collision_util_dealloc_snapshot(self) + !! author: David A. Minton + !! + !! Finalizer will deallocate all allocatables + implicit none + ! Arguments + class(collision_snapshot), intent(inout) :: self !! Collsion snapshot object + + if (allocated(self%collider)) then + call self%collider%dealloc() + deallocate(self%collider) + end if + + call self%encounter_snapshot%dealloc() + + return + end subroutine collision_util_dealloc_snapshot + + module subroutine collision_util_dealloc_impactors(self) !! author: David A. Minton !! !! Resets the collider object variables to 0 and deallocates the index and mass distributions @@ -302,10 +315,10 @@ module subroutine collision_util_reset_impactors(self) self%rbimp(:) = 0.0_DP return - end subroutine collision_util_reset_impactors + end subroutine collision_util_dealloc_impactors - module subroutine collision_util_reset_fragments(self) + module subroutine collision_util_dealloc_fragments(self) !! author: David A. Minton !! !! Deallocates all allocatables @@ -332,16 +345,26 @@ module subroutine collision_util_reset_fragments(self) self%n_unit(:,:) = 0.0_DP return - end subroutine collision_util_reset_fragments + end subroutine collision_util_dealloc_fragments - module subroutine collision_util_reset_system(self) + module subroutine collision_util_dealloc_basic(self) !! author: David A. Minton !! !! Resets the collider nbody_system and deallocates all allocatables implicit none ! Arguments - class(collision_basic), intent(inout) :: self !! Collision system object + class(collision_basic), intent(inout) :: self !! Collision system object + + if (allocated(self%impactors)) then + call self%impactors%dealloc() + deallocate(self%impactors) + end if + + if (allocated(self%fragments)) then + call self%fragments%dealloc() + deallocate(self%fragments) + end if select type(before => self%before) class is (swiftest_nbody_system) @@ -362,15 +385,57 @@ module subroutine collision_util_reset_system(self) self%pe(:) = 0.0_DP self%te(:) = 0.0_DP - if (allocated(self%impactors)) call self%impactors%reset() - if (allocated(self%fragments)) deallocate(self%fragments) + self%dscale = 1.0_DP + self%mscale = 1.0_DP + self%tscale = 1.0_DP + self%vscale = 1.0_DP + self%Escale = 1.0_DP + self%Lscale = 1.0_DP + return - end subroutine collision_util_reset_system + end subroutine collision_util_dealloc_basic - module subroutine collision_util_set_coordinate_collider(self) + module subroutine collision_util_reset_fragments(self) !! author: David A. Minton + !! + !! Resets all position and velocity-dependent fragment quantities in order to do a fresh calculation (does not reset mass, radius, or other values that get set prior to the call to fraggle_generate) + implicit none + ! Arguments + class(collision_fragments(*)), intent(inout) :: self + + self%rc(:,:) = 0.0_DP + self%vc(:,:) = 0.0_DP + self%rh(:,:) = 0.0_DP + self%vh(:,:) = 0.0_DP + self%rb(:,:) = 0.0_DP + self%vb(:,:) = 0.0_DP + self%rot(:,:) = 0.0_DP + self%r_unit(:,:) = 0.0_DP + self%t_unit(:,:) = 0.0_DP + self%n_unit(:,:) = 0.0_DP + + self%rmag(:) = 0.0_DP + self%vmag(:) = 0.0_DP + self%rotmag(:) = 0.0_DP + + self%L_orbit_tot(:) = 0.0_DP + self%L_spin_tot(:) = 0.0_DP + self%L_orbit(:,:) = 0.0_DP + self%L_spin(:,:) = 0.0_DP + self%ke_orbit_tot = 0.0_DP + self%ke_spin_tot = 0.0_DP + self%pe = 0.0_DP + self%be = 0.0_DP + self%ke_orbit(:) = 0.0_DP + self%ke_spin(:) = 0.0_DP + + return + end subroutine collision_util_reset_fragments + + module subroutine collision_util_set_coordinate_collider(self) + !! !! Defines the collisional coordinate nbody_system, including the unit vectors of both the nbody_system and individual fragments. implicit none @@ -422,7 +487,7 @@ module subroutine collision_util_set_coordinate_impactors(self) !! Defines the collisional coordinate nbody_system, including the unit vectors of both the nbody_system and individual fragments. implicit none ! Arguments - class(collision_impactors), intent(inout) :: self !! Collisional nbody_system + class(collision_impactors), intent(inout) :: self !! Collisional nbody_system ! Internals real(DP), dimension(NDIM) :: delta_r, delta_v, L_total real(DP) :: L_mag, mtot @@ -484,7 +549,7 @@ module subroutine collision_util_setup_collider(self, nbody_system) !! but not fragments. Those are setup later when the number of fragments is known. implicit none ! Arguments - class(collision_basic), intent(inout) :: self !! Encounter collision system object + class(collision_basic), intent(inout) :: self !! Encounter collision system object class(base_nbody_system), intent(in) :: nbody_system !! Current nbody system. Used as a mold for the before/after snapshots call self%setup_impactors() @@ -582,51 +647,6 @@ module subroutine collision_util_shift_vector_to_origin(m_frag, vec_frag) end subroutine collision_util_shift_vector_to_origin - subroutine collision_util_save_snapshot(collision_history, snapshot) - !! author: David A. Minton - !! - !! Checks the current size of the encounter storage against the required size and extends it by a factor of 2 more than requested if it is too small. - !! Note: The reason to extend it by a factor of 2 is for performance. When there are many enounters per step, resizing every time you want to add an - !! encounter takes significant computational effort. Resizing by a factor of 2 is a tradeoff between performance (fewer resize calls) and memory managment - !! Memory usage grows by a factor of 2 each time it fills up, but no more. - implicit none - ! Arguments - class(collision_storage(*)), allocatable, intent(inout) :: collision_history !! Collision history object - class(encounter_snapshot), intent(in) :: snapshot !! Encounter snapshot object - ! Internals - type(collision_storage(nframes=:)), allocatable :: tmp - integer(I4B) :: i, nnew, nold, nbig - - ! Advance the snapshot frame counter - collision_history%iframe = collision_history%iframe + 1 - - ! Check to make sure the current encounter_history object is big enough. If not, grow it by a factor of 2 - nnew = collision_history%iframe - nold = collision_history%nframes - - if (nnew > nold) then - nbig = nold - do while (nbig < nnew) - nbig = nbig * 2 - end do - allocate(collision_storage(nbig) :: tmp) - tmp%iframe = collision_history%iframe - call move_alloc(collision_history%nc, tmp%nc) - - do i = 1, nold - if (allocated(collision_history%frame(i)%item)) call move_alloc(collision_history%frame(i)%item, tmp%frame(i)%item) - end do - deallocate(collision_history) - call move_alloc(tmp,collision_history) - nnew = nbig - end if - - collision_history%frame(nnew) = snapshot - - return - end subroutine collision_util_save_snapshot - - module subroutine collision_util_snapshot(self, param, nbody_system, t, arg) !! author: David A. Minton !! @@ -634,11 +654,11 @@ module subroutine collision_util_snapshot(self, param, nbody_system, t, arg) !! can be played back through the encounter implicit none ! Internals - class(collision_storage(*)), intent(inout) :: self !! Swiftest storage object - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object to store - real(DP), intent(in), optional :: t !! Time of snapshot if different from nbody_system time - character(*), intent(in), optional :: arg !! "before": takes a snapshot just before the collision. "after" takes the snapshot just after the collision. + class(collision_storage), intent(inout) :: self !! Swiftest storage object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters + class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object to store + real(DP), intent(in), optional :: t !! Time of snapshot if different from nbody_system time + character(*), intent(in), optional :: arg !! "before": takes a snapshot just before the collision. "after" takes the snapshot just after the collision. ! Arguments class(collision_snapshot), allocatable, save :: snapshot character(len=:), allocatable :: stage @@ -700,8 +720,7 @@ module subroutine collision_util_snapshot(self, param, nbody_system, t, arg) end select ! Save the snapshot for posterity - call collision_util_save_snapshot(nbody_system%collision_history,snapshot) - deallocate(snapshot) + call self%save(snapshot) case default write(*,*) "collision_util_snapshot requies either 'before' or 'after' passed to 'arg'" end select diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index 7c76cd4ca..f62ced49a 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -17,7 +17,7 @@ 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(encounter_storage), intent(inout) :: self !! Encounter storage object class(base_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i @@ -145,7 +145,7 @@ module subroutine encounter_io_netcdf_initialize_output(self, param) 667 continue write(*,*) "Error creating encounter output file. " // trim(adjustl(errmsg)) - call util_exit(FAILURE) + call base_util_exit(FAILURE) end subroutine encounter_io_netcdf_initialize_output @@ -227,7 +227,7 @@ module subroutine encounter_io_netcdf_write_frame_snapshot(self, history, param) implicit none ! Arguments class(encounter_snapshot), intent(in) :: self !! Swiftest encounter structure - class(encounter_storage(*)), intent(inout) :: history !! Encounter storage object + class(encounter_storage), intent(inout) :: history !! Encounter storage object class(base_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals diff --git a/src/encounter/encounter_module.f90 b/src/encounter/encounter_module.f90 index ed7ae4dad..720627485 100644 --- a/src/encounter/encounter_module.f90 +++ b/src/encounter/encounter_module.f90 @@ -38,24 +38,25 @@ module encounter real(DP), dimension(:,:), allocatable :: v2 !! the velocity of body 2 in the encounter integer(I4B), dimension(:), allocatable :: level !! Recursion level (used in SyMBA) contains - procedure :: setup => encounter_util_setup_list !! A constructor that sets the number of encounters and allocates and initializes all arrays - procedure :: append => encounter_util_append_list !! Appends elements from one structure to another - procedure :: copy => encounter_util_copy_list !! Copies elements from the source encounter list into self. - procedure :: dealloc => encounter_util_dealloc_list !! Deallocates all allocatables - procedure :: spill => encounter_util_spill_list !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) - procedure :: resize => encounter_util_resize_list !! Checks the current size of the encounter list against the required size and extends it by a factor of 2 more than requested if it is too small. + procedure :: setup => encounter_util_setup_list !! A constructor that sets the number of encounters and allocates and initializes all arrays + procedure :: append => encounter_util_append_list !! Appends elements from one structure to another + procedure :: copy => encounter_util_copy_list !! Copies elements from the source encounter list into self. + procedure :: dealloc => encounter_util_dealloc_list !! Deallocates all allocatables + procedure :: spill => encounter_util_spill_list !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) + procedure :: resize => encounter_util_resize_list !! Checks the current size of the encounter list against the required size and extends it by a factor of 2 more than requested if it is too small. end type encounter_list - type :: encounter_snapshot + type, extends(base_object) :: encounter_snapshot !! A simplified version of a SyMBA nbody system object for storing minimal snapshots of the system state during encounters class(base_object), allocatable :: pl !! Massive body data structure class(base_object), allocatable :: tp !! Test particle data structure real(DP) :: t !! Simulation time when snapshot was taken integer(I8B) :: iloop !! Loop number at time of snapshot contains - procedure :: write_frame => encounter_io_netcdf_write_frame_snapshot !! Writes a frame of encounter data to file - procedure :: get_idvals => encounter_util_get_idvalues_snapshot !! Gets an array of all id values saved in this snapshot + procedure :: write_frame => encounter_io_netcdf_write_frame_snapshot !! Writes a frame of encounter data to file + procedure :: dealloc => encounter_util_dealloc_snapshot !! Deallocates all allocatables + procedure :: get_idvals => encounter_util_get_idvalues_snapshot !! Gets an array of all id values saved in this snapshot final :: encounter_final_snapshot end type encounter_snapshot @@ -73,6 +74,7 @@ module encounter class(encounter_netcdf_parameters), allocatable :: nc !! NetCDF object attached to this storage object contains procedure :: dump => encounter_io_netcdf_dump !! Dumps contents of encounter history to file + procedure :: dealloc => encounter_util_dealloc_storage !! Deallocates all allocatables procedure :: get_index_values => encounter_util_get_vals_storage !! Gets the unique values of the indices of a storage object (i.e. body id or time value) procedure :: make_index_map => encounter_util_index_map !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id procedure :: take_snapshot => encounter_util_snapshot !! Take a minimal snapshot of the system through an encounter @@ -95,10 +97,12 @@ module encounter type encounter_bounding_box type(encounter_bounding_box_1D), dimension(SWEEPDIM) :: aabb contains - procedure :: setup => encounter_util_setup_aabb !! Setup a new axis-aligned bounding box structure + procedure :: dealloc => encounter_util_dealloc_bounding_box !! Deallocates all allocatables + procedure :: setup => encounter_util_setup_aabb !! Setup a new axis-aligned bounding box structure procedure :: sweep_single => encounter_check_sweep_aabb_single_list !! Sweeps the sorted bounding box extents and returns the encounter candidates procedure :: sweep_double => encounter_check_sweep_aabb_double_list !! Sweeps the sorted bounding box extents and returns the encounter candidates generic :: sweep => sweep_single, sweep_double + final :: encounter_final_bounding_box end type @@ -122,7 +126,7 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, rplm, vplm, rplt, nenc, index1, index2, lvdotr) use base, only: base_parameters implicit none - class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s + class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s integer(I4B), intent(in) :: nplm !! Total number of fully interacting massive bodies integer(I4B), intent(in) :: nplt !! Total number of partially interacting masive bodies (GM < GMTINY) real(DP), dimension(:,:), intent(in) :: rplm !! Position vectors of fully interacting massive bodies @@ -141,7 +145,7 @@ end subroutine encounter_check_all_plplm module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) use base, only: base_parameters implicit none - class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s + class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s integer(I4B), intent(in) :: npl !! Total number of massive bodies integer(I4B), intent(in) :: ntp !! Total number of test particles real(DP), dimension(:,:), intent(in) :: rpl !! Position vectors of massive bodies @@ -216,14 +220,14 @@ end subroutine encounter_check_sweep_aabb_single_list module subroutine encounter_io_netcdf_dump(self, param) implicit none - class(encounter_storage(*)), intent(inout) :: self !! Encounter storage object + class(encounter_storage), intent(inout) :: self !! Encounter storage object class(base_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine encounter_io_netcdf_dump module subroutine encounter_io_netcdf_initialize_output(self, param) implicit none class(encounter_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset - class(base_parameters), intent(in) :: param + class(base_parameters), intent(in) :: param end subroutine encounter_io_netcdf_initialize_output module subroutine encounter_io_netcdf_open(self, param, readonly) @@ -235,9 +239,9 @@ end subroutine encounter_io_netcdf_open module subroutine encounter_io_netcdf_write_frame_snapshot(self, history, param) implicit none - class(encounter_snapshot), intent(in) :: self !! Swiftest encounter structure - class(encounter_storage(*)), intent(inout) :: history !! Encounter storage object - class(base_parameters), intent(inout) :: param !! Current run configuration parameters + class(encounter_snapshot), intent(in) :: self !! Swiftest encounter structure + class(encounter_storage), intent(inout) :: history !! Encounter storage object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine encounter_io_netcdf_write_frame_snapshot module subroutine encounter_util_setup_aabb(self, n, n_last) @@ -271,11 +275,26 @@ module subroutine encounter_util_dealloc_aabb(self) class(encounter_bounding_box_1D), intent(inout) :: self !!Bounding box structure along a single dimension end subroutine encounter_util_dealloc_aabb + module subroutine encounter_util_dealloc_bounding_box(self) + implicit none + class(encounter_bounding_box), intent(inout) :: self !! Bounding box structure + end subroutine encounter_util_dealloc_bounding_box + module subroutine encounter_util_dealloc_list(self) implicit none class(encounter_list), intent(inout) :: self !! Swiftest encounter list object end subroutine encounter_util_dealloc_list + module subroutine encounter_util_dealloc_snapshot(self) + implicit none + class(encounter_snapshot), intent(inout) :: self !! Encounter shapshot object + end subroutine encounter_util_dealloc_snapshot + + module subroutine encounter_util_dealloc_storage(self) + implicit none + class(encounter_storage), intent(inout) :: self !! Swiftest storage object + end subroutine encounter_util_dealloc_storage + module subroutine encounter_util_get_idvalues_snapshot(self, idvals) implicit none class(encounter_snapshot), intent(in) :: self !! Encounter snapshot object @@ -283,14 +302,14 @@ module subroutine encounter_util_get_idvalues_snapshot(self, idvals) end subroutine encounter_util_get_idvalues_snapshot module subroutine encounter_util_get_vals_storage(self, idvals, tvals) - class(encounter_storage(*)), intent(in) :: self !! Encounter storages object + class(encounter_storage), intent(in) :: self !! Encounter storages object integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values in all snapshots real(DP), dimension(:), allocatable, intent(out) :: tvals !! Array of all time values in all snapshots end subroutine encounter_util_get_vals_storage module subroutine encounter_util_index_map(self) implicit none - class(encounter_storage(*)), intent(inout) :: self !! Encounter storage object + class(encounter_storage), intent(inout) :: self !! Encounter storage object end subroutine encounter_util_index_map module subroutine encounter_util_resize_list(self, nnew) @@ -301,11 +320,11 @@ end subroutine encounter_util_resize_list module subroutine encounter_util_snapshot(self, param, nbody_system, t, arg) implicit none - class(encounter_storage(*)), intent(inout) :: self !! Swiftest storage object - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object to store - real(DP), intent(in), optional :: t !! Time of snapshot if different from system time - character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) + class(encounter_storage), intent(inout) :: self !! Swiftest storage object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters + class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object to store + real(DP), intent(in), optional :: t !! Time of snapshot if different from system time + character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) end subroutine encounter_util_snapshot module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestructive) @@ -333,6 +352,19 @@ subroutine encounter_final_aabb(self) return end subroutine encounter_final_aabb + subroutine encounter_final_bounding_box(self) + !! author: David A. Minton + !! + !! Finalize the bounding box object + implicit none + ! Arguments + type(encounter_bounding_box), intent(inout) :: self + + call self%dealloc() + + return + end subroutine encounter_final_bounding_box + subroutine encounter_final_netcdf_parameters(self) !! author: David A. Minton @@ -370,15 +402,9 @@ subroutine encounter_final_storage(self) !! Deallocates allocatable arrays in an encounter snapshot implicit none ! Arguments - type(encounter_storage(*)), intent(inout) :: self !! Encounter storage object - ! Internals - integer(I4B) :: i - - do i = 1, self%nframes - if (allocated(self%frame(i)%item)) deallocate(self%frame(i)%item) - end do + type(encounter_storage), intent(inout) :: self !! Encounter storage object - return + call self%dealloc() return end subroutine encounter_final_storage diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index 0377c290c..94ea7d942 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -85,6 +85,7 @@ module subroutine encounter_util_dealloc_aabb(self) ! Arguments class(encounter_bounding_box_1D), intent(inout) :: self + self%n = 0 if (allocated(self%ind)) deallocate(self%ind) if (allocated(self%ibeg)) deallocate(self%ibeg) if (allocated(self%iend)) deallocate(self%iend) @@ -93,6 +94,24 @@ module subroutine encounter_util_dealloc_aabb(self) end subroutine encounter_util_dealloc_aabb + module subroutine encounter_util_dealloc_bounding_box(self) + !! author: David A. Minton + !! + !! Deallocates all allocatables + implicit none + ! Arguments + class(encounter_bounding_box), intent(inout) :: self !! Bounding box structure + ! Internals + integer(I4B) :: i + + do i = 1,NDIM + call self%aabb(i)%dealloc() + end do + + return + end subroutine encounter_util_dealloc_bounding_box + + module subroutine encounter_util_dealloc_list(self) !! author: David A. Minton !! @@ -101,6 +120,8 @@ module subroutine encounter_util_dealloc_list(self) ! Arguments class(encounter_list), intent(inout) :: self + self%nenc = 0 + if (allocated(self%tcollision)) deallocate(self%tcollision) if (allocated(self%lclosest)) deallocate(self%lclosest) if (allocated(self%lvdotr)) deallocate(self%lvdotr) @@ -119,6 +140,38 @@ module subroutine encounter_util_dealloc_list(self) end subroutine encounter_util_dealloc_list + module subroutine encounter_util_dealloc_snapshot(self) + !! author: David A. Minton + !! + !! Deallocates all allocatables + implicit none + ! Arguments + class(encounter_snapshot), intent(inout) :: self !! Encounter shapshot object + + if (allocated(self%pl)) deallocate(self%pl) + if (allocated(self%tp)) deallocate(self%tp) + + return + end subroutine encounter_util_dealloc_snapshot + + + module subroutine encounter_util_dealloc_storage(self) + !! author: David A. Minton + !! + !! Resets a storage object by deallocating all items and resetting the frame counter to 0 + use base, only : base_util_dealloc_storage + implicit none + ! Arguments + class(encounter_storage), intent(inout) :: self !! Swiftest storage object + + if (allocated(self%nc)) deallocate(self%nc) + + call base_util_dealloc_storage(self) + + return + end subroutine encounter_util_dealloc_storage + + module subroutine encounter_util_get_idvalues_snapshot(self, idvals) !! author: David A. Minton !! @@ -164,7 +217,7 @@ module subroutine encounter_util_get_vals_storage(self, idvals, tvals) !! !! Gets the id values in a self object, regardless of whether it is encounter of collision ! Argument - class(encounter_storage(*)), intent(in) :: self !! Encounter storages object + class(encounter_storage), intent(in) :: self !! Encounter storages object integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values in all snapshots real(DP), dimension(:), allocatable, intent(out) :: tvals !! Array of all time values in all snapshots ! Internals @@ -224,7 +277,7 @@ module subroutine encounter_util_index_map(self) !! Basically this will make a unique list of ids that exist in all of the saved snapshots implicit none ! Arguments - class(encounter_storage(*)), intent(inout) :: self !! Swiftest storage object + class(encounter_storage), intent(inout) :: self !! Swiftest storage object ! Internals integer(I4B), dimension(:), allocatable :: idvals real(DP), dimension(:), allocatable :: tvals @@ -412,65 +465,19 @@ module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestru end subroutine encounter_util_spill_list - subroutine encounter_util_save_snapshot(encounter_history, snapshot) - !! author: David A. Minton - !! - !! Checks the current size of the encounter storage against the required size and extends it by a factor of 2 more than requested if it is too small. - !! Note: The reason to extend it by a factor of 2 is for performance. When there are many enounters per step, resizing every time you want to add an - !! encounter takes significant computational effort. Resizing by a factor of 2 is a tradeoff between performance (fewer resize calls) and memory managment - !! Memory usage grows by a factor of 2 each time it fills up, but no more. - implicit none - ! Arguments - class(encounter_storage(*)), allocatable, intent(inout) :: encounter_history !! SyMBA encounter storage object - class(encounter_snapshot), intent(in) :: snapshot !! Encounter snapshot object - ! Internals - type(encounter_storage(nframes=:)), allocatable :: tmp - integer(I4B) :: i, nnew, nold, nbig - - ! Advance the snapshot frame counter - encounter_history%iframe = encounter_history%iframe + 1 - - ! Check to make sure the current encounter_history object is big enough. If not, grow it by a factor of 2 - nnew = encounter_history%iframe - nold = encounter_history%nframes - - if (nnew > nold) then - nbig = nold - do while (nbig < nnew) - nbig = nbig * 2 - end do - allocate(encounter_storage(nbig) :: tmp) - tmp%iframe = encounter_history%iframe - call move_alloc(encounter_history%nc, tmp%nc) - - do i = 1, nold - if (allocated(encounter_history%frame(i)%item)) call move_alloc(encounter_history%frame(i)%item, tmp%frame(i)%item) - end do - deallocate(encounter_history) - call move_alloc(tmp,encounter_history) - end if - - ! Find out which time slot this belongs in by searching for an existing slot - ! with the same value of time or the first available one - encounter_history%frame(nnew) = snapshot - - return - end subroutine encounter_util_save_snapshot - - module subroutine encounter_util_snapshot(self, param, nbody_system, t, arg) !! author: David A. Minton !! !! Takes a minimal snapshot of the state of the system during an encounter so that the trajectories !! can be played back through the encounter - use symba + use symba, only : symba_pl, symba_tp, symba_nbody_system implicit none ! Internals - class(encounter_storage(*)), intent(inout) :: self !! Swiftest storage object - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object to store - real(DP), intent(in), optional :: t !! Time of snapshot if different from system time - character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) + class(encounter_storage), intent(inout) :: self !! Swiftest storage object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters + class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object to store + real(DP), intent(in), optional :: t !! Time of snapshot if different from system time + character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) ! Arguments class(encounter_snapshot), allocatable :: snapshot integer(I4B) :: i, pi, pj, k, npl_snap, ntp_snap, iflag @@ -586,11 +593,8 @@ module subroutine encounter_util_snapshot(self, param, nbody_system, t, arg) end if ! Save the snapshot - select type (encounter_history => nbody_system%encounter_history) - class is (encounter_storage(*)) - encounter_history%nid = encounter_history%nid + ntp_snap + npl_snap - call encounter_util_save_snapshot(nbody_system%encounter_history,snapshot) - end select + self%nid = self%nid + ntp_snap + npl_snap + call self%save(snapshot) case("closest") associate(plpl_encounter => nbody_system%plpl_encounter, pltp_encounter => nbody_system%pltp_encounter) if (any(plpl_encounter%lclosest(:))) then @@ -655,7 +659,7 @@ module subroutine encounter_util_snapshot(self, param, nbody_system, t, arg) pl_snap%vh(:,2) = vb(:,2) + vcom(:) call pl_snap%sort("id", ascending=.true.) - call encounter_util_save_snapshot(nbody_system%encounter_history,snapshot) + call self%save(snapshot) end if end do @@ -683,5 +687,4 @@ module subroutine encounter_util_snapshot(self, param, nbody_system, t, arg) end subroutine encounter_util_snapshot - end submodule s_encounter_util diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 9e712c81e..227b295d7 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -46,7 +46,7 @@ module subroutine fraggle_generate(self, nbody_system, param, t) message = "Supercatastrophic disruption between" case default write(*,*) "Error in swiftest_collision, unrecognized collision regime" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select call self%set_mass_dist(param) call self%disrupt(nbody_system, param, t) @@ -94,7 +94,6 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur logical :: lk_plpl, lfailure_local logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes real(DP) :: dE, dL - integer(I4B) :: i character(len=STRMAX) :: message real(DP), parameter :: fail_scale_initial = 1.001_DP @@ -108,9 +107,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur class is (swiftest_nbody_system) select type(param) class is (swiftest_parameters) - select type(fragments => self%fragments) - class is (fraggle_fragments(*)) - associate(impactors => self%impactors, nfrag => fragments%nbody, pl => nbody_system%pl) + associate(impactors => self%impactors, fragments => self%fragments, nfrag => self%fragments%nbody, pl => nbody_system%pl) write(message,*) nfrag call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle generating " // trim(adjustl(message)) // " fragments.") @@ -124,7 +121,6 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur call ieee_set_flag(ieee_all, .false.) ! Set all fpe flags to quiet call self%set_natural_scale() - call fragments%reset() lfailure_local = .false. call self%get_energy_and_momentum(nbody_system, param, phase="before") self%fail_scale = fail_scale_initial @@ -149,8 +145,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur end associate end select end select - end select - call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily + call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) ! Restore the original halting modes return end subroutine fraggle_generate_disrupt diff --git a/src/fraggle/fraggle_module.f90 b/src/fraggle/fraggle_module.f90 index b0eb38a76..84c88ba6b 100644 --- a/src/fraggle/fraggle_module.f90 +++ b/src/fraggle/fraggle_module.f90 @@ -15,15 +15,6 @@ module fraggle implicit none public - !> Class definition for the variables that describe a collection of fragments by Fraggle barycentric coordinates - type, extends(collision_fragments) :: fraggle_fragments - contains - - procedure :: reset => fraggle_util_reset_fragments !! Resets all position and velocity-dependent fragment quantities in order to do a fresh calculation (does not reset mass, radius, or other values that get set prior to the call to fraggle_generate) - final :: fraggle_final_fragments !! Finalizer will deallocate all allocatables - end type fraggle_fragments - - 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 @@ -32,8 +23,7 @@ module fraggle procedure :: hitandrun => fraggle_generate_hitandrun !! Generates either a pure hit and run, or one in which the runner is disrupted procedure :: set_mass_dist => fraggle_util_set_mass_dist !! Sets the distribution of mass among the fragments depending on the regime type procedure :: setup_fragments => fraggle_util_setup_fragments_system !! Initializer for the fragments of the collision system. - procedure :: reset => fraggle_util_reset_system !! Deallocates all allocatables - final :: fraggle_final_system !! Finalizer will deallocate all allocatables + !procedure :: reset => fraggle_util_dealloc_system !! Deallocates all allocatables end type collision_fraggle interface @@ -88,64 +78,11 @@ module subroutine fraggle_util_setup_fragments_system(self, nfrag) integer(I4B), intent(in) :: nfrag !! Number of fragments to create end subroutine fraggle_util_setup_fragments_system - module subroutine fraggle_util_reset_fragments(self) - implicit none - class(fraggle_fragments(*)), intent(inout) :: self - end subroutine fraggle_util_reset_fragments - - module subroutine fraggle_util_reset_system(self) - implicit none - class(collision_fraggle), intent(inout) :: self !! Collision system object - end subroutine fraggle_util_reset_system - module subroutine fraggle_util_set_mass_dist(self, param) implicit none class(collision_fraggle), intent(inout) :: self !! Fraggle collision object class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters end subroutine fraggle_util_set_mass_dist end interface - - contains - - subroutine fraggle_final_fragments(self) - !! author: David A. Minton - !! - !! Finalizer will deallocate all allocatables - implicit none - ! Arguments - type(fraggle_fragments(*)), intent(inout) :: self !! Fraggle encountar storage object - - if (allocated(self%info)) deallocate(self%info) - - return - end subroutine fraggle_final_fragments - - - subroutine fraggle_final_impactors(self) - !! author: David A. Minton - !! - !! Finalizer will deallocate all allocatables - implicit none - ! Arguments - type(collision_impactors), intent(inout) :: self !! Fraggle impactors object - call self%reset() - return - end subroutine fraggle_final_impactors - - - subroutine fraggle_final_system(self) - !! author: David A. Minton - !! - !! Finalizer will deallocate all allocatables - implicit none - ! Arguments - type(collision_fraggle), intent(inout) :: self !! Collision impactors storage object - - call self%reset() - if (allocated(self%impactors)) deallocate(self%impactors) - if (allocated(self%fragments)) deallocate(self%fragments) - - return - end subroutine fraggle_final_system - + end module fraggle \ No newline at end of file diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 667ab6a15..e9e2dc9bf 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -11,53 +11,6 @@ use swiftest contains - module subroutine fraggle_util_reset_fragments(self) - !! author: David A. Minton - !! - !! Resets all position and velocity-dependent fragment quantities in order to do a fresh calculation (does not reset mass, radius, or other values that get set prior to the call to fraggle_generate) - implicit none - ! Arguments - class(fraggle_fragments(*)), intent(inout) :: self - - self%rc(:,:) = 0.0_DP - self%vc(:,:) = 0.0_DP - self%rh(:,:) = 0.0_DP - self%vh(:,:) = 0.0_DP - self%rb(:,:) = 0.0_DP - self%vb(:,:) = 0.0_DP - self%rot(:,:) = 0.0_DP - self%r_unit(:,:) = 0.0_DP - self%t_unit(:,:) = 0.0_DP - self%n_unit(:,:) = 0.0_DP - - self%rmag(:) = 0.0_DP - self%rotmag(:) = 0.0_DP - - return - end subroutine fraggle_util_reset_fragments - - - module subroutine fraggle_util_reset_system(self) - !! author: David A. Minton - !! - !! Resets the collider system and deallocates all allocatables - implicit none - ! Arguments - class(collision_fraggle), intent(inout) :: self !! Collision system object - - self%dscale = 1.0_DP - self%mscale = 1.0_DP - self%tscale = 1.0_DP - self%vscale = 1.0_DP - self%Escale = 1.0_DP - self%Lscale = 1.0_DP - - call collision_util_reset_system(self) - - return - end subroutine fraggle_util_reset_system - - module subroutine fraggle_util_set_mass_dist(self, param) !! author: David A. Minton !! @@ -146,8 +99,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) write(*,*) "collision_util_set_mass_dist_fragments error: Unrecognized regime code",impactors%regime end select - select type(fragments => self%fragments) - class is (collision_fragments(*)) + associate(fragments => self%fragments) fragments%mtot = mtot allocate(mass, mold=fragments%mass) @@ -223,7 +175,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) fragments%origin_body(3:nfrag) = 2 end if - end select + end associate end associate @@ -242,7 +194,7 @@ module subroutine fraggle_util_setup_fragments_system(self, nfrag) integer(I4B), intent(in) :: nfrag !! Number of fragments to create if (allocated(self%fragments)) deallocate(self%fragments) - allocate(fraggle_fragments(nbody=nfrag) :: self%fragments) + allocate(collision_fragments(nbody=nfrag) :: self%fragments) self%fragments%nbody = nfrag return diff --git a/src/globals/globals_module.f90 b/src/globals/globals_module.f90 index 0f03b4200..fa6b18582 100644 --- a/src/globals/globals_module.f90 +++ b/src/globals/globals_module.f90 @@ -41,7 +41,7 @@ module globals integer(I4B), parameter :: LOWERCASE_END = iachar('z') !! ASCII character set parameter for lower to upper conversion - end of lowercase integer(I4B), parameter :: UPPERCASE_OFFSET = iachar('A') - iachar('a') !! ASCII character set parameter for lower to upper conversion - offset between upper and lower - real(SP), parameter :: VERSION_NUMBER = 1.0_SP !! swiftest version + real(SP), parameter :: VERSION_NUMBER = 1.0_SP !! Swiftest version !> Symbolic name for integrator types character(*), parameter :: UNKNOWN_INTEGRATOR = "UKNOWN INTEGRATOR" diff --git a/src/helio/helio_module.f90 b/src/helio/helio_module.f90 index a9d41b1a9..cd5581837 100644 --- a/src/helio/helio_module.f90 +++ b/src/helio/helio_module.f90 @@ -22,7 +22,6 @@ module helio contains procedure :: step => helio_step_system !! Advance the Helio nbody system forward in time by one step procedure :: initialize => helio_util_setup_initialize_system !! Performs Helio-specific initilization steps, including converting to DH coordinates - final :: helio_final_system !! Finalizes the Helio nbody_system object - deallocates all allocatables end type helio_nbody_system @@ -44,7 +43,6 @@ module helio procedure :: accel => helio_kick_getacch_pl !! Compute heliocentric accelerations of massive bodies procedure :: kick => helio_kick_vb_pl !! Kicks the barycentric velocities procedure :: step => helio_step_pl !! Steps the body forward one stepsize - final :: helio_final_pl !! Finalizes the Helio massive body object - deallocates all allocatables end type helio_pl @@ -58,7 +56,6 @@ module helio procedure :: accel => helio_kick_getacch_tp !! Compute heliocentric accelerations of massive bodies procedure :: kick => helio_kick_vb_tp !! Kicks the barycentric velocities procedure :: step => helio_step_tp !! Steps the body forward one stepsize - final :: helio_final_tp !! Finalizes the Helio test particle object - deallocates all allocatables end type helio_tp interface @@ -201,47 +198,4 @@ module subroutine helio_step_tp(self, nbody_system, param, t, dt) end subroutine helio_step_tp end interface - contains - - subroutine helio_final_pl(self) - !! author: David A. Minton - !! - !! Finalize the Helio massive body object - deallocates all allocatables - implicit none - ! Arguments - type(helio_pl), intent(inout) :: self !! Helio massive body object - - call self%dealloc() - - return - end subroutine helio_final_pl - - - subroutine helio_final_system(self) - !! author: David A. Minton - !! - !! Finalize the Helio nbody system object - deallocates all allocatables - implicit none - ! Arguments - type(helio_nbody_system), intent(inout) :: self !! Helio nbody system object - - call whm_final_system(self%whm_nbody_system) - - return - end subroutine helio_final_system - - - subroutine helio_final_tp(self) - !! author: David A. Minton - !! - !! Finalize the Helio test particle object - deallocates all allocatables - implicit none - ! Arguments - type(helio_tp), intent(inout) :: self !! Helio test particle object - - call self%dealloc() - - return - end subroutine helio_final_tp - end module helio diff --git a/src/misc/solver_module.f90 b/src/misc/solver_module.f90 index 4fbf7b7b3..d8fc85c89 100644 --- a/src/misc/solver_module.f90 +++ b/src/misc/solver_module.f90 @@ -244,7 +244,7 @@ module function solve_rkf45(f, y0in, t1, dt0, tol) result(y1) if (s >= 1.0_DP) exit ! Good step! if (i == MAXREDUX) then write(*,*) "Something has gone wrong in util_solve_rkf45!! Step size reduction has gone too far this time!" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end if end do diff --git a/src/netcdf_io/netcdf_io_implementations.f90 b/src/netcdf_io/netcdf_io_implementations.f90 index f917058f1..e279ef840 100644 --- a/src/netcdf_io/netcdf_io_implementations.f90 +++ b/src/netcdf_io/netcdf_io_implementations.f90 @@ -23,7 +23,7 @@ module subroutine netcdf_io_check(status, call_identifier) if(status /= NF90_NOERR) then if (present(call_identifier)) write(*,*) "NetCDF error in ",trim(call_identifier) write(*,*) trim(nf90_strerror(status)) - call util_exit(FAILURE) + call base_util_exit(FAILURE) end if return diff --git a/src/operator/operator_cross.f90 b/src/operator/operator_cross.f90 index cec60d23b..f784a1c98 100644 --- a/src/operator/operator_cross.f90 +++ b/src/operator/operator_cross.f90 @@ -8,6 +8,7 @@ !! If not, see: https://www.gnu.org/licenses. submodule(operators) s_operator_cross + use, intrinsic :: ieee_exceptions use swiftest !! author: David A. Minton !! @@ -21,6 +22,8 @@ pure module function operator_cross_sp(A, B) result(C) implicit none real(SP), dimension(:), intent(in) :: A, B real(SP), dimension(3) :: C + + call ieee_set_halting_mode(ieee_underflow, .false.) C(1) = A(2) * B(3) - A(3) * B(2) C(2) = A(3) * B(1) - A(1) * B(3) C(3) = A(1) * B(2) - A(2) * B(1) @@ -31,6 +34,8 @@ pure module function operator_cross_dp(A, B) result(C) implicit none real(DP), dimension(:), intent(in) :: A, B real(DP), dimension(3) :: C + + call ieee_set_halting_mode(ieee_underflow, .false.) C(1) = A(2) * B(3) - A(3) * B(2) C(2) = A(3) * B(1) - A(1) * B(3) C(3) = A(1) * B(2) - A(2) * B(1) @@ -41,6 +46,8 @@ pure module function operator_cross_qp(A, B) result(C) implicit none real(QP), dimension(:), intent(in) :: A, B real(QP), dimension(3) :: C + + call ieee_set_halting_mode(ieee_underflow, .false.) C(1) = A(2) * B(3) - A(3) * B(2) C(2) = A(3) * B(1) - A(1) * B(3) C(3) = A(1) * B(2) - A(2) * B(1) @@ -51,6 +58,7 @@ pure module function operator_cross_i1b(A, B) result(C) implicit none integer(I1B), dimension(:), intent(in) :: A, B integer(I1B), dimension(3) :: C + C(1) = A(2) * B(3) - A(3) * B(2) C(2) = A(3) * B(1) - A(1) * B(3) C(3) = A(1) * B(2) - A(2) * B(1) @@ -61,6 +69,7 @@ pure module function operator_cross_i2b(A, B) result(C) implicit none integer(I2B), dimension(:), intent(in) :: A, B integer(I2B), dimension(3) :: C + C(1) = A(2) * B(3) - A(3) * B(2) C(2) = A(3) * B(1) - A(1) * B(3) C(3) = A(1) * B(2) - A(2) * B(1) diff --git a/src/operator/operator_mag.f90 b/src/operator/operator_mag.f90 index cdbd2b773..8bf6bff5b 100644 --- a/src/operator/operator_mag.f90 +++ b/src/operator/operator_mag.f90 @@ -8,18 +8,20 @@ !! If not, see: https://www.gnu.org/licenses. submodule(operators) s_operator_mag + use, intrinsic :: ieee_exceptions !! author: David A. Minton !! !! Contains implementations for the .mag. operator for all defined real types !! Computes the magnitude of a vector or array of vectors using norm2 !! Single vector implementations: B = .mag. A(1:3) !! Vector list implementations: B(:) = .mag. A(1:3, :) - contains +contains pure module function operator_mag_sp(A) result(B) implicit none real(SP), dimension(:), intent(in) :: A real(SP) :: B + call ieee_set_halting_mode(ieee_underflow, .false.) B = norm2(A(:)) return end function operator_mag_sp @@ -28,6 +30,7 @@ pure module function operator_mag_dp(A) result(B) implicit none real(DP), dimension(:), intent(in) :: A real(DP) :: B + call ieee_set_halting_mode(ieee_underflow, .false.) B = norm2(A(:)) return end function operator_mag_dp @@ -40,6 +43,7 @@ pure module function operator_mag_el_sp(A) result(B) n = size(A, 2) if (allocated(B)) deallocate(B) allocate(B(n)) + call ieee_set_halting_mode(ieee_underflow, .false.) do concurrent (i=1:n) B(i) = norm2(A(:, i)) end do @@ -54,6 +58,7 @@ pure module function operator_mag_el_dp(A) result(B) n = size(A, 2) if (allocated(B)) deallocate(B) allocate(B(n)) + call ieee_set_halting_mode(ieee_underflow, .false.) do concurrent (i=1:n) B(i) = norm2(A(:, i)) end do @@ -68,6 +73,7 @@ pure module function operator_mag_el_qp(A) result(B) n = size(A, 2) if (allocated(B)) deallocate(B) allocate(B(n)) + call ieee_set_halting_mode(ieee_underflow, .false.) do concurrent (i=1:n) B(i) = norm2(A(:, i)) end do diff --git a/src/operator/operator_unit.f90 b/src/operator/operator_unit.f90 index 2b75e3851..39f3fce53 100644 --- a/src/operator/operator_unit.f90 +++ b/src/operator/operator_unit.f90 @@ -8,13 +8,14 @@ !! If not, see: https://www.gnu.org/licenses. submodule(operators) s_operator_unit + use, intrinsic :: ieee_exceptions !! author: David A. Minton !! !! Contains implementations for the .unit. operator for all defined real types !! Returns a unit vector or array of unit vectors from an input vector or array of vectors !! Single vector implementations: B = .unit. A(1:NDIM) !! Vector list implementations: B(:) = .unit. A(1:NDIM, :) - contains +contains pure module function operator_unit_sp(A) result(B) implicit none @@ -24,6 +25,7 @@ pure module function operator_unit_sp(A) result(B) ! Internals real(SP) :: Amag + call ieee_set_halting_mode(ieee_underflow, .false.) Amag = norm2(A(:)) if (Amag > tiny(1._SP)) then B(:) = A(:) / Amag @@ -43,6 +45,7 @@ pure module function operator_unit_dp(A) result(B) ! Internals real(DP) :: Amag + call ieee_set_halting_mode(ieee_underflow, .false.) Amag = norm2(A(:)) if (Amag > tiny(1._DP)) then B(:) = A(:) / Amag @@ -62,6 +65,7 @@ pure module function operator_unit_qp(A) result(B) ! Internals real(QP) :: Amag + call ieee_set_halting_mode(ieee_underflow, .false.) Amag = norm2(A(:)) if (Amag > tiny(1._QP)) then B(:) = A(:) / Amag diff --git a/src/rmvs/rmvs_module.f90 b/src/rmvs/rmvs_module.f90 index 3417ea0e2..733a81e60 100644 --- a/src/rmvs/rmvs_module.f90 +++ b/src/rmvs/rmvs_module.f90 @@ -32,9 +32,9 @@ module rmvs real(DP), dimension(:,:), allocatable :: vbeg !! Planet velocities at beginning ot step contains !> Replace the abstract procedures with concrete ones + procedure :: dealloc => rmvs_util_dealloc_system !! Performs RMVS-specific deallocation procedure :: initialize => rmvs_util_setup_initialize_system !! Performs RMVS-specific initilization steps, including generating the close encounter planetocentric structures - procedure :: step => rmvs_step_system !! Advance the RMVS nbody system forward in time by one step - final :: rmvs_final_system !! Finalizes the RMVS nbody system object - deallocates all allocatables + procedure :: step => rmvs_step_system !! Advance the RMVS nbody system forward in time by one step end type rmvs_nbody_system type, private :: rmvs_interp @@ -187,6 +187,11 @@ module subroutine rmvs_util_dealloc_pl(self) class(rmvs_pl), intent(inout) :: self !! RMVS massive body object end subroutine rmvs_util_dealloc_pl + module subroutine rmvs_util_dealloc_system(self) + implicit none + class(rmvs_nbody_system), intent(inout) :: self + end subroutine rmvs_util_dealloc_system + module subroutine rmvs_util_dealloc_tp(self) implicit none class(rmvs_tp), intent(inout) :: self !! RMVS test particle object @@ -314,21 +319,6 @@ subroutine rmvs_final_pl(self) end subroutine rmvs_final_pl - subroutine rmvs_final_system(self) - !! author: David A. Minton - !! - !! Finalize the RMVS nbody system object - deallocates all allocatables - implicit none - ! Arguments - type(rmvs_nbody_system), intent(inout) :: self !! RMVS nbody system object - - if (allocated(self%vbeg)) deallocate(self%vbeg) - call whm_final_system(self%whm_nbody_system) - - return - end subroutine rmvs_final_system - - subroutine rmvs_final_tp(self) !! author: David A. Minton !! diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 19b8779f6..bcb38797a 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -119,7 +119,7 @@ subroutine rmvs_interp_out(cb, pl, dt) vtmp(:,i),new_line('a'), & " STOPPING " call swiftest_io_log_one_message(COLLISION_LOG_OUT,message) - call util_exit(FAILURE) + call base_util_exit(FAILURE) end if end do end if @@ -142,7 +142,7 @@ subroutine rmvs_interp_out(cb, pl, dt) vtmp(:,i), new_line('a'), & " STOPPING " call swiftest_io_log_one_message(COLLISION_LOG_OUT,message) - call util_exit(FAILURE) + call base_util_exit(FAILURE) end if end do end if @@ -289,7 +289,7 @@ subroutine rmvs_interp_in(cb, pl, nbody_system, param, dt, outer_index) vtmp(:,i), new_line('a'), & " STOPPING " call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) - call util_exit(failure) + call base_util_exit(failure) end if end do end if @@ -313,7 +313,7 @@ subroutine rmvs_interp_in(cb, pl, nbody_system, param, dt, outer_index) write(*, *) xtmp(:,i) write(*, *) vtmp(:,i) write(*, *) " STOPPING " - call util_exit(failure) + call base_util_exit(failure) end if end do end if diff --git a/src/rmvs/rmvs_util.f90 b/src/rmvs/rmvs_util.f90 index 07e4f9a51..b9462840b 100644 --- a/src/rmvs/rmvs_util.f90 +++ b/src/rmvs/rmvs_util.f90 @@ -39,7 +39,7 @@ module subroutine rmvs_util_append_pl(self, source, lsource_mask) end associate class default write(*,*) "Invalid object passed to the append method. Source must be of class rmvs_pl or its descendents!" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select return @@ -68,7 +68,7 @@ module subroutine rmvs_util_append_tp(self, source, lsource_mask) end associate class default write(*,*) "Invalid object passed to the append method. Source must be of class rmvs_tp or its descendents!" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select return @@ -126,6 +126,22 @@ module subroutine rmvs_util_dealloc_pl(self) end subroutine rmvs_util_dealloc_pl + module subroutine rmvs_util_dealloc_system(self) + !! author: David A. Minton + !! + !! Deallocates all allocatables and resets all values to defaults. Acts as a base for a finalizer + implicit none + ! Arguments + class(rmvs_nbody_system), intent(inout) :: self + + self%lplanetocentric = .false. + if (allocated(self%vbeg)) deallocate(self%vbeg) + call self%whm_nbody_system%dealloc() + + return + end subroutine rmvs_util_dealloc_system + + module subroutine rmvs_util_dealloc_tp(self) !! author: David A. Minton !! @@ -174,7 +190,7 @@ module subroutine rmvs_util_fill_pl(self, inserts, lfill_list) call whm_util_fill_pl(keeps, inserts, lfill_list) class default write(*,*) "Invalid object passed to the fill method. Source must be of class rmvs_pl or its descendents!" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select end associate @@ -204,7 +220,7 @@ module subroutine rmvs_util_fill_tp(self, inserts, lfill_list) call swiftest_util_fill_tp(keeps, inserts, lfill_list) ! Note: whm_tp does not have its own fill method, so we skip back to the base class class default write(*,*) "Invalid object passed to the fill method. Source must be of class rmvs_tp or its descendents!" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select end associate @@ -564,7 +580,7 @@ module subroutine rmvs_util_spill_pl(self, discards, lspill_list, ldestructive) call whm_util_spill_pl(keeps, discards, lspill_list, ldestructive) class default write(*,*) "Invalid object passed to the spill method. Source must be of class rmvs_pl or its descendents!" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select end associate @@ -595,7 +611,7 @@ module subroutine rmvs_util_spill_tp(self, discards, lspill_list, ldestructive) call swiftest_util_spill_tp(keeps, discards, lspill_list, ldestructive) class default write(*,*) "Invalid object passed to the spill method. Source must be of class rmvs_tp or its descendents!" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select end associate diff --git a/src/swiftest/swiftest_driver.f90 b/src/swiftest/swiftest_driver.f90 index d75b8a561..8c6a63f6d 100644 --- a/src/swiftest/swiftest_driver.f90 +++ b/src/swiftest/swiftest_driver.f90 @@ -125,5 +125,5 @@ program swiftest_driver end associate end associate - call util_exit(SUCCESS) + call base_util_exit(SUCCESS) end program swiftest_driver diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index da3ccb2d2..66d22816a 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -178,7 +178,7 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal) ! Save the frame of data to the bin file in the slot just after the present one for diagnostics call self%write_frame(nc, param) call nc%close() - call util_exit(FAILURE) + call base_util_exit(FAILURE) end if end if end associate @@ -187,7 +187,7 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal) 667 continue write(*,*) "Error writing energy and momentum tracking file: " // trim(adjustl(errmsg)) - call util_exit(FAILURE) + call base_util_exit(FAILURE) end subroutine swiftest_io_conservation_report @@ -286,7 +286,7 @@ module subroutine swiftest_io_dump_param(self, param_file_name) 667 continue write(*,*) "Error opening parameter dump file " // trim(adjustl(errmsg)) - call util_exit(FAILURE) + call base_util_exit(FAILURE) end subroutine swiftest_io_dump_param @@ -339,7 +339,7 @@ module subroutine swiftest_io_dump_storage(self, param) !! cadence is not divisible by the total number of loops). implicit none ! Arguments - class(swiftest_storage(*)), intent(inout) :: self !! Swiftest simulation history storage object + class(swiftest_storage), intent(inout) :: self !! Swiftest simulation history storage object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i @@ -386,18 +386,18 @@ module subroutine swiftest_io_get_args(integrator, param_file_name, display_styl do i = 1,narg call get_command_argument(i, arg(i), status = ierr(i)) end do - if (any(ierr /= 0)) call util_exit(USAGE) + if (any(ierr /= 0)) call base_util_exit(USAGE) else - call util_exit(USAGE) + call base_util_exit(USAGE) end if if (narg == 1) then if (arg(1) == '-v' .or. arg(1) == '--version') then call swiftest_util_version() else if (arg(1) == '-h' .or. arg(1) == '--help') then - call util_exit(HELP) + call base_util_exit(HELP) else - call util_exit(USAGE) + call base_util_exit(USAGE) end if else if (narg >= 2) then call swiftest_io_toupper(arg(1)) @@ -421,7 +421,7 @@ module subroutine swiftest_io_get_args(integrator, param_file_name, display_styl case default integrator = UNKNOWN_INTEGRATOR write(*,*) trim(adjustl(arg(1))) // ' is not a valid integrator.' - call util_exit(USAGE) + call base_util_exit(USAGE) end select param_file_name = trim(adjustl(arg(2))) end if @@ -432,7 +432,7 @@ module subroutine swiftest_io_get_args(integrator, param_file_name, display_styl call swiftest_io_toupper(arg(3)) display_style = trim(adjustl(arg(3))) else - call util_exit(USAGE) + call base_util_exit(USAGE) end if return @@ -804,7 +804,7 @@ module subroutine swiftest_io_netcdf_initialize_output(self, param) 667 continue write(*,*) "Error creating NetCDF output file. " // trim(adjustl(errmsg)) - call util_exit(FAILURE) + call base_util_exit(FAILURE) end subroutine swiftest_io_netcdf_initialize_output @@ -945,7 +945,7 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier ! Return integer(I4B) :: ierr !! Error code: returns 0 if the read is successful ! Internals - integer(I4B) :: i, idmax, npl_check, ntp_check, nplm_check, str_max, status + integer(I4B) :: i, idmax, npl_check, ntp_check, nplm_check, str_max, status, npl, ntp real(DP), dimension(:), allocatable :: rtemp real(DP), dimension(:,:), allocatable :: vectemp integer(I4B), dimension(:), allocatable :: itemp @@ -955,11 +955,13 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier call nc%open(param, readonly=.true.) call nc%find_tslot(self%t) call self%read_hdr(nc, param) - associate(cb => self%cb, pl => self%pl, tp => self%tp, npl => self%pl%nbody, ntp => self%tp%nbody, tslot => nc%tslot) + associate(cb => self%cb, pl => self%pl, tp => self%tp,tslot => nc%tslot) + ! Save these values as variables as they get reset by the setup method + npl = pl%nbody + ntp = tp%nbody call pl%setup(npl, param) call tp%setup(ntp, param) - call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%name_dimid, len=idmax), "netcdf_io_read_frame_system nf90_inquire_dimension name_dimid" ) allocate(rtemp(idmax)) allocate(vectemp(NDIM,idmax)) @@ -971,7 +973,7 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier if (tslot > nc%max_tslot) then write(*,*) write(*,*) "Error in reading frame from NetCDF file. Requested time index value (",tslot,") exceeds maximum time in file (",nc%max_tslot,"). " - call util_exit(FAILURE) + call base_util_exit(FAILURE) end if call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%str_dimid, len=str_max), "netcdf_io_read_frame_system nf90_inquire_dimension str_dimid" ) @@ -997,19 +999,19 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier if (npl_check /= npl) then write(*,*) "Error reading in NetCDF file: The recorded value of npl does not match the number of active massive bodies" - call util_exit(failure) + call base_util_exit(failure) end if if (ntp_check /= ntp) then write(*,*) "Error reading in NetCDF file: The recorded value of ntp does not match the number of active test particles" - call util_exit(failure) + call base_util_exit(failure) end if if (param%lmtiny_pl) then nplm_check = count(pack(rtemp,plmask) > param%GMTINY ) if (nplm_check /= pl%nplm) then write(*,*) "Error reading in NetCDF file: The recorded value of nplm does not match the number of active fully interacting massive bodies" - call util_exit(failure) + call base_util_exit(failure) end if end if @@ -2675,7 +2677,7 @@ module subroutine swiftest_io_read_in_cb(self, param) 667 continue write(*,*) "Error reading central body file: " // trim(adjustl(errmsg)) - call util_exit(FAILURE) + call base_util_exit(FAILURE) end subroutine swiftest_io_read_in_cb @@ -2714,7 +2716,7 @@ module subroutine swiftest_io_read_in_system(self, param) end if ierr = self%read_frame(tmp_param%system_history%nc, tmp_param) deallocate(tmp_param) - if (ierr /=0) call util_exit(FAILURE) + if (ierr /=0) call base_util_exit(FAILURE) end if param%loblatecb = ((self%cb%j2rp2 /= 0.0_DP) .or. (self%cb%j4rp4 /= 0.0_DP)) @@ -2828,7 +2830,7 @@ module function swiftest_io_read_frame_body(self, iu, param) result(ierr) class default write(*,*) "Error reading body file: " // trim(adjustl(errmsg)) end select - call util_exit(FAILURE) + call base_util_exit(FAILURE) end function swiftest_io_read_frame_body @@ -2859,7 +2861,7 @@ module subroutine swiftest_io_read_in_param(self, param_file_name) 667 continue write(self%display_unit,*) "Error reading parameter file: " // trim(adjustl(errmsg)) - call util_exit(FAILURE) + call base_util_exit(FAILURE) end subroutine swiftest_io_read_in_param @@ -2891,7 +2893,7 @@ module subroutine swiftest_io_set_display_param(self, display_style) self%log_output = .true. case default write(*,*) display_style, " is an unknown display style" - call util_exit(USAGE) + call base_util_exit(USAGE) end select self%display_style = display_style @@ -2900,7 +2902,7 @@ module subroutine swiftest_io_set_display_param(self, display_style) 667 continue write(*,*) "Error opening swiftest log file: " // trim(adjustl(errmsg)) - call util_exit(FAILURE) + call base_util_exit(FAILURE) end subroutine swiftest_io_set_display_param @@ -3002,7 +3004,7 @@ module subroutine swiftest_io_write_frame_system(self, param) 667 continue write(*,*) "Error writing nbody_system frame: " // trim(adjustl(errmsg)) - call util_exit(FAILURE) + call base_util_exit(FAILURE) end subroutine swiftest_io_write_frame_system end submodule s_swiftest_io diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index c5b1602d5..33d620022 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -61,6 +61,7 @@ module swiftest class(swiftest_netcdf_parameters), allocatable :: nc !! NetCDF object attached to this storage object contains procedure :: dump => swiftest_io_dump_storage !! Dumps storage object contents to file + procedure :: dealloc => swiftest_util_dealloc_storage !! Resets a storage object by deallocating all items and resetting the frame counter to 0 procedure :: get_index_values => swiftest_util_get_vals_storage !! Gets the unique values of the indices of a storage object (i.e. body id or time value) procedure :: make_index_map => swiftest_util_index_map_storage !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id procedure :: take_snapshot => swiftest_util_snapshot_system !! Takes a snapshot of the nbody_system for later file storage @@ -70,13 +71,15 @@ module swiftest ! The following extended types or their children should be used, where possible, as the base of any types defined in additional modules, such as new integrators. type, extends(base_parameters) :: swiftest_parameters - type(swiftest_storage(nframes=:)), allocatable :: system_history + class(swiftest_storage), allocatable :: system_history contains + procedure :: dealloc => swiftest_util_dealloc_param procedure :: dump => swiftest_io_dump_param procedure :: reader => swiftest_io_param_reader procedure :: writer => swiftest_io_param_writer procedure :: read_in => swiftest_io_read_in_param procedure :: set_display => swiftest_io_set_display_param + final :: swiftest_final_param end type swiftest_parameters @@ -87,7 +90,7 @@ module swiftest integer(I4B), dimension(:), allocatable :: child !! Index of children particles contains procedure :: dealloc => swiftest_util_dealloc_kin !! Deallocates all allocatable arrays - final :: swiftest_final_kin !! Finalizes the Swiftest kinship object - deallocates all allocatables + final :: swiftest_final_kin !! Finalizes the Swiftest kinship object - deallocates all allocatables end type swiftest_kinship @@ -209,7 +212,8 @@ module swiftest real(DP) :: R0 = 0.0_DP !! Initial radius of the central body real(DP) :: dR = 0.0_DP !! Change in the radius of the central body contains - procedure :: read_in => swiftest_io_read_in_cb !! Read in central body initial conditions from an ASCII file + procedure :: dealloc => swiftest_util_dealloc_cb !! Deallocates all allocatables and resets all values to defaults + procedure :: read_in => swiftest_io_read_in_cb !! Read in central body initial conditions from an ASCII file procedure :: write_frame => swiftest_io_netcdf_write_frame_cb !! I/O routine for writing out a single frame of time-series data for all bodies in the system in NetCDF format procedure :: write_info => swiftest_io_netcdf_write_info_cb !! Dump contents of particle information metadata to file end type swiftest_cb @@ -244,11 +248,11 @@ module swiftest contains ! Massive body-specific concrete methods ! These are concrete because they are the same implemenation for all integrators - procedure :: make_impactors => swiftest_make_impactors_pl !! Make impactors out of the current kinship relationships + procedure :: make_impactors => swiftest_util_make_impactors_pl !! Make impactors out of the current kinship relationships procedure :: discard => swiftest_discard_pl !! Placeholder method for discarding massive bodies procedure :: accel_int => swiftest_kick_getacch_int_pl !! Compute direct cross (third) term heliocentric accelerations of massive bodies procedure :: accel_obl => swiftest_obl_acc_pl !! Compute the barycentric accelerations of bodies due to the oblateness of the central body - procedure :: setup => swiftest_util_setup_pl !! A base constructor that sets the number of bodies and allocates and initializes all arrays + procedure :: setup => swiftest_util_setup_pl !! A base constructor that sets the number of bodies and allocates and initializes all arrays ! procedure :: accel_tides => tides_kick_getacch_pl !! Compute the accelerations of bodies due to tidal interactions with the central body procedure :: append => swiftest_util_append_pl !! Appends elements from one structure to another procedure :: h2b => swiftest_util_coord_h2b_pl !! Convert massive bodies from heliocentric to barycentric coordinates (position and velocity) @@ -309,21 +313,21 @@ module swiftest !! This superclass contains a minimial nbody_system of a set of test particles (tp), massive bodies (pl), and a central body (cb) !! The full swiftest_nbody_system type that is used as the parent class of all integrators is defined in collision - class(swiftest_cb), allocatable :: cb !! Central body data structure - class(swiftest_pl), allocatable :: pl !! Massive body data structure - class(swiftest_tp), allocatable :: tp !! Test particle data structure + class(swiftest_cb), allocatable :: cb !! Central body data structure + class(swiftest_pl), allocatable :: pl !! Massive body data structure + class(swiftest_tp), allocatable :: tp !! Test particle data structure - class(swiftest_tp), allocatable :: tp_discards !! Discarded test particle data structure - class(swiftest_pl), allocatable :: pl_discards !! Discarded massive body particle data structure - class(swiftest_pl), allocatable :: pl_adds !! List of added bodies in mergers or collisions - class(swiftest_tp), allocatable :: tp_adds !! List of added bodies in mergers or collisions - class(encounter_list), allocatable :: pltp_encounter !! List of massive body-test particle encounters in a single step - class(encounter_list), allocatable :: plpl_encounter !! List of massive body-massive body encounters in a single step - class(collision_list_plpl), allocatable :: plpl_collision !! List of massive body-massive body collisions in a single step - class(collision_list_plpl), allocatable :: pltp_collision !! List of massive body-massive body collisions in a single step - class(collision_basic), allocatable :: collider !! Collision system object - class(encounter_storage(nframes=:)), allocatable :: encounter_history !! Stores encounter history for later retrieval and saving to file - class(collision_storage(nframes=:)), allocatable :: collision_history !! Stores encounter history for later retrieval and saving to file + class(swiftest_tp), allocatable :: tp_discards !! Discarded test particle data structure + class(swiftest_pl), allocatable :: pl_discards !! Discarded massive body particle data structure + class(swiftest_pl), allocatable :: pl_adds !! List of added bodies in mergers or collisions + class(swiftest_tp), allocatable :: tp_adds !! List of added bodies in mergers or collisions + class(encounter_list), allocatable :: pltp_encounter !! List of massive body-test particle encounters in a single step + class(encounter_list), allocatable :: plpl_encounter !! List of massive body-massive body encounters in a single step + class(collision_list_plpl), allocatable :: plpl_collision !! List of massive body-massive body collisions in a single step + class(collision_list_plpl), allocatable :: pltp_collision !! List of massive body-massive body collisions in a single step + class(collision_basic), allocatable :: collider !! Collision system object + class(encounter_storage), allocatable :: encounter_history !! Stores encounter history for later retrieval and saving to file + class(collision_storage), allocatable :: collision_history !! Stores encounter history for later retrieval and saving to file real(DP) :: t = -1.0_DP !! Integration current time real(DP) :: GMtot = 0.0_DP !! Total nbody_system mass - used for barycentric coordinate conversion @@ -351,20 +355,20 @@ module swiftest real(DP) :: E_untracked = 0.0_DP !! Energy gained from nbody_system due to escaped bodies ! Energy, momentum, and mass errors (used in error reporting) - real(DP) :: ke_orbit_error = 0.0_DP - real(DP) :: ke_spin_error = 0.0_DP - real(DP) :: pe_error = 0.0_DP - real(DP) :: be_error = 0.0_DP + real(DP) :: ke_orbit_error = 0.0_DP + real(DP) :: ke_spin_error = 0.0_DP + real(DP) :: pe_error = 0.0_DP + real(DP) :: be_error = 0.0_DP real(DP) :: E_orbit_error = 0.0_DP - real(DP) :: Ecoll_error = 0.0_DP + real(DP) :: Ecoll_error = 0.0_DP real(DP) :: E_untracked_error = 0.0_DP - real(DP) :: te_error = 0.0_DP + real(DP) :: te_error = 0.0_DP real(DP) :: L_orbit_error = 0.0_DP real(DP) :: L_spin_error = 0.0_DP real(DP) :: L_escape_error = 0.0_DP - real(DP) :: L_total_error = 0.0_DP - real(DP) :: Mtot_error = 0.0_DP - real(DP) :: Mescape_error = 0.0_DP + real(DP) :: L_total_error = 0.0_DP + real(DP) :: Mtot_error = 0.0_DP + real(DP) :: Mescape_error = 0.0_DP logical :: lbeg !! True if this is the beginning of a step. This is used so that test particle steps can be calculated !! separately from massive bodies. Massive body variables are saved at half steps, and passed to @@ -379,25 +383,26 @@ module swiftest procedure :: conservation_report => swiftest_io_conservation_report !! Compute energy and momentum and print out the change with time procedure :: display_run_information => swiftest_io_display_run_information !! Displays helpful information about the run procedure :: dump => swiftest_io_dump_system !! Dump the state of the nbody_system to a file - procedure :: get_t0_values => swiftest_io_netcdf_get_t0_values_system !! Validates the dump file to check whether the dump file initial conditions duplicate the last frame of the netcdf output. + procedure :: get_t0_values => swiftest_io_netcdf_get_t0_values_system !! Validates the dump file to check whether the dump file initial conditions duplicate the last frame of the netcdf output. procedure :: read_frame => swiftest_io_netcdf_read_frame_system !! Read in a frame of input data from file procedure :: write_frame_netcdf => swiftest_io_netcdf_write_frame_system !! Write a frame of input data from file - procedure :: write_frame_system => swiftest_io_write_frame_system !! Write a frame of input data from file procedure :: read_hdr => swiftest_io_netcdf_read_hdr_system !! Read a header for an output frame in NetCDF format procedure :: write_hdr => swiftest_io_netcdf_write_hdr_system !! Write a header for an output frame in NetCDF format - procedure :: read_in => swiftest_io_read_in_system !! Reads the initial conditions for an nbody system procedure :: read_particle_info => swiftest_io_netcdf_read_particle_info_system !! Read in particle metadata from file + procedure :: read_in => swiftest_io_read_in_system !! Reads the initial conditions for an nbody system + procedure :: write_discard => swiftest_io_write_discard !! Write out information about discarded and merged planets and test particles in SyMBA + procedure :: write_frame_system => swiftest_io_write_frame_system !! Write a frame of input data from file procedure :: obl_pot => swiftest_obl_pot_system !! Compute the contribution to the total gravitational potential due solely to the oblateness of the central body + procedure :: dealloc => swiftest_util_dealloc_system !! Deallocates all allocatables and resets all values to defaults. Acts as a base for a finalizer + procedure :: get_energy_and_momentum => swiftest_util_get_energy_and_momentum_system !! Calculates the total nbody_system energy and momentum + procedure :: get_idvals => swiftest_util_get_idvalues_system !! Returns an array of all id values in use in the nbody_system + procedure :: rescale => swiftest_util_rescale_system !! Rescales the nbody_system into a new set of units procedure :: initialize => swiftest_util_setup_initialize_system !! Initialize the nbody_system from input files procedure :: init_particle_info => swiftest_util_setup_initialize_particle_info_system !! Initialize the nbody_system from input files - ! procedure :: step_spin => tides_step_spin_system !! Steps the spins of the massive & central bodies due to tides. - procedure :: set_msys => swiftest_util_set_msys !! Sets the value of msys from the masses of nbody_system bodies. - procedure :: get_energy_and_momentum => swiftest_util_get_energy_and_momentum_system !! Calculates the total nbody_system energy and momentum - procedure :: get_idvals => swiftest_util_get_idvalues_system !! Returns an array of all id values in use in the nbody_system - procedure :: rescale => swiftest_util_rescale_system !! Rescales the nbody_system into a new set of units - procedure :: validate_ids => swiftest_util_valid_id_system !! Validate the numerical ids passed to the nbody_system and save the maximum value - procedure :: write_discard => swiftest_io_write_discard !! Write out information about discarded and merged planets and test particles in SyMBA - generic :: write_frame => write_frame_system, write_frame_netcdf !! Generic method call for reading a frame of output data + ! procedure :: step_spin => tides_step_spin_system !! Steps the spins of the massive & central bodies due to tides. + procedure :: set_msys => swiftest_util_set_msys !! Sets the value of msys from the masses of nbody_system bodies. + procedure :: validate_ids => swiftest_util_valid_id_system !! Validate the numerical ids passed to the nbody_system and save the maximum value + generic :: write_frame => write_frame_system, write_frame_netcdf !! Generic method call for reading a frame of output data end type swiftest_nbody_system @@ -597,7 +602,7 @@ end subroutine swiftest_io_dump_system module subroutine swiftest_io_dump_storage(self, param) implicit none - class(swiftest_storage(*)), intent(inout) :: self !! Swiftest simulation history storage object + class(swiftest_storage), intent(inout) :: self !! Swiftest simulation history storage object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine swiftest_io_dump_storage @@ -1137,7 +1142,6 @@ module subroutine swiftest_util_append_body(self, source, lsource_mask) logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine swiftest_util_append_body - module subroutine swiftest_util_append_pl(self, source, lsource_mask) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object @@ -1235,11 +1239,31 @@ module subroutine swiftest_util_dealloc_kin(self) class(swiftest_kinship), intent(inout) :: self !! Swiftest kinship object end subroutine swiftest_util_dealloc_kin + module subroutine swiftest_util_dealloc_param(self) + implicit none + class(swiftest_parameters),intent(inout) :: self !! Collection of parameters + end subroutine swiftest_util_dealloc_param + + module subroutine swiftest_util_dealloc_cb(self) + implicit none + class(swiftest_cb), intent(inout) :: self !! Swiftest central body object + end subroutine swiftest_util_dealloc_cb + module subroutine swiftest_util_dealloc_pl(self) implicit none class(swiftest_pl), intent(inout) :: self end subroutine swiftest_util_dealloc_pl + module subroutine swiftest_util_dealloc_storage(self) + implicit none + class(swiftest_storage), intent(inout) :: self !! Swiftest storage object + end subroutine swiftest_util_dealloc_storage + + module subroutine swiftest_util_dealloc_system(self) + implicit none + class(swiftest_nbody_system), intent(inout) :: self + end subroutine swiftest_util_dealloc_system + module subroutine swiftest_util_dealloc_tp(self) implicit none class(swiftest_tp), intent(inout) :: self @@ -1390,7 +1414,7 @@ end subroutine swiftest_util_get_potential_energy_triangular interface module subroutine swiftest_util_get_vals_storage(self, idvals, tvals) - class(swiftest_storage(*)), intent(in) :: self !! Swiftest storage object + class(swiftest_storage), intent(in) :: self !! Swiftest storage object integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values in all snapshots real(DP), dimension(:), allocatable, intent(out) :: tvals !! Array of all time values in all snapshots end subroutine swiftest_util_get_vals_storage @@ -1403,9 +1427,26 @@ end subroutine swiftest_util_index_array module subroutine swiftest_util_index_map_storage(self) implicit none - class(swiftest_storage(*)), intent(inout) :: self !! Swiftest storage object + class(swiftest_storage), intent(inout) :: self !! Swiftest storage object end subroutine swiftest_util_index_map_storage + module subroutine swiftest_util_make_impactors_pl(self, idx) + implicit none + class(swiftest_pl), intent(inout) :: self !! Massive body object + integer(I4B), dimension(:), intent(in) :: idx !! Array holding the indices of the two bodies involved in the collision) + end subroutine swiftest_util_make_impactors_pl + + module subroutine swiftest_util_peri(n,m, r, v, atp, q, isperi) + implicit none + integer(I4B), intent(in) :: n !! Number of bodies + real(DP), dimension(:), intent(in) :: m !! Mass term (mu for HELIO coordinates, and Gmtot for BARY) + real(DP), dimension(:,:), intent(in) :: r !! Position vectors (rh for HELIO coordinates, rb for BARY) + real(DP), dimension(:,:), intent(in) :: v !! Position vectors (vh for HELIO coordinates, rb for BARY) + real(DP), dimension(:), intent(out) :: atp !! Semimajor axis + real(DP), dimension(:), intent(out) :: q !! Periapsis + integer(I4B), dimension(:), intent(inout) :: isperi !! Periapsis passage flag + end subroutine swiftest_util_peri + module subroutine swiftest_util_peri_body(self, nbody_system, param) implicit none class(swiftest_body), intent(inout) :: self !! SyMBA massive body object @@ -1440,10 +1481,6 @@ module subroutine swiftest_util_reset_kinship_pl(self, idx) integer(I4B), dimension(:), intent(in) :: idx !! Index array of bodies to reset end subroutine swiftest_util_reset_kinship_pl - module subroutine swiftest_util_reset_storage(self) - implicit none - class(swiftest_storage(*)), intent(inout) :: self !! Swiftest storage object - end subroutine swiftest_util_reset_storage end interface @@ -1504,6 +1541,14 @@ module subroutine swiftest_util_resize_pl(self, nnew) integer(I4B), intent(in) :: nnew !! New size neded end subroutine swiftest_util_resize_pl + module subroutine swiftest_util_resize_storage(storage, nold, nnew) + use base, only : base_storage + implicit none + class(base_storage), allocatable, intent(inout) :: storage !! Original storage object + integer(I4B), intent(in) :: nold !! Old size + integer(I4B), intent(in) :: nnew !! New size + end subroutine swiftest_util_resize_storage + module subroutine swiftest_util_resize_tp(self, nnew) implicit none class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object @@ -1582,9 +1627,16 @@ module subroutine swiftest_util_set_rhill_approximate(self,cb) class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object end subroutine swiftest_util_set_rhill_approximate + module subroutine swiftest_util_snapshot_save(storage, snapshot) + use base, only : base_storage + implicit none + class(base_storage), allocatable, intent(inout) :: storage !! Storage ncounter storage object + class(*), intent(in) :: snapshot !! Object to snapshot + end subroutine swiftest_util_snapshot_save + module subroutine swiftest_util_snapshot_system(self, param, nbody_system, t, arg) implicit none - class(swiftest_storage(*)), intent(inout) :: self !! Swiftest storage object + class(swiftest_storage), intent(inout) :: self !! Swiftest storage object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object to store real(DP), intent(in), optional :: t !! Time of snapshot if different from nbody_system time @@ -1869,32 +1921,32 @@ end subroutine swiftest_util_version end interface contains - subroutine swiftest_make_impactors_pl(self, idx) + subroutine swiftest_final_kin(self) !! author: David A. Minton !! - !! This is a simple wrapper function that is used to make a type-bound procedure using a subroutine whose interface is in the collision module, which must be defined first + !! Finalize the swiftest kinship object - deallocates all allocatables implicit none - class(swiftest_pl), intent(inout) :: self !! Massive body object - integer(I4B), dimension(:), intent(in) :: idx !! Array holding the indices of the two bodies involved in the collision) - - call collision_resolve_make_impactors_pl(self, idx) - + ! Argument + type(swiftest_kinship), intent(inout) :: self !! SyMBA kinship object + + call self%dealloc() + return - end subroutine swiftest_make_impactors_pl + end subroutine swiftest_final_kin - subroutine swiftest_final_kin(self) + subroutine swiftest_final_param(self) !! author: David A. Minton !! - !! Finalize the swiftest kinship object - deallocates all allocatables + !! Finalize the Swiftest parameter object - deallocates all allocatables implicit none ! Argument - type(swiftest_kinship), intent(inout) :: self !! SyMBA kinship object + type(swiftest_parameters), intent(inout) :: self !! SyMBA kinship object call self%dealloc() return - end subroutine swiftest_final_kin + end subroutine swiftest_final_param subroutine swiftest_final_netcdf_parameters(self) @@ -1917,34 +1969,12 @@ subroutine swiftest_final_storage(self) !! Finalizer for the storage data type implicit none ! Arguments - type(swiftest_storage(*)) :: self - ! Internals - integer(I4B) :: i - - do i = 1, self%nframes - if (allocated(self%frame(i)%item)) deallocate(self%frame(i)%item) - end do + type(swiftest_storage) :: self + + call self%dealloc() return end subroutine swiftest_final_storage - - subroutine swiftest_final_system(self) - !! author: David A. Minton - !! - !! Finalize the swiftest nbody system object - deallocates all allocatables - implicit none - ! Argument - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object - - if (allocated(self%cb)) deallocate(self%cb) - if (allocated(self%pl)) deallocate(self%pl) - if (allocated(self%tp)) deallocate(self%tp) - if (allocated(self%tp_discards)) deallocate(self%tp_discards) - if (allocated(self%pl_discards)) deallocate(self%pl_discards) - - return - end subroutine swiftest_final_system - end module swiftest diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 093756af3..062dd1ff9 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -307,7 +307,7 @@ module subroutine swiftest_util_append_pl(self, source, lsource_mask) end associate class default write(*,*) "Invalid object passed to the append method. Source must be of class swiftest_pl or its descendents" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select return @@ -334,7 +334,7 @@ module subroutine swiftest_util_append_tp(self, source, lsource_mask) end associate class default write(*,*) "Invalid object passed to the append method. Source must be of class swiftest_tp or its descendents" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select return @@ -707,18 +707,21 @@ end subroutine swiftest_util_copy_particle_info_arr module subroutine swiftest_util_dealloc_body(self) !! author: David A. Minton !! - !! Finalize the swiftest body object - deallocates all allocatables + !! Finalize the Swiftest body object - deallocates all allocatables implicit none ! Argument class(swiftest_body), intent(inout) :: self - if (allocated(self%info)) deallocate(self%info) + self%nbody = 0 + self%lfirst = .true. + if (allocated(self%id)) deallocate(self%id) + if (allocated(self%info)) deallocate(self%info) if (allocated(self%status)) deallocate(self%status) + if (allocated(self%lmask)) deallocate(self%lmask) if (allocated(self%ldiscard)) deallocate(self%ldiscard) if (allocated(self%lcollision)) deallocate(self%lcollision) if (allocated(self%lencounter)) deallocate(self%lencounter) - if (allocated(self%lmask)) deallocate(self%lmask) if (allocated(self%mu)) deallocate(self%mu) if (allocated(self%rh)) deallocate(self%rh) if (allocated(self%vh)) deallocate(self%vh) @@ -729,9 +732,11 @@ module subroutine swiftest_util_dealloc_body(self) if (allocated(self%agr)) deallocate(self%agr) if (allocated(self%atide)) deallocate(self%atide) if (allocated(self%ir3h)) deallocate(self%ir3h) + if (allocated(self%isperi)) deallocate(self%isperi) + if (allocated(self%peri)) deallocate(self%peri) + if (allocated(self%atp)) deallocate(self%atp) if (allocated(self%a)) deallocate(self%a) if (allocated(self%e)) deallocate(self%e) - if (allocated(self%e)) deallocate(self%e) if (allocated(self%inc)) deallocate(self%inc) if (allocated(self%capom)) deallocate(self%capom) if (allocated(self%omega)) deallocate(self%omega) @@ -741,6 +746,48 @@ module subroutine swiftest_util_dealloc_body(self) end subroutine swiftest_util_dealloc_body + module subroutine swiftest_util_dealloc_cb(self) + !! author: David A. Minton + !! + !! Finalize the Swiftest central body object - deallocates all allocatables + implicit none + ! Arguments + class(swiftest_cb), intent(inout) :: self !! Swiftest central body object + + if (allocated(self%info)) deallocate(self%info) + + self%id = 0 + self%mass = 0.0_DP + self%Gmass = 0.0_DP + self%radius = 0.0_DP + self%density = 1.0_DP + self%j2rp2 = 0.0_DP + self%j4rp4 = 0.0_DP + self%aobl = 0.0_DP + self%atide = 0.0_DP + self%aoblbeg = 0.0_DP + self%aoblend = 0.0_DP + self%atidebeg = 0.0_DP + self%atideend = 0.0_DP + self%rb = 0.0_DP + self%vb = 0.0_DP + self%agr = 0.0_DP + self%Ip = 0.0_DP + self%rot = 0.0_DP + self%k2 = 0.0_DP + self%Q = 0.0_DP + self%tlag = 0.0_DP + self%L0 = 0.0_DP + self%dL = 0.0_DP + self%GM0 = 0.0_DP + self%dGM = 0.0_DP + self%R0 = 0.0_DP + self%dR = 0.0_DP + + return + end subroutine swiftest_util_dealloc_cb + + module subroutine swiftest_util_dealloc_kin(self) !! author: David A. Minton !! @@ -755,10 +802,29 @@ module subroutine swiftest_util_dealloc_kin(self) end subroutine swiftest_util_dealloc_kin + module subroutine swiftest_util_dealloc_param(self) + !! author: David A. Minton + !! + !! Deallocates all allocatables + implicit none + ! Arguments + class(swiftest_parameters),intent(inout) :: self !! Collection of parameters + + if (allocated(self%system_history)) then + call self%system_history%dealloc() + deallocate(self%system_history) + end if + + call base_util_dealloc_param(self) + + return + end subroutine swiftest_util_dealloc_param + + module subroutine swiftest_util_dealloc_pl(self) !! author: David A. Minton !! - !! Finalize the swiftest massive body object - deallocates all allocatables + !! Finalize the Swiftest massive body object - deallocates all allocatables implicit none ! Argument class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object @@ -771,8 +837,11 @@ module subroutine swiftest_util_dealloc_pl(self) if (allocated(self%renc)) deallocate(self%renc) if (allocated(self%radius)) deallocate(self%radius) if (allocated(self%density)) deallocate(self%density) - if (allocated(self%rot)) deallocate(self%rot) + if (allocated(self%rbeg)) deallocate(self%rbeg) + if (allocated(self%rend)) deallocate(self%rend) + if (allocated(self%vbeg)) deallocate(self%vbeg) if (allocated(self%Ip)) deallocate(self%Ip) + if (allocated(self%rot)) deallocate(self%rot) if (allocated(self%k2)) deallocate(self%k2) if (allocated(self%Q)) deallocate(self%Q) if (allocated(self%tlag)) deallocate(self%tlag) @@ -795,16 +864,99 @@ module subroutine swiftest_util_dealloc_pl(self) end subroutine swiftest_util_dealloc_pl + module subroutine swiftest_util_dealloc_storage(self) + !! author: David A. Minton + !! + !! Resets a storage object by deallocating all items and resetting the frame counter to 0 + use base, only : base_util_dealloc_storage + implicit none + ! Arguments + class(swiftest_storage), intent(inout) :: self !! Swiftest storage object + + if (allocated(self%nc)) deallocate(self%nc) + call base_util_dealloc_storage(self) + + return + end subroutine swiftest_util_dealloc_storage + + + module subroutine swiftest_util_dealloc_system(self) + !! author: David A. Minton + !! + !! Deallocates all allocatables and resets all values to defaults. Acts as a base for a finalizer + implicit none + ! Arguments + class(swiftest_nbody_system), intent(inout) :: self + + if (allocated(self%cb)) deallocate(self%cb) + if (allocated(self%pl)) deallocate(self%pl) + if (allocated(self%tp)) deallocate(self%tp) + if (allocated(self%tp_discards)) deallocate(self%tp_discards) + if (allocated(self%pl_discards)) deallocate(self%pl_discards) + if (allocated(self%pl_adds)) deallocate(self%pl_adds) + if (allocated(self%tp_adds)) deallocate(self%tp_adds) + if (allocated(self%pltp_encounter)) deallocate(self%pltp_encounter) + if (allocated(self%plpl_encounter)) deallocate(self%plpl_encounter) + if (allocated(self%plpl_collision)) deallocate(self%plpl_collision) + if (allocated(self%pltp_collision)) deallocate(self%pltp_collision) + if (allocated(self%collider)) deallocate(self%collider) + if (allocated(self%encounter_history)) deallocate(self%encounter_history) + if (allocated(self%collision_history)) deallocate(self%collision_history) + + self%t = -1.0_DP + self%GMtot = 0.0_DP + self%ke_orbit = 0.0_DP + self%ke_spin = 0.0_DP + self%pe = 0.0_DP + self%be = 0.0_DP + self%te = 0.0_DP + self%oblpot = 0.0_DP + self%L_orbit = 0.0_DP + self%L_spin = 0.0_DP + self%L_total = 0.0_DP + self%ke_orbit_orig = 0.0_DP + self%ke_spin_orig = 0.0_DP + self%pe_orig = 0.0_DP + self%be_orig = 0.0_DP + self%E_orbit_orig = 0.0_DP + self%GMtot_orig = 0.0_DP + self%L_total_orig = 0.0_DP + self%L_orbit_orig = 0.0_DP + self%L_spin_orig = 0.0_DP + self%L_escape = 0.0_DP + self%GMescape = 0.0_DP + self%E_collisions = 0.0_DP + self%E_untracked = 0.0_DP + + self%ke_orbit_error = 0.0_DP + self%ke_spin_error = 0.0_DP + self%pe_error = 0.0_DP + self%be_error = 0.0_DP + self%E_orbit_error = 0.0_DP + self%Ecoll_error = 0.0_DP + self%E_untracked_error = 0.0_DP + self%te_error = 0.0_DP + self%L_orbit_error = 0.0_DP + self%L_spin_error = 0.0_DP + self%L_escape_error = 0.0_DP + self%L_total_error = 0.0_DP + self%Mtot_error = 0.0_DP + self%Mescape_error = 0.0_DP + + return + end subroutine swiftest_util_dealloc_system + + module subroutine swiftest_util_dealloc_tp(self) !! author: David A. Minton !! - !! Finalize the swiftest test particle object - deallocates all allocatables + !! Finalize the Swiftest test particle object - deallocates all allocatables implicit none ! Argument class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object - if (allocated(self%nplenc)) deallocate(self%nplenc) if (allocated(self%k_pltp)) deallocate(self%k_pltp) + if (allocated(self%nplenc)) deallocate(self%nplenc) call swiftest_util_dealloc_body(self) @@ -1468,7 +1620,7 @@ module subroutine swiftest_util_get_vals_storage(self, idvals, tvals) !! !! Gets the id values in a storage object, regardless of whether it is encounter of collision ! Argument - class(swiftest_storage(*)), intent(in) :: self !! Swiftest storage object + class(swiftest_storage), intent(in) :: self !! Swiftest storage object integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values in all snapshots real(DP), dimension(:), allocatable, intent(out) :: tvals !! Array of all time values in all snapshots ! Internals @@ -1560,7 +1712,7 @@ module subroutine swiftest_util_index_map_storage(self) !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id implicit none ! Arguments - class(swiftest_storage(*)), intent(inout) :: self !! Swiftest storage object + class(swiftest_storage), intent(inout) :: self !! Swiftest storage object ! Internals integer(I4B), dimension(:), allocatable :: idvals real(DP), dimension(:), allocatable :: tvals @@ -1576,7 +1728,22 @@ module subroutine swiftest_util_index_map_storage(self) return end subroutine swiftest_util_index_map_storage - subroutine swiftest_util_peri(n,m, r, v, atp, q, isperi) + + module subroutine swiftest_util_make_impactors_pl(self, idx) + !! author: David A. Minton + !! + !! This is a simple wrapper function that is used to make a type-bound procedure using a subroutine whose interface is in the collision module, which must be defined first + implicit none + class(swiftest_pl), intent(inout) :: self !! Massive body object + integer(I4B), dimension(:), intent(in) :: idx !! Array holding the indices of the two bodies involved in the collision) + + call collision_resolve_make_impactors_pl(self, idx) + + return + end subroutine swiftest_util_make_impactors_pl + + + module subroutine swiftest_util_peri(n,m, r, v, atp, q, isperi) !! author: David A. Minton !! !! Helper function that does the pericenter passage computation for any body @@ -1685,13 +1852,6 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param) deallocate(lmask) end if - ! Store the original plplenc list so we don't remove any of the original encounters - nenc_old = nbody_system%plpl_encounter%nenc - if (nenc_old > 0) then - allocate(plplenc_old, source=nbody_system%plpl_encounter) - call plplenc_old%copy(nbody_system%plpl_encounter) - end if - ! Add in any new bodies if (nadd > 0) then ! Append the adds to the main pl object @@ -1724,83 +1884,92 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param) ! Reset the kinship trackers call pl%reset_kinship([(i, i=1, npl)]) - ! Re-build the encounter list - ! Be sure to get the level info if this is a SyMBA nbody_system - select type(nbody_system) - class is (symba_nbody_system) - select type(pl) - class is (symba_pl) - select type(tp) - class is (symba_tp) - lencounter = pl%encounter_check(param, nbody_system, param%dt, nbody_system%irec) - if (tp%nbody > 0) then - lencounter = tp%encounter_check(param, nbody_system, param%dt, nbody_system%irec) + if (allocated(nbody_system%plpl_encounter)) then + ! Store the original plplenc list so we don't remove any of the original encounters + nenc_old = nbody_system%plpl_encounter%nenc + if (nenc_old > 0) then + allocate(plplenc_old, source=nbody_system%plpl_encounter) + call plplenc_old%copy(nbody_system%plpl_encounter) end if - end select - end select - end select - ! Re-index the encounter list as the index values may have changed - if (nenc_old > 0) then - nencmin = min(nbody_system%plpl_encounter%nenc, plplenc_old%nenc) - nbody_system%plpl_encounter%nenc = nencmin - do k = 1, nencmin - idnew1 = nbody_system%plpl_encounter%id1(k) - idnew2 = nbody_system%plpl_encounter%id2(k) - idold1 = plplenc_old%id1(k) - idold2 = plplenc_old%id2(k) - if ((idnew1 == idold1) .and. (idnew2 == idold2)) then - ! This is an encounter we already know about, so save the old information - nbody_system%plpl_encounter%lvdotr(k) = plplenc_old%lvdotr(k) - nbody_system%plpl_encounter%lclosest(k) = plplenc_old%lclosest(k) - nbody_system%plpl_encounter%status(k) = plplenc_old%status(k) - nbody_system%plpl_encounter%r1(:,k) = plplenc_old%r1(:,k) - nbody_system%plpl_encounter%r2(:,k) = plplenc_old%r2(:,k) - nbody_system%plpl_encounter%v1(:,k) = plplenc_old%v1(:,k) - nbody_system%plpl_encounter%v2(:,k) = plplenc_old%v2(:,k) - nbody_system%plpl_encounter%tcollision(k) = plplenc_old%tcollision(k) - nbody_system%plpl_encounter%level(k) = plplenc_old%level(k) - else if (((idnew1 == idold2) .and. (idnew2 == idold1))) then - ! This is an encounter we already know about, but with the order reversed, so save the old information - nbody_system%plpl_encounter%lvdotr(k) = plplenc_old%lvdotr(k) - nbody_system%plpl_encounter%lclosest(k) = plplenc_old%lclosest(k) - nbody_system%plpl_encounter%status(k) = plplenc_old%status(k) - nbody_system%plpl_encounter%r1(:,k) = plplenc_old%r2(:,k) - nbody_system%plpl_encounter%r2(:,k) = plplenc_old%r1(:,k) - nbody_system%plpl_encounter%v1(:,k) = plplenc_old%v2(:,k) - nbody_system%plpl_encounter%v2(:,k) = plplenc_old%v1(:,k) - nbody_system%plpl_encounter%tcollision(k) = plplenc_old%tcollision(k) - nbody_system%plpl_encounter%level(k) = plplenc_old%level(k) + ! Re-build the encounter list + ! Be sure to get the level info if this is a SyMBA nbody_system + select type(nbody_system) + class is (symba_nbody_system) + select type(pl) + class is (symba_pl) + select type(tp) + class is (symba_tp) + lencounter = pl%encounter_check(param, nbody_system, param%dt, nbody_system%irec) + if (tp%nbody > 0) then + lencounter = tp%encounter_check(param, nbody_system, param%dt, nbody_system%irec) end if - nbody_system%plpl_encounter%index1(k) = findloc(pl%id(1:npl), nbody_system%plpl_encounter%id1(k), dim=1) - nbody_system%plpl_encounter%index2(k) = findloc(pl%id(1:npl), nbody_system%plpl_encounter%id2(k), dim=1) - end do - if (allocated(lmask)) deallocate(lmask) - allocate(lmask(nencmin)) - nenc_old = nencmin - if (any(nbody_system%plpl_encounter%index1(1:nencmin) == 0) .or. any(nbody_system%plpl_encounter%index2(1:nencmin) == 0)) then - lmask(:) = nbody_system%plpl_encounter%index1(1:nencmin) /= 0 .and. nbody_system%plpl_encounter%index2(1:nencmin) /= 0 - else - return - end if - nencmin = count(lmask(:)) - nbody_system%plpl_encounter%nenc = nencmin - if (nencmin > 0) then - nbody_system%plpl_encounter%index1(1:nencmin) = pack(nbody_system%plpl_encounter%index1(1:nenc_old), lmask(1:nenc_old)) - nbody_system%plpl_encounter%index2(1:nencmin) = pack(nbody_system%plpl_encounter%index2(1:nenc_old), lmask(1:nenc_old)) - nbody_system%plpl_encounter%id1(1:nencmin) = pack(nbody_system%plpl_encounter%id1(1:nenc_old), lmask(1:nenc_old)) - nbody_system%plpl_encounter%id2(1:nencmin) = pack(nbody_system%plpl_encounter%id2(1:nenc_old), lmask(1:nenc_old)) - nbody_system%plpl_encounter%lvdotr(1:nencmin) = pack(nbody_system%plpl_encounter%lvdotr(1:nenc_old), lmask(1:nenc_old)) - nbody_system%plpl_encounter%lclosest(1:nencmin) = pack(nbody_system%plpl_encounter%lclosest(1:nenc_old), lmask(1:nenc_old)) - nbody_system%plpl_encounter%status(1:nencmin) = pack(nbody_system%plpl_encounter%status(1:nenc_old), lmask(1:nenc_old)) - nbody_system%plpl_encounter%tcollision(1:nencmin) = pack(nbody_system%plpl_encounter%tcollision(1:nenc_old), lmask(1:nenc_old)) - nbody_system%plpl_encounter%level(1:nencmin) = pack(nbody_system%plpl_encounter%level(1:nenc_old), lmask(1:nenc_old)) - do i = 1, NDIM - nbody_system%plpl_encounter%r1(i, 1:nencmin) = pack(nbody_system%plpl_encounter%r1(i, 1:nenc_old), lmask(1:nenc_old)) - nbody_system%plpl_encounter%r2(i, 1:nencmin) = pack(nbody_system%plpl_encounter%r2(i, 1:nenc_old), lmask(1:nenc_old)) - nbody_system%plpl_encounter%v1(i, 1:nencmin) = pack(nbody_system%plpl_encounter%v1(i, 1:nenc_old), lmask(1:nenc_old)) - nbody_system%plpl_encounter%v2(i, 1:nencmin) = pack(nbody_system%plpl_encounter%v2(i, 1:nenc_old), lmask(1:nenc_old)) + end select + end select + end select + + ! Re-index the encounter list as the index values may have changed + if (nenc_old > 0) then + nencmin = min(nbody_system%plpl_encounter%nenc, plplenc_old%nenc) + nbody_system%plpl_encounter%nenc = nencmin + do k = 1, nencmin + idnew1 = nbody_system%plpl_encounter%id1(k) + idnew2 = nbody_system%plpl_encounter%id2(k) + idold1 = plplenc_old%id1(k) + idold2 = plplenc_old%id2(k) + if ((idnew1 == idold1) .and. (idnew2 == idold2)) then + ! This is an encounter we already know about, so save the old information + nbody_system%plpl_encounter%lvdotr(k) = plplenc_old%lvdotr(k) + nbody_system%plpl_encounter%lclosest(k) = plplenc_old%lclosest(k) + nbody_system%plpl_encounter%status(k) = plplenc_old%status(k) + nbody_system%plpl_encounter%r1(:,k) = plplenc_old%r1(:,k) + nbody_system%plpl_encounter%r2(:,k) = plplenc_old%r2(:,k) + nbody_system%plpl_encounter%v1(:,k) = plplenc_old%v1(:,k) + nbody_system%plpl_encounter%v2(:,k) = plplenc_old%v2(:,k) + nbody_system%plpl_encounter%tcollision(k) = plplenc_old%tcollision(k) + nbody_system%plpl_encounter%level(k) = plplenc_old%level(k) + else if (((idnew1 == idold2) .and. (idnew2 == idold1))) then + ! This is an encounter we already know about, but with the order reversed, so save the old information + nbody_system%plpl_encounter%lvdotr(k) = plplenc_old%lvdotr(k) + nbody_system%plpl_encounter%lclosest(k) = plplenc_old%lclosest(k) + nbody_system%plpl_encounter%status(k) = plplenc_old%status(k) + nbody_system%plpl_encounter%r1(:,k) = plplenc_old%r2(:,k) + nbody_system%plpl_encounter%r2(:,k) = plplenc_old%r1(:,k) + nbody_system%plpl_encounter%v1(:,k) = plplenc_old%v2(:,k) + nbody_system%plpl_encounter%v2(:,k) = plplenc_old%v1(:,k) + nbody_system%plpl_encounter%tcollision(k) = plplenc_old%tcollision(k) + nbody_system%plpl_encounter%level(k) = plplenc_old%level(k) + end if + nbody_system%plpl_encounter%index1(k) = findloc(pl%id(1:npl), nbody_system%plpl_encounter%id1(k), dim=1) + nbody_system%plpl_encounter%index2(k) = findloc(pl%id(1:npl), nbody_system%plpl_encounter%id2(k), dim=1) end do + if (allocated(lmask)) deallocate(lmask) + allocate(lmask(nencmin)) + nenc_old = nencmin + if (any(nbody_system%plpl_encounter%index1(1:nencmin) == 0) .or. any(nbody_system%plpl_encounter%index2(1:nencmin) == 0)) then + lmask(:) = nbody_system%plpl_encounter%index1(1:nencmin) /= 0 .and. nbody_system%plpl_encounter%index2(1:nencmin) /= 0 + else + return + end if + nencmin = count(lmask(:)) + nbody_system%plpl_encounter%nenc = nencmin + if (nencmin > 0) then + nbody_system%plpl_encounter%index1(1:nencmin) = pack(nbody_system%plpl_encounter%index1(1:nenc_old), lmask(1:nenc_old)) + nbody_system%plpl_encounter%index2(1:nencmin) = pack(nbody_system%plpl_encounter%index2(1:nenc_old), lmask(1:nenc_old)) + nbody_system%plpl_encounter%id1(1:nencmin) = pack(nbody_system%plpl_encounter%id1(1:nenc_old), lmask(1:nenc_old)) + nbody_system%plpl_encounter%id2(1:nencmin) = pack(nbody_system%plpl_encounter%id2(1:nenc_old), lmask(1:nenc_old)) + nbody_system%plpl_encounter%lvdotr(1:nencmin) = pack(nbody_system%plpl_encounter%lvdotr(1:nenc_old), lmask(1:nenc_old)) + nbody_system%plpl_encounter%lclosest(1:nencmin) = pack(nbody_system%plpl_encounter%lclosest(1:nenc_old), lmask(1:nenc_old)) + nbody_system%plpl_encounter%status(1:nencmin) = pack(nbody_system%plpl_encounter%status(1:nenc_old), lmask(1:nenc_old)) + nbody_system%plpl_encounter%tcollision(1:nencmin) = pack(nbody_system%plpl_encounter%tcollision(1:nenc_old), lmask(1:nenc_old)) + nbody_system%plpl_encounter%level(1:nencmin) = pack(nbody_system%plpl_encounter%level(1:nenc_old), lmask(1:nenc_old)) + do i = 1, NDIM + nbody_system%plpl_encounter%r1(i, 1:nencmin) = pack(nbody_system%plpl_encounter%r1(i, 1:nenc_old), lmask(1:nenc_old)) + nbody_system%plpl_encounter%r2(i, 1:nencmin) = pack(nbody_system%plpl_encounter%r2(i, 1:nenc_old), lmask(1:nenc_old)) + nbody_system%plpl_encounter%v1(i, 1:nencmin) = pack(nbody_system%plpl_encounter%v1(i, 1:nenc_old), lmask(1:nenc_old)) + nbody_system%plpl_encounter%v2(i, 1:nencmin) = pack(nbody_system%plpl_encounter%v2(i, 1:nenc_old), lmask(1:nenc_old)) + end do + end if end if end if end associate @@ -1871,7 +2040,6 @@ module subroutine swiftest_util_reset_kinship_pl(self, idx) ! Internals integer(I4B) :: i, j - self%kin(idx(:))%parent = idx(:) self%kin(idx(:))%nchild = 0 do j = 1, size(idx(:)) @@ -2601,7 +2769,7 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param) write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled' case default write(*,*) 'Unkown integrator',param%integrator - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select allocate(swiftest_particle_info :: nbody_system%cb%info) @@ -2657,9 +2825,13 @@ module subroutine swiftest_util_setup_initialize_system(self, param) class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody_system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - allocate(swiftest_storage(param%dump_cadence) :: param%system_history) + if (allocated(param%system_history)) then + call param%system_history%dealloc() + deallocate(param%system_history) + end if + allocate(swiftest_storage :: param%system_history) + call param%system_history%setup(param%dump_cadence) allocate(swiftest_netcdf_parameters :: param%system_history%nc) - call param%system_history%reset() associate(nbody_system => self, cb => self%cb, pl => self%pl, tp => self%tp, nc => param%system_history%nc) call nbody_system%read_in(param) @@ -2874,7 +3046,7 @@ module subroutine swiftest_util_setup_tp(self, n, param) return end subroutine swiftest_util_setup_tp - + module subroutine swiftest_util_snapshot_system(self, param, nbody_system, t, arg) !! author: David A. Minton @@ -2882,19 +3054,68 @@ module subroutine swiftest_util_snapshot_system(self, param, nbody_system, t, ar !! Takes a snapshot of the nbody_system for later file storage implicit none ! Arguments - class(swiftest_storage(*)), intent(inout) :: self !! Swiftest storage object + class(swiftest_storage), intent(inout) :: self !! Swiftest storage object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object to store real(DP), intent(in), optional :: t !! Time of snapshot if different from nbody_system time character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) - - self%iframe = self%iframe + 1 + ! Internals + class(swiftest_nbody_system), allocatable :: snapshot + + ! Take a minimal snapshot wihout all of the extra storage objects + allocate(snapshot, mold=nbody_system) + allocate(snapshot%cb, source=nbody_system%cb ) + allocate(snapshot%pl, source=nbody_system%pl ) + allocate(snapshot%tp, source=nbody_system%tp ) + + snapshot%t = nbody_system%t + snapshot%GMtot = nbody_system%GMtot + snapshot%ke_orbit = nbody_system%ke_orbit + snapshot%ke_spin = nbody_system%ke_spin + snapshot%pe = nbody_system%pe + snapshot%be = nbody_system%be + snapshot%te = nbody_system%te + snapshot%oblpot = nbody_system%oblpot + snapshot%L_orbit = nbody_system%L_orbit + snapshot%L_spin = nbody_system%L_spin + snapshot%L_total = nbody_system%L_total + snapshot%ke_orbit_orig = nbody_system%ke_orbit_orig + snapshot%ke_spin_orig = nbody_system%ke_spin_orig + snapshot%pe_orig = nbody_system%pe_orig + snapshot%be_orig = nbody_system%be_orig + snapshot%E_orbit_orig = nbody_system%E_orbit_orig + snapshot%GMtot_orig = nbody_system%GMtot_orig + snapshot%L_total_orig = nbody_system%L_total_orig + snapshot%L_orbit_orig = nbody_system%L_orbit_orig + snapshot%L_spin_orig = nbody_system%L_spin_orig + snapshot%L_escape = nbody_system%L_escape + snapshot%GMescape = nbody_system%GMescape + snapshot%E_collisions = nbody_system%E_collisions + snapshot%E_untracked = nbody_system%E_untracked + snapshot%ke_orbit_error = nbody_system%ke_orbit_error + snapshot%ke_spin_error = nbody_system%ke_spin_error + snapshot%pe_error = nbody_system%pe_error + snapshot%be_error = nbody_system%be_error + snapshot%E_orbit_error = nbody_system%E_orbit_error + snapshot%Ecoll_error = nbody_system%Ecoll_error + snapshot%E_untracked_error = nbody_system%E_untracked_error + snapshot%te_error = nbody_system%te_error + snapshot%L_orbit_error = nbody_system%L_orbit_error + snapshot%L_spin_error = nbody_system%L_spin_error + snapshot%L_escape_error = nbody_system%L_escape_error + snapshot%L_total_error = nbody_system%L_total_error + snapshot%Mtot_error = nbody_system%Mtot_error + snapshot%Mescape_error = nbody_system%Mescape_error + snapshot%lbeg = nbody_system%lbeg + + + ! Store a snapshot of the nbody_system for posterity + call base_util_snapshot_save(self, snapshot) self%nt = self%iframe - self%frame(self%iframe) = nbody_system ! Store a snapshot of the nbody_system for posterity self%nid = self%nid + 1 ! Central body if (allocated(nbody_system%pl)) self%nid = self%nid + nbody_system%pl%nbody if (allocated(nbody_system%tp)) self%nid = self%nid + nbody_system%tp%nbody - + return end subroutine swiftest_util_snapshot_system @@ -4509,7 +4730,7 @@ module subroutine swiftest_util_valid_id_system(self, param) if (idarr(i) == idarr(i+1)) then write(*, *) "Swiftest error:" write(*, *) " more than one body/particle has id = ", idarr(i) - call util_exit(FAILURE) + call base_util_exit(FAILURE) end if end do param%maxid = max(param%maxid, maxval(idarr)) diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index c454fffb0..ca7761ff0 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -19,7 +19,7 @@ module function symba_encounter_check_pl(self, param, nbody_system, dt, irec) re implicit none ! Arguments class(symba_pl), intent(inout) :: self !! SyMBA test particle object - class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters class(symba_nbody_system), intent(inout) :: nbody_system !! SyMBA nbody system object real(DP), intent(in) :: dt !! step size integer(I4B), intent(in) :: irec !! Current recursion level @@ -90,7 +90,7 @@ end function symba_encounter_check_pl module function symba_encounter_check_list_plpl(self, param, nbody_system, dt, irec) result(lany_encounter) implicit none class(symba_list_plpl), intent(inout) :: self !! SyMBA pl-pl encounter list object - class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters class(symba_nbody_system), intent(inout) :: nbody_system !! SyMBA nbody system object real(DP), intent(in) :: dt !! step size integer(I4B), intent(in) :: irec !! Current recursion level @@ -158,7 +158,7 @@ end function symba_encounter_check_list_plpl module function symba_encounter_check_list_pltp(self, param, nbody_system, dt, irec) result(lany_encounter) implicit none class(symba_list_pltp), intent(inout) :: self !! SyMBA pl-tp encounter list object - class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters class(symba_nbody_system), intent(inout) :: nbody_system !! SyMBA nbody system object real(DP), intent(in) :: dt !! step size integer(I4B), intent(in) :: irec !! Current recursion level @@ -234,7 +234,7 @@ module function symba_encounter_check_tp(self, param, nbody_system, dt, irec) re implicit none ! Arguments class(symba_tp), intent(inout) :: self !! SyMBA test particle object - class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters class(symba_nbody_system), intent(inout) :: nbody_system !! SyMBA nbody system object real(DP), intent(in) :: dt !! step size integer(I4B), intent(in) :: irec !! Current recursion level diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 5ab8480d1..05256309a 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -21,7 +21,7 @@ module subroutine symba_kick_getacch_int_pl(self, param) implicit none ! Arguments class(symba_pl), intent(inout) :: self !! SyMBA massive body object - class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameter + class(swiftest_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter ! Internals ! type(interaction_timer), save :: itimer ! logical, save :: lfirst = .true. diff --git a/src/symba/symba_module.f90 b/src/symba/symba_module.f90 index bdc7011e3..efa8dda8f 100644 --- a/src/symba/symba_module.f90 +++ b/src/symba/symba_module.f90 @@ -48,7 +48,6 @@ module symba procedure :: sort => symba_util_sort_pl !! Sorts body arrays by a sortable componen procedure :: rearrange => symba_util_sort_rearrange_pl !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods procedure :: spill => symba_util_spill_pl !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) - final :: symba_final_pl !! Finalizes the SyMBA massive body object - deallocates all allocatables end type symba_pl @@ -69,7 +68,6 @@ module symba procedure :: sort => symba_util_sort_tp !! Sorts body arrays by a sortable componen procedure :: rearrange => symba_util_sort_rearrange_tp !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods procedure :: spill => symba_util_spill_tp !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) - final :: symba_final_tp !! Finalizes the SyMBA test particle object - deallocates all allocatables end type symba_tp @@ -78,7 +76,6 @@ module symba contains procedure :: encounter_check => symba_encounter_check_list_plpl !! Checks if massive bodies are going through close encounters with each other procedure :: kick => symba_kick_list_plpl !! Kick barycentric velocities of active massive bodies within SyMBA recursion - final :: symba_final_list_plpl !! Finalizes the SyMBA test particle object - deallocates all allocatables end type symba_list_plpl @@ -87,20 +84,19 @@ module symba contains procedure :: encounter_check => symba_encounter_check_list_pltp !! Checks if massive bodies are going through close encounters with test particles procedure :: kick => symba_kick_list_pltp !! Kick barycentric velocities of active test particles within SyMBA recursion - final :: symba_final_list_pltp !! Finalizes the SyMBA test particle object - deallocates all allocatables end type symba_list_pltp type, extends(helio_nbody_system) :: symba_nbody_system integer(I4B) :: irec = -1 !! nbody_system recursion level contains - procedure :: initialize => symba_util_setup_initialize_system !! Performs SyMBA-specific initilization steps + procedure :: dealloic => symba_util_dealloc_system !! Deallocates all allocatables + procedure :: initialize => symba_util_setup_initialize_system !! Performs SyMBA-specific initilization steps procedure :: step => symba_step_system !! Advance the SyMBA nbody system forward in time by one step procedure :: interp => symba_step_interp_system !! Perform an interpolation step on the SymBA nbody system procedure :: set_recur_levels => symba_step_set_recur_levels_system !! Sets recursion levels of bodies and encounter lists to the current nbody_system level procedure :: recursive_step => symba_step_recur_system !! Step interacting planets and active test particles ahead in democratic heliocentric coordinates at the current recursion level, if applicable, and descend to the next deeper level if necessary procedure :: reset => symba_step_reset_system !! Resets pl, tp,and encounter structures at the start of a new step - final :: symba_final_system !! Finalizes the SyMBA nbody system object - deallocates all allocatables end type symba_nbody_system interface @@ -130,7 +126,7 @@ end subroutine symba_drift_tp module function symba_encounter_check_pl(self, param, nbody_system, dt, irec) result(lany_encounter) implicit none class(symba_pl), intent(inout) :: self !! SyMBA test particle object - class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters class(symba_nbody_system), intent(inout) :: nbody_system !! SyMBA nbody system object real(DP), intent(in) :: dt !! step size integer(I4B), intent(in) :: irec !! Current recursion level @@ -140,7 +136,7 @@ end function symba_encounter_check_pl module function symba_encounter_check_list_plpl(self, param, nbody_system, dt, irec) result(lany_encounter) implicit none class(symba_list_plpl), intent(inout) :: self !! SyMBA pl-pl encounter list object - class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters class(symba_nbody_system), intent(inout) :: nbody_system !! SyMBA nbody system object real(DP), intent(in) :: dt !! step size integer(I4B), intent(in) :: irec !! Current recursion level @@ -150,7 +146,7 @@ end function symba_encounter_check_list_plpl module function symba_encounter_check_list_pltp(self, param, nbody_system, dt, irec) result(lany_encounter) implicit none class(symba_list_pltp), intent(inout) :: self !! SyMBA pl-tp encounter list object - class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters class(symba_nbody_system), intent(inout) :: nbody_system !! SyMBA nbody system object real(DP), intent(in) :: dt !! step size integer(I4B), intent(in) :: irec !! Current recursion level @@ -160,7 +156,7 @@ end function symba_encounter_check_list_pltp module function symba_encounter_check_tp(self, param, nbody_system, dt, irec) result(lany_encounter) implicit none class(symba_tp), intent(inout) :: self !! SyMBA test particle object - class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters class(symba_nbody_system), intent(inout) :: nbody_system !! SyMBA nbody system object real(DP), intent(in) :: dt !! step size integer(I4B), intent(in) :: irec !! Current recursion level @@ -203,7 +199,7 @@ end subroutine symba_io_param_writer module subroutine symba_kick_getacch_int_pl(self, param) implicit none class(symba_pl), intent(inout) :: self !! SyMBA massive body object - class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters end subroutine symba_kick_getacch_int_pl module subroutine symba_kick_getacch_pl(self, nbody_system, param, t, lbeg) @@ -242,6 +238,11 @@ module subroutine symba_kick_list_pltp(self, nbody_system, dt, irec, sgn) integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration end subroutine symba_kick_list_pltp + module subroutine symba_util_dealloc_system(self) + implicit none + class(symba_nbody_system), intent(inout) :: self + end subroutine symba_util_dealloc_system + module subroutine symba_util_setup_initialize_system(self, param) implicit none class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody_system object @@ -407,79 +408,4 @@ module subroutine symba_util_spill_tp(self, discards, lspill_list, ldestructive) end subroutine symba_util_spill_tp end interface - - contains - - - subroutine symba_final_list_plpl(self) - !! author: David A. Minton - !! - !! Finalize the pl-tp list - deallocates all allocatables - implicit none - type(symba_list_plpl), intent(inout) :: self !! SyMBA encounter list object - - call self%dealloc() - - return - end subroutine symba_final_list_plpl - - - subroutine symba_final_list_pltp(self) - !! author: David A. Minton - !! - !! Finalize the pl-tp list - deallocates all allocatables - implicit none - type(symba_list_pltp), intent(inout) :: self !! SyMBA encounter list object - - call self%dealloc() - - return - end subroutine symba_final_list_pltp - - - subroutine symba_final_pl(self) - !! author: David A. Minton - !! - !! Finalize the SyMBA massive body object - deallocates all allocatables - implicit none - ! Argument - type(symba_pl), intent(inout) :: self !! SyMBA massive body object - - call self%dealloc() - - return - end subroutine symba_final_pl - - - subroutine symba_final_system(self) - !! author: David A. Minton - !! - !! Finalize the SyMBA nbody system object - deallocates all allocatables - implicit none - ! Argument - type(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - - if (allocated(self%pl_adds)) deallocate(self%pl_adds) - if (allocated(self%pltp_encounter)) deallocate(self%pltp_encounter) - if (allocated(self%plpl_encounter)) deallocate(self%plpl_encounter) - if (allocated(self%plpl_collision)) deallocate(self%plpl_collision) - - call helio_final_system(self%helio_nbody_system) - - return - end subroutine symba_final_system - - - subroutine symba_final_tp(self) - !! author: David A. Minton - !! - !! Finalize the SyMBA test particleobject - deallocates all allocatables - implicit none - ! Argument - type(symba_tp), intent(inout) :: self !! SyMBA test particle object - - call self%dealloc() - - return - end subroutine symba_final_tp end module symba \ No newline at end of file diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 5d2737234..c811b1968 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -198,7 +198,7 @@ recursive module subroutine symba_step_recur_system(self, param, t, ireci) 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) + call base_util_exit(FAILURE) END IF irecp = ireci + 1 if (ireci == 0) then diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index c3dc1be18..ccda74e01 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -33,7 +33,7 @@ module subroutine symba_util_append_pl(self, source, lsource_mask) end associate class default write(*,*) "Invalid object passed to the append method. Source must be of class symba_pl or its descendents!" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select return @@ -61,7 +61,7 @@ module subroutine symba_util_append_tp(self, source, lsource_mask) end associate class default write(*,*) "Invalid object passed to the append method. Source must be of class symba_tp or its descendents!" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select return @@ -87,6 +87,21 @@ module subroutine symba_util_dealloc_pl(self) end subroutine symba_util_dealloc_pl + module subroutine symba_util_dealloc_system(self) + !! author: David A. Minton + !! + !! Deallocates all allocatables and resets all values to defaults. Acts as a base for a finalizer + implicit none + ! Arguments + class(symba_nbody_system), intent(inout) :: self + + self%irec = -1 + call self%helio_nbody_system%dealloc() + + return + end subroutine symba_util_dealloc_system + + module subroutine symba_util_dealloc_tp(self) !! author: David A. Minton !! @@ -125,7 +140,7 @@ module subroutine symba_util_fill_pl(self, inserts, lfill_list) call swiftest_util_fill_pl(keeps, inserts, lfill_list) ! Note: helio_pl does not have its own fill method, so we skip back to the base class class default write(*,*) "Invalid object passed to the fill method. Source must be of class symba_pl or its descendents!" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select end associate @@ -155,7 +170,7 @@ module subroutine symba_util_fill_tp(self, inserts, lfill_list) call swiftest_util_fill_tp(keeps, inserts, lfill_list) ! Note: helio_tp does not have its own fill method, so we skip back to the base class class default write(*,*) "Invalid object passed to the fill method. Source must be of class symba_tp or its descendents!" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select end associate @@ -233,7 +248,6 @@ module subroutine symba_util_resize_tp(self, nnew) return end subroutine symba_util_resize_tp - module subroutine symba_util_set_renc(self, scale) !! author: David A. Minton !! @@ -258,6 +272,7 @@ module subroutine symba_util_set_renc(self, scale) return end subroutine symba_util_set_renc + module subroutine symba_util_setup_initialize_system(self, param) !! author: David A. Minton !! @@ -272,6 +287,8 @@ module subroutine symba_util_setup_initialize_system(self, param) type(encounter_storage) :: encounter_history type(collision_storage) :: collision_history + call encounter_history%setup(4096) + call collision_history%setup(4096) ! Call parent method associate(nbody_system => self) call helio_util_setup_initialize_system(nbody_system, param) @@ -281,7 +298,6 @@ module subroutine symba_util_setup_initialize_system(self, param) if (param%lenc_save_trajectory .or. param%lenc_save_closest) then allocate(encounter_netcdf_parameters :: encounter_history%nc) - call encounter_history%reset() select type(nc => encounter_history%nc) class is (encounter_netcdf_parameters) nc%file_name = ENCOUNTER_OUTFILE @@ -292,9 +308,8 @@ module subroutine symba_util_setup_initialize_system(self, param) end select allocate(nbody_system%encounter_history, source=encounter_history) end if - + allocate(collision_netcdf_parameters :: collision_history%nc) - call collision_history%reset() select type(nc => collision_history%nc) class is (collision_netcdf_parameters) nc%file_name = COLLISION_OUTFILE @@ -522,7 +537,7 @@ module subroutine symba_util_spill_pl(self, discards, lspill_list, ldestructive) call swiftest_util_spill_pl(keeps, discards, lspill_list, ldestructive) class default write(*,*) "Invalid object passed to the spill method. Source must be of class symba_pl or its descendents!" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select end associate @@ -554,7 +569,7 @@ module subroutine symba_util_spill_tp(self, discards, lspill_list, ldestructive) call swiftest_util_spill_tp(keeps, discards, lspill_list, ldestructive) class default write(*,*) "Invalid object passed to the spill method. Source must be of class symba_tp or its descendents!" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select end associate diff --git a/src/whm/whm_drift.f90 b/src/whm/whm_drift.f90 index 5d507c42b..1161919f2 100644 --- a/src/whm/whm_drift.f90 +++ b/src/whm/whm_drift.f90 @@ -50,7 +50,7 @@ module subroutine whm_drift_pl(self, nbody_system, param, dt) call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) end if end do - call util_exit(FAILURE) + call base_util_exit(FAILURE) end if end associate diff --git a/src/whm/whm_module.f90 b/src/whm/whm_module.f90 index d902065f6..d2d10f971 100644 --- a/src/whm/whm_module.f90 +++ b/src/whm/whm_module.f90 @@ -272,8 +272,6 @@ end subroutine whm_util_sort_rearrange_pl contains - - subroutine whm_final_pl(self) !! author: David A. Minton !! @@ -296,7 +294,7 @@ subroutine whm_final_system(self) ! Arguments type(whm_nbody_system), intent(inout) :: self !! WHM nbody system object - call swiftest_final_system(self) + call self%dealloc() return end subroutine whm_final_system diff --git a/src/whm/whm_util.f90 b/src/whm/whm_util.f90 index a9377536d..751e97d2a 100644 --- a/src/whm/whm_util.f90 +++ b/src/whm/whm_util.f90 @@ -35,7 +35,7 @@ module subroutine whm_util_append_pl(self, source, lsource_mask) end associate class default write(*,*) "Invalid object passed to the append method. Source must be of class whm_pl or its descendents" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select return @@ -87,7 +87,7 @@ module subroutine whm_util_fill_pl(self, inserts, lfill_list) call swiftest_util_fill_pl(keeps, inserts, lfill_list) class default write(*,*) "Invalid object passed to the fill method. Inserts must be of class whm_pl or its descendents!" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select end associate @@ -324,7 +324,7 @@ module subroutine whm_util_spill_pl(self, discards, lspill_list, ldestructive) call swiftest_util_spill_pl(keeps, discards, lspill_list, ldestructive) class default write(*,*) "Invalid object passed to the spill method. Source must be of class whm_pl or its descendents!" - call util_exit(FAILURE) + call base_util_exit(FAILURE) end select end associate