From f7afdf34fa032de0d28ac1e899ad086d329d8974 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 7 Dec 2022 19:32:17 -0500 Subject: [PATCH] Improved the stability of the encounter snapshot --- examples/Fragmentation/Fragmentation_Movie.py | 6 ++-- src/modules/symba_classes.f90 | 17 ++++++++--- src/symba/symba_io.f90 | 10 +++---- src/symba/symba_util.f90 | 30 +++++++++++++++---- 4 files changed, 45 insertions(+), 18 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 329a6f343..ab2e2b34d 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -41,8 +41,8 @@ movie_titles = dict(zip(available_movie_styles, movie_title_list)) # These initial conditions were generated by trial and error -pos_vectors = {"disruption_headon" : [np.array([1.0, -5.0e-05, 0.0]), - np.array([1.0, 5.0e-05 ,0.0])], +pos_vectors = {"disruption_headon" : [np.array([1.0, -5.0e-04, 0.0]), + np.array([1.0, 5.0e-04 ,0.0])], "supercatastrophic_off_axis": [np.array([1.0, -4.2e-05, 0.0]), np.array([1.0, 4.2e-05, 0.0])], "hitandrun" : [np.array([1.0, -2.0e-05, 0.0]), @@ -203,7 +203,7 @@ def data_stream(self, frame=0): minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades sim.set_parameter(fragmentation=True, fragmentation_save="TRAJECTORY", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) - sim.run(dt=1e-5, tstop=4.0e-3, istep_out=10, dump_cadence=0) + sim.run(dt=1e-4, tstop=1.0e-4, istep_out=1, dump_cadence=0) print("Generating animation") anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1) \ No newline at end of file diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 6c3aaae12..4c1b6b673 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -229,11 +229,15 @@ module symba_classes end type symba_nbody_system - type, extends(symba_nbody_system) :: symba_encounter_snapshot - integer(I4B) :: tslot !! The index for the time array in the final NetCDF file + type :: symba_encounter_snapshot + type(symba_pl) :: pl !! Massive body data structure + type(symba_tp) :: tp !! Test particle data structure + real(DP) :: t = 0.0_DP !! Time at the snapshot + integer(I4B) :: tslot = 0 !! The index for the time array in the final NetCDF file contains - procedure :: write_encounter_frame => symba_io_encounter_write_frame !! Writes a frame of encounter data to file - generic :: write_frame => write_encounter_frame + procedure :: write_encounter_frame => symba_io_encounter_write_frame !! Writes a frame of encounter data to file + procedure :: dealloc => symba_util_dealloc_snapshot !! Deallocates all allocatable arrays + generic :: write_frame => write_encounter_frame !! Writes a snaphot frame to file final :: symba_util_final_encounter_snapshot end type symba_encounter_snapshot @@ -657,6 +661,11 @@ module subroutine symba_util_dealloc_merger(self) class(symba_merger), intent(inout) :: self !! SyMBA body merger object end subroutine symba_util_dealloc_merger + module subroutine symba_util_dealloc_snapshot(self) + implicit none + class(symba_encounter_snapshot), intent(inout) :: self !! SyMBA nbody system object + end subroutine symba_util_dealloc_snapshot + module subroutine symba_util_dealloc_system(self) implicit none class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index bc0f1b985..a4e5a9922 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -95,7 +95,8 @@ module subroutine symba_io_encounter_initialize_output(self, param) call check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], nc%ptype_varid), "symba_io_encounter_initialize_output nf90_def_var ptype_varid" ) call check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%rh_varid), "symba_io_encounter_initialize_output nf90_def_var rh_varid" ) call check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%vh_varid), "symba_io_encounter_initialize_output nf90_def_var vh_varid" ) - call check( nf90_def_var(nc%id, nc%gmass_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%Gmass_varid), "symba_io_encounter_initialize_output nf90_def_var Gmass_varid" ) + call check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%Gmass_varid), "symba_io_encounter_initialize_output nf90_def_var Gmass_varid" ) + call check( nf90_def_var(nc%id, nc%level_varname, NF90_INT, [nc%id_dimid, nc%time_dimid], nc%level_varid), "symba_io_encounter_initialize_output nf90_def_var level_varid" ) if (param%lclose) then call check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%radius_varid), "symba_io_encounter_initialize_output nf90_def_var radius_varid" ) end if @@ -153,8 +154,7 @@ module subroutine symba_io_encounter_write_frame(self, nc, param) tslot = self%tslot call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "symba_io_encounter_write_frame nf90_put_var time_varid" ) - select type(pl => self%pl) - class is (symba_pl) + associate(pl => self%pl, tp => self%tp) npl = pl%nbody do i = 1, npl idslot = pl%id(i) @@ -162,6 +162,7 @@ module subroutine symba_io_encounter_write_frame(self, nc, param) call check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var pl rh_varid" ) call check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var pl vh_varid" ) call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[idslot, tslot]), "symba_io_encounter_write_frame nf90_put_var pl Gmass_varid" ) + call check( nf90_put_var(nc%id, nc%level_varid, pl%levelg(i), start=[idslot, tslot]), "symba_io_encounter_write_frame nf90_put_var pl level_varid" ) if (param%lclose) call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[idslot, tslot]), "symba_io_encounter_write_frame nf90_put_var pl radius_varid" ) @@ -175,10 +176,7 @@ module subroutine symba_io_encounter_write_frame(self, nc, param) charstring = trim(adjustl(pl%info(i)%particle_type)) call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "symba_io_encounter_write_frame nf90_put_var pl particle_type_varid" ) end do - - end select - associate(tp => self%tp) ntp = tp%nbody do i = 1, ntp idslot = tp%id(i) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index b96270fc9..b5b74d93e 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -238,10 +238,28 @@ module subroutine symba_util_dealloc_merger(self) end subroutine symba_util_dealloc_merger + module subroutine symba_util_dealloc_snapshot(self) + !! author: David A. Minton + !! + !! Deallocates allocatable arrays in an encounter snapshot + implicit none + ! Arguments + class(symba_encounter_snapshot), intent(inout) :: self !! SyMBA nbody system object + + + call self%pl%dealloc() + call self%tp%dealloc() + self%t = 0.0_DP + self%tslot = 0 + + return + end subroutine symba_util_dealloc_snapshot + + module subroutine symba_util_dealloc_system(self) !! author: David A. Minton !! - !! Deallocates all allocatabale arrays + !! Deallocates all allocatabale arrays in a SyMBA system object implicit none ! Arguments class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object @@ -937,6 +955,8 @@ module subroutine symba_util_resize_storage(self, nnew) nbig = nbig * 2 end do allocate(symba_encounter_storage(nbig) :: tmp) + tmp%tvals(1:nold) = self%encounter_history%tvals(1:nold) + tmp%tvals(nold+1:nbig) = huge(1.0_DP) if (lmalloc) then do i = 1, nold if (allocated(self%encounter_history%frame(i)%item)) call move_alloc(self%encounter_history%frame(i)%item, tmp%frame(i)%item) @@ -1328,8 +1348,6 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) allocate(symba_encounter_snapshot :: snapshot) snapshot%t = t - allocate(symba_pl :: snapshot%pl) - allocate(symba_tp :: snapshot%tp) if (npl + ntp == 0) return npl_snap = npl ntp_snap = ntp @@ -1353,6 +1371,8 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) allocate(snapshot%pl%id(npl_snap)) allocate(snapshot%pl%info(npl_snap)) allocate(snapshot%pl%Gmass(npl_snap)) + allocate(snapshot%pl%levelg(npl_snap)) + snapshot%pl%levelg(:) = pack(pl%levelg(1:npl), pl%lmask(1:npl)) snapshot%pl%id(:) = pack(pl%id(1:npl), pl%lmask(1:npl)) snapshot%pl%info(:) = pack(pl%info(1:npl), pl%lmask(1:npl)) snapshot%pl%Gmass(:) = pack(pl%Gmass(1:npl), pl%lmask(1:npl)) @@ -1394,8 +1414,6 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) end if ! Save the snapshot - self%encounter_history%iframe = self%encounter_history%iframe + 1 - call self%resize_storage(self%encounter_history%iframe) ! 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 @@ -1403,6 +1421,8 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) if (t <= self%encounter_history%tvals(i)) then snapshot%tslot = i self%encounter_history%tvals(i) = t + self%encounter_history%iframe = i + call self%resize_storage(i) exit end if end do