diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 4e2ada041..a664ef49e 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -28,6 +28,7 @@ import swiftest import numpy as np +import xarray as xr import matplotlib.pyplot as plt import matplotlib.animation as animation from pathlib import Path @@ -75,21 +76,37 @@ for k,v in body_Gmass.items(): body_radius[k] = [((Gmass/GU)/(4./3.*np.pi*density))**(1./3.) for Gmass in v] - # ---------------------------------------------------------------------------------------------------------------------- # Define the animation class that will generate the movies of the fragmentation outcomes # ---------------------------------------------------------------------------------------------------------------------- -figsize = (4,4) + + +def encounter_combiner(sim): + """ + Combines simulation data with encounter data to produce a dataset that contains the position, + mass, radius, etc. of both. It will interpolate over empty time values to fill in gaps. + """ + + # Only keep a minimal subset of necessary data from the simulation and encounter datasets + keep_vars = ['rh','Gmass','radius'] + data = sim.data[keep_vars] + enc = sim.enc[keep_vars].load() + + # Remove any encounter data at the same time steps that appear in the data to prevent duplicates + t_not_duplicate = ~enc['time'].isin(data['time']) + enc = enc.where(t_not_duplicate,drop=True) + + # The following will combine the two datasets along the time dimension, sort the time dimension, and then fill in any time gaps with interpolation + ds = xr.combine_nested([data,enc],concat_dim='time').sortby("time").interpolate_na(dim="time") + + return ds + class AnimatedScatter(object): """An animated scatter plot using matplotlib.animations.FuncAnimation.""" def __init__(self, sim, animfile, title, style, nskip=1): - if style == "disruption_headon" or style == "hitandrun": - tgood = sim.enc.where(sim.enc['Gmass'] > 0.8 * body_Gmass[style][0]).time - else: - tgood = sim.enc.where(sim.enc['Gmass'] > 0.01 * body_Gmass[style][1]).time - self.ds = sim.enc[['rh','vh','Gmass','radius']].sel(time=tgood).load().interpolate_na(dim="time") + self.ds = encounter_combiner(sim) nframes = int(self.ds['time'].size) self.sim = sim @@ -101,6 +118,7 @@ def __init__(self, sim, animfile, title, style, nskip=1): 'Central body': 'xkcd:almost black'} # Set up the figure and axes... + self.figsize = (4,4) self.fig, self.ax = self.setup_plot() # Then setup FuncAnimation. @@ -110,7 +128,7 @@ def __init__(self, sim, animfile, title, style, nskip=1): print(f"Finished writing {animfile}") def setup_plot(self): - fig = plt.figure(figsize=figsize, dpi=300) + fig = plt.figure(figsize=self.figsize, dpi=300) plt.tight_layout(pad=0) # Calculate the distance along the y-axis between the colliding bodies at the start of the simulation. @@ -120,7 +138,7 @@ def setup_plot(self): scale_frame = abs(rhy1) + abs(rhy2) ax = plt.Axes(fig, [0.1, 0.1, 0.8, 0.8]) - self.ax_pt_size = figsize[0] * 0.8 * 72 / (2 * scale_frame) + self.ax_pt_size = self.figsize[0] * 0.8 * 72 / (2 * scale_frame) ax.set_xlim(-scale_frame, scale_frame) ax.set_ylim(-scale_frame, scale_frame) ax.set_xticks([]) @@ -185,6 +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-4, tstop=4.0e-3, istep_out=1, dump_cadence=0) + sim.run(dt=1e-4, tstop=1.0e-3, 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/CMakeLists.txt b/src/CMakeLists.txt index 7f45ddd5a..eb9fb20d5 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -30,6 +30,7 @@ SET(FAST_MATH_FILES ${SRC}/encounter/encounter_check.f90 ${SRC}/encounter/encounter_setup.f90 ${SRC}/encounter/encounter_util.f90 + ${SRC}/encounter/encounter_io.f90 ${SRC}/fraggle/fraggle_generate.f90 ${SRC}/fraggle/fraggle_io.f90 ${SRC}/fraggle/fraggle_placeholder.f90 @@ -74,6 +75,7 @@ SET(FAST_MATH_FILES ${SRC}/util/util_dealloc.f90 ${SRC}/util/util_exit.f90 ${SRC}/util/util_fill.f90 + ${SRC}/util/util_final.f90 ${SRC}/util/util_flatten.f90 ${SRC}/util/util_get_energy_momentum.f90 ${SRC}/util/util_index_array.f90 diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 new file mode 100644 index 000000000..0cc07b009 --- /dev/null +++ b/src/encounter/encounter_io.f90 @@ -0,0 +1,207 @@ +!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh +!! This file is part of Swiftest. +!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +!! You should have received a copy of the GNU General Public License along with Swiftest. +!! If not, see: https://www.gnu.org/licenses. + +submodule (encounter_classes) s_encounter_io + use swiftest +contains + + + module subroutine encounter_io_dump(self, param) + !! author: David A. Minton + !! + !! Dumps the time history of an encounter to file. + implicit none + ! Arguments + class(encounter_storage(*)), intent(inout) :: self !! Encounter storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: i + + ! Most of this is just temporary test code just to get something working. Eventually this should get cleaned up. + + do i = 1, self%nframes + if (allocated(self%frame(i)%item)) then + select type(snapshot => self%frame(i)%item) + class is (encounter_snapshot) + param%ioutput = self%tslot(i) + call snapshot%write_frame(self%nc,param) + end select + else + exit + end if + end do + + + return + end subroutine encounter_io_dump + + + module subroutine encounter_io_initialize(self, param) + !! author: David A. Minton + !! + !! Initialize a NetCDF encounter file system. This is a simplified version of the main simulation output NetCDF file, but with fewer variables. + use, intrinsic :: ieee_arithmetic + use netcdf + implicit none + ! Arguments + class(encounter_io_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: nvar, varid, vartype + real(DP) :: dfill + real(SP) :: sfill + logical :: fileExists + character(len=STRMAX) :: errmsg + integer(I4B) :: ndims + + associate(nc => self) + dfill = ieee_value(dfill, IEEE_QUIET_NAN) + sfill = ieee_value(sfill, IEEE_QUIET_NAN) + + select case (param%out_type) + case("NETCDF_FLOAT") + self%out_type = NF90_FLOAT + case("NETCDF_DOUBLE") + self%out_type = NF90_DOUBLE + end select + + ! Check if the file exists, and if it does, delete it + inquire(file=nc%enc_file, exist=fileExists) + if (fileExists) then + open(unit=LUN, file=nc%enc_file, status="old", err=667, iomsg=errmsg) + close(unit=LUN, status="delete") + end if + + call check( nf90_create(nc%enc_file, NF90_NETCDF4, nc%id), "encounter_io_initialize nf90_create" ) + + ! Dimensions + call check( nf90_def_dim(nc%id, nc%time_dimname, nc%time_dimsize, nc%time_dimid), "encounter_io_initialize nf90_def_dim time_dimid" ) ! Simulation time dimension + call check( nf90_def_dim(nc%id, nc%space_dimname, NDIM , nc%space_dimid), "encounter_io_initialize nf90_def_dim space_dimid" ) ! 3D space dimension + call check( nf90_def_dim(nc%id, nc%id_dimname, param%maxid, nc%id_dimid), "encounter_io_initialize nf90_def_dim id_dimid" ) ! dimension to store particle id numbers + call check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "encounter_io_initialize nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) + + ! Dimension coordinates + call check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, nc%time_dimid, nc%time_varid), "encounter_io_initialize nf90_def_var time_varid" ) + call check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "encounter_io_initialize nf90_def_var space_varid" ) + call check( nf90_def_var(nc%id, nc%id_dimname, NF90_INT, nc%id_dimid, nc%id_varid), "encounter_io_initialize nf90_def_var id_varid" ) + + ! Variables + call check( nf90_def_var(nc%id, nc%name_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], nc%name_varid), "encounter_io_initialize nf90_def_var name_varid" ) + call check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], nc%ptype_varid), "encounter_io_initialize 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), "encounter_io_initialize 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), "encounter_io_initialize 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), "encounter_io_initialize 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), "encounter_io_initialize 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), "encounter_io_initialize nf90_def_var radius_varid" ) + end if + if (param%lrotation) then + call check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%Ip_varid), "encounter_io_initialize nf90_def_var Ip_varid" ) + call check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%rot_varid), "encounter_io_initialize nf90_def_var rot_varid" ) + end if + + call check( nf90_inquire(nc%id, nVariables=nvar), "encounter_io_initialize nf90_inquire nVariables" ) + do varid = 1, nvar + call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "encounter_io_initialize nf90_inquire_variable" ) + select case(vartype) + case(NF90_INT) + call check( nf90_def_var_fill(nc%id, varid, 0, NF90_FILL_INT), "encounter_io_initialize nf90_def_var_fill NF90_INT" ) + case(NF90_FLOAT) + call check( nf90_def_var_fill(nc%id, varid, 0, sfill), "encounter_io_initialize nf90_def_var_fill NF90_FLOAT" ) + case(NF90_DOUBLE) + call check( nf90_def_var_fill(nc%id, varid, 0, dfill), "encounter_io_initialize nf90_def_var_fill NF90_DOUBLE" ) + case(NF90_CHAR) + call check( nf90_def_var_fill(nc%id, varid, 0, 0), "encounter_io_initialize nf90_def_var_fill NF90_CHAR" ) + end select + end do + + ! Take the file out of define mode + call check( nf90_enddef(nc%id), "encounter_io_initialize nf90_enddef" ) + + ! Add in the space dimension coordinates + call check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "encounter_io_initialize nf90_put_var space" ) + end associate + + return + + 667 continue + write(*,*) "Error creating encounter output file. " // trim(adjustl(errmsg)) + call util_exit(FAILURE) + end subroutine encounter_io_initialize + + + module subroutine encounter_io_write_frame(self, nc, param) + !! author: David A. Minton + !! + !! Write a frame of output of an encounter trajectory. + use netcdf + implicit none + ! Arguments + class(encounter_snapshot), intent(in) :: self !! Swiftest encounter structure + class(encounter_io_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: i, tslot, idslot, old_mode, npl, ntp + character(len=NAMELEN) :: charstring + + tslot = param%ioutput + + call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "encounter_io_write_frame nf90_set_fill" ) + + call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "encounter_io_write_frame nf90_put_var time_varid" ) + + associate(pl => self%pl, tp => self%tp) + npl = pl%nbody + do i = 1, npl + idslot = pl%id(i) + call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[idslot]), "encounter_io_write_frame nf90_put_var pl id_varid" ) + call check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_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]), "encounter_io_write_frame nf90_put_var pl vh_varid" ) + call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl Gmass_varid" ) + select type(pl) + class is (symba_pl) + call check( nf90_put_var(nc%id, nc%level_varid, pl%levelg(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl level_varid" ) + end select + + if (param%lclose) call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl radius_varid" ) + + if (param%lrotation) then + call check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl Ip_varid" ) + call check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl rotx_varid" ) + end if + + charstring = trim(adjustl(pl%info(i)%name)) + call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "encounter_io_write_frame nf90_put_var pl name_varid" ) + 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]), "encounter_io_write_frame nf90_put_var pl particle_type_varid" ) + end do + + ntp = tp%nbody + do i = 1, ntp + idslot = tp%id(i) + call check( nf90_put_var(nc%id, nc%id_varid, tp%id(i), start=[idslot]), "encounter_io_write_frame nf90_put_var tp id_varid" ) + call check( nf90_put_var(nc%id, nc%rh_varid, tp%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var tp rh_varid" ) + call check( nf90_put_var(nc%id, nc%vh_varid, tp%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var tp vh_varid" ) + + charstring = trim(adjustl(tp%info(i)%name)) + call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "encounter_io_write_frame nf90_put_var tp name_varid" ) + charstring = trim(adjustl(tp%info(i)%particle_type)) + call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "encounter_io_write_frame nf90_put_var tp particle_type_varid" ) + end do + end associate + + call check( nf90_set_fill(nc%id, old_mode, old_mode) ) + + return + end subroutine encounter_io_write_frame + + + + +end submodule s_encounter_io \ No newline at end of file diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index d0c801ab4..0d3a66d62 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -136,6 +136,36 @@ module subroutine encounter_util_final_list(self) end subroutine encounter_util_final_list + module subroutine encounter_util_final_snapshot(self) + !! author: David A. Minton + !! + !! Deallocates allocatable arrays in an encounter snapshot + implicit none + ! Arguments + type(encounter_snapshot), intent(inout) :: self !! Encounter storage object + + if (allocated(self%pl)) deallocate(self%pl) + if (allocated(self%tp)) deallocate(self%tp) + self%t = 0.0_DP + + return + end subroutine encounter_util_final_snapshot + + + module subroutine encounter_util_final_storage(self) + !! author: David A. Minton + !! + !! Deallocates allocatable arrays in an encounter snapshot + implicit none + ! Arguments + type(encounter_storage(*)), intent(inout) :: self !! Encounter storage object + + call util_final_storage(self%swiftest_storage) + + return + end subroutine encounter_util_final_storage + + module subroutine encounter_util_resize_list(self, nnew) !! author: David A. Minton !! diff --git a/src/fraggle/fraggle_io.f90 b/src/fraggle/fraggle_io.f90 index 87d2d7d11..9d00bfb0f 100644 --- a/src/fraggle/fraggle_io.f90 +++ b/src/fraggle/fraggle_io.f90 @@ -12,6 +12,26 @@ contains + + module subroutine fraggle_io_encounter_dump(self, param) + implicit none + class(fraggle_storage(*)), intent(inout) :: self !! Encounter storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine fraggle_io_encounter_dump + + module subroutine fraggle_io_encounter_initialize_output(self, param) + implicit none + class(fraggle_io_encounter_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(in) :: param + end subroutine fraggle_io_encounter_initialize_output + + module subroutine fraggle_io_encounter_write_frame(self, nc, param) + implicit none + class(fraggle_encounter_snapshot), intent(in) :: self !! Swiftest encounter structure + class(encounter_io_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine fraggle_io_encounter_write_frame + module subroutine fraggle_io_log_generate(frag) !! author: David A. Minton !! @@ -36,31 +56,6 @@ module subroutine fraggle_io_log_generate(frag) write(LUN,fmtlabel) ' Qloss |', -frag%Qloss / abs(frag%Etot_before) write(LUN,fmtlabel) ' dE - Qloss |', (frag%Etot_after - frag%Etot_before + frag%Qloss) / abs(frag%Etot_before) write(LUN, "(' -------------------------------------------------------------------------------------')") - write(LUN, *) "Individual fragment values (collisional system natural units)" - write(LUN, *) "mass" - do i = 1, frag%nbody - write(LUN, *) i, frag%mass(i) - end do - write(LUN, *) "x_coll" - do i = 1, frag%nbody - write(LUN, *) i, frag%x_coll(:,i) - end do - write(LUN, *) "v_coll" - do i = 1, frag%nbody - write(LUN, *) i, frag%v_coll(:,i) - end do - write(LUN, *) "xb" - do i = 1, frag%nbody - write(LUN, *) i, frag%xb(:,i) - end do - write(LUN, *) "vb" - do i = 1, frag%nbody - write(LUN, *) i, frag%vb(:,i) - end do - write(LUN, *) "rot" - do i = 1, frag%nbody - write(LUN, *) i, frag%rot(:,i) - end do close(LUN) @@ -82,69 +77,6 @@ module subroutine fraggle_io_log_pl(pl, param) integer(I4B) :: i character(STRMAX) :: errmsg - open(unit=LUN, file=FRAGGLE_LOG_OUT, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) - write(LUN, *, err = 667, iomsg = errmsg) - - write(LUN, *) "--------------------------------------------------------------------" - write(LUN, *) " Fraggle fragment final body properties" - write(LUN, *) "--------------------------------------------------------------------" - write(LUN, *) "id, name" - do i = 1, pl%nbody - write(LUN, *) i, pl%id(i), pl%info(i)%name - end do - write(LUN, *) "mass, Gmass" - do i = 1, pl%nbody - write(LUN, *) i, pl%mass(i), pl%Gmass(i) - end do - write(LUN, *) "radius" - do i = 1, pl%nbody - write(LUN, *) i, pl%radius(i) - end do - write(LUN, *) "xb" - do i = 1, pl%nbody - write(LUN, *) i, pl%xb(:,i) - end do - write(LUN, *) "vb" - do i = 1, pl%nbody - write(LUN, *) i, pl%vb(:,i) - end do - write(LUN, *) "rh" - do i = 1, pl%nbody - write(LUN, *) i, pl%rh(:,i) - end do - write(LUN, *) "vh" - do i = 1, pl%nbody - write(LUN, *) i, pl%vh(:,i) - end do - - if (param%lrotation) then - write(LUN, *) "rot" - do i = 1, pl%nbody - write(LUN, *) i, pl%rot(:,i) - end do - write(LUN, *) "Ip" - do i = 1, pl%nbody - write(LUN, *) i, pl%Ip(:,i) - end do - end if - - ! if (param%ltides) then - ! write(LUN, *) "Q" - ! do i = 1, pl%nbody - ! write(LUN, *) i, pl%Q(i) - ! end do - ! write(LUN, *) "k2" - ! do i = 1, pl%nbody - ! write(LUN, *) i, pl%k2(i) - ! end do - ! write(LUN, *) "tlag" - ! do i = 1, pl%nbody - ! write(LUN, *) i, pl%tlag(i) - ! end do - ! end if - - close(LUN) - return 667 continue write(*,*) "Error writing Fraggle message to log file: " // trim(adjustl(errmsg)) @@ -167,46 +99,20 @@ module subroutine fraggle_io_log_regime(colliders, frag) write(LUN, *) "--------------------------------------------------------------------" write(LUN, *) " Fraggle collisional regime determination results" write(LUN, *) "--------------------------------------------------------------------" - write(LUN, *) "----------------------- Collider information -----------------------" write(LUN, *) "True number of colliders : ",colliders%ncoll write(LUN, *) "Index list of true colliders : ",colliders%idx(1:colliders%ncoll) - write(LUN, *) "-------------------- Two-body equialent values ---------------------" - write(LUN, *) "mass1 : ",colliders%mass(1) - write(LUN, *) "radius1 : ",colliders%radius(1) - write(LUN, *) "xb1 : ",colliders%xb(:,1) - write(LUN, *) "vb1 : ",colliders%vb(:,1) - write(LUN, *) "rot1 : ",colliders%rot(:,1) - write(LUN, *) "Ip1 : ",colliders%Ip(:,1) - write(LUN, *) "L_spin1 : ",colliders%L_spin(:,1) - write(LUN, *) "L_orbit1 : ",colliders%L_orbit(:,1) - write(LUN, *) "mass2 : ",colliders%mass(2) - write(LUN, *) "radius2 : ",colliders%radius(2) - write(LUN, *) "xb2 : ",colliders%xb(:,2) - write(LUN, *) "vb2 : ",colliders%vb(:,2) - write(LUN, *) "rot2 : ",colliders%rot(:,2) - write(LUN, *) "Ip2 : ",colliders%Ip(:,2) - write(LUN, *) "L_spin2 : ",colliders%L_spin(:,2) - write(LUN, *) "L_orbit2 : ",colliders%L_orbit(:,2) - write(LUN, *) "------------------------------ Regime -----------------------------" select case(frag%regime) case(COLLRESOLVE_REGIME_MERGE) - write(LUN, *) "Merge" + write(LUN, *) "Regime: Merge" case(COLLRESOLVE_REGIME_DISRUPTION) - write(LUN, *) "Disruption" + write(LUN, *) "Regime: Disruption" case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - write(LUN, *) "Supercatastrophic disruption" + write(LUN, *) "Regime: Supercatastrophic disruption" case(COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - write(LUN, *) "Graze and merge" + write(LUN, *) "Regime: Graze and merge" case(COLLRESOLVE_REGIME_HIT_AND_RUN) - write(LUN, *) "Hit and run" + write(LUN, *) "Regime: Hit and run" end select - write(LUN, *) "----------------------- Fragment information ----------------------" - write(LUN, *) "Total mass of fragments : ", frag%mtot - write(LUN, *) "Largest fragment mass : ", frag%mass_dist(1) - write(LUN, *) "Second-largest fragment mass : ", frag%mass_dist(2) - write(LUN, *) "Remaining fragment mass : ", frag%mass_dist(3) - write(LUN, *) "Center of mass position : ", frag%xbcom(:) - write(LUN, *) "Center of mass velocity : ", frag%vbcom(:) write(LUN, *) "Energy loss : ", frag%Qloss write(LUN, *) "--------------------------------------------------------------------" close(LUN) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 5c0803912..8688ec2d9 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -131,6 +131,70 @@ module subroutine fraggle_util_construct_temporary_system(frag, system, param, t end subroutine fraggle_util_construct_temporary_system + + module subroutine fraggle_util_final_colliders(self) + !! author: David A. Minton + !! + !! Finalizer will deallocate all allocatables + implicit none + ! Arguments + type(fraggle_colliders), intent(inout) :: self !! Fraggle encountar storage object + + if (allocated(self%idx)) deallocate(self%idx) + + return + end subroutine fraggle_util_final_colliders + + + module subroutine fraggle_util_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%x_coll)) deallocate(self%x_coll) + if (allocated(self%v_coll)) deallocate(self%v_coll) + if (allocated(self%v_r_unit)) deallocate(self%v_r_unit) + if (allocated(self%v_t_unit)) deallocate(self%v_t_unit) + if (allocated(self%v_n_unit)) deallocate(self%v_n_unit) + if (allocated(self%rmag)) deallocate(self%rmag) + if (allocated(self%rotmag)) deallocate(self%rotmag) + if (allocated(self%v_r_mag)) deallocate(self%v_r_mag) + if (allocated(self%v_t_mag)) deallocate(self%v_t_mag) + + return + end subroutine fraggle_util_final_fragments + + + module subroutine fraggle_util_final_storage(self) + !! author: David A. Minton + !! + !! Finalizer will deallocate all allocatables + implicit none + ! Arguments + type(fraggle_storage(*)), intent(inout) :: self !! Fraggle encountar storage object + + call util_final_storage(self%swiftest_storage) + + return + end subroutine fraggle_util_final_storage + + module subroutine fraggle_util_final_snapshot(self) + !! author: David A. Minton + !! + !! Finalizer will deallocate all allocatables + implicit none + ! Arguments + type(fraggle_encounter_snapshot), intent(inout) :: self !! Fraggle encountar storage object + + call encounter_util_final_snapshot(self%encounter_snapshot) + + return + end subroutine fraggle_util_final_snapshot + + module subroutine fraggle_util_get_energy_momentum(self, colliders, system, param, lbefore) !! Author: David A. Minton !! diff --git a/src/helio/helio_util.f90 b/src/helio/helio_util.f90 index 698bb4849..73a88d58c 100644 --- a/src/helio/helio_util.f90 +++ b/src/helio/helio_util.f90 @@ -33,7 +33,7 @@ module subroutine helio_util_final_system(self) ! Arguments type(helio_nbody_system), intent(inout) :: self !! Helio nbody system object - call self%dealloc() + call whm_util_final_system(self%whm_nbody_system) return end subroutine helio_util_final_system diff --git a/src/io/io.f90 b/src/io/io.f90 index 92924d458..d14f0a694 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -122,9 +122,7 @@ module subroutine io_conservation_report(self, param, lterminal) ! Internals real(DP), dimension(NDIM) :: Ltot_now, Lorbit_now, Lspin_now real(DP) :: ke_orbit_now, ke_spin_now, pe_now, Eorbit_now - real(DP) :: Eorbit_error, Etot_error, Ecoll_error real(DP) :: GMtot_now - real(DP) :: Lerror, Merror character(len=STRMAX) :: errmsg integer(I4B), parameter :: EGYIU = 72 character(len=*), parameter :: EGYTERMFMT = '(" DL/L0 = ", ES12.5 & diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index 4ff06d4a8..6d5abce79 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -28,7 +28,7 @@ program swiftest_driver integer(I4B) :: iout !! Output cadence counter integer(I4B) :: idump !! Dump cadence counter type(walltimer) :: integration_timer !! Object used for computing elapsed wall time - real(DP) :: tfrac + real(DP) :: tfrac !! Fraction of total simulation time completed type(progress_bar) :: pbar !! Object used to print out a progress bar character(*), parameter :: statusfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // & '"; Number of active pl, tp = ", I6, ", ", I6)' @@ -127,7 +127,7 @@ program swiftest_driver if (iout == istep_out) then iout = 0 idump = idump + 1 - system_history%frame(idump) = nbody_system + system_history%frame(idump) = nbody_system ! Store a snapshot of the system for posterity if (idump == dump_cadence) then idump = 0 diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index b04fca1a1..c7edc5e8a 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -38,18 +38,48 @@ module encounter_classes 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. - final :: encounter_util_final_list !! Finalize the encounter list - deallocates all allocatables + final :: encounter_util_final_list !! Finalize the encounter list - deallocates all allocatables end type encounter_list + type :: encounter_snapshot + !! A simplified version of a SyMBA nbody system object for storing minimal snapshots of the system state during encounters + class(swiftest_pl), allocatable :: pl !! Massive body data structure + class(swiftest_tp), allocatable :: tp !! Test particle data structure + real(DP) :: t !! Simulation time when snapshot was taken + contains + procedure :: write_frame => encounter_io_write_frame !! Writes a frame of encounter data to file + final :: encounter_util_final_snapshot + end type encounter_snapshot + + !> NetCDF dimension and variable names for the enounter save object + type, extends(netcdf_parameters) :: encounter_io_parameters + integer(I4B) :: ienc_frame = 1 !! Current frame number for the encounter history + character(STRMAX) :: enc_file !! Encounter output file name + character(NAMELEN) :: level_varname = "level" !! Recursion depth + integer(I4B) :: level_varid !! ID for the recursion level variable + integer(I4B) :: time_dimsize = 0 !! Number of time values in snapshot + integer(I4B) :: id_dimsize = 0 !! Number of potential id values in snapshot + contains + procedure :: initialize => encounter_io_initialize !! Initialize a set of parameters used to identify a NetCDF output object + end type encounter_io_parameters + + !> A class that that is used to store simulation history data between file output + type, extends(swiftest_storage) :: encounter_storage + type(encounter_io_parameters) :: nc !! NetCDF parameter object containing the details about the file attached to this storage object + contains + procedure :: dump => encounter_io_dump !! Dumps contents of encounter history to file + final :: encounter_util_final_storage + end type encounter_storage + type encounter_bounding_box_1D integer(I4B) :: n !! Number of bodies with extents integer(I4B), dimension(:), allocatable :: ind !! Sorted minimum/maximum extent indices (value > n indicates an ending index) integer(I8B), dimension(:), allocatable :: ibeg !! Beginning index for box integer(I8B), dimension(:), allocatable :: iend !! Ending index for box contains - procedure :: sort => encounter_check_sort_aabb_1D !! Sorts the bounding box extents along a single dimension prior to the sweep phase - procedure :: dealloc => encounter_util_dealloc_aabb !! Deallocates all allocatables - final :: encounter_util_final_aabb !! Finalize the axis-aligned bounding box (1D) - deallocates all allocatables + procedure :: sort => encounter_check_sort_aabb_1D !! Sorts the bounding box extents along a single dimension prior to the sweep phase + procedure :: dealloc => encounter_util_dealloc_aabb !! Deallocates all allocatables + final :: encounter_util_final_aabb !! Finalize the axis-aligned bounding box (1D) - deallocates all allocatables end type type encounter_bounding_box @@ -173,6 +203,25 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical array indicating which pairs are approaching end subroutine encounter_check_sweep_aabb_single_list + module subroutine encounter_io_dump(self, param) + implicit none + class(encounter_storage(*)), intent(inout) :: self !! Encounter storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine encounter_io_dump + + module subroutine encounter_io_initialize(self, param) + implicit none + class(encounter_io_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(in) :: param + end subroutine encounter_io_initialize + + module subroutine encounter_io_write_frame(self, nc, param) + implicit none + class(encounter_snapshot), intent(in) :: self !! Swiftest encounter structure + class(encounter_io_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine encounter_io_write_frame + module subroutine encounter_setup_aabb(self, n, n_last) implicit none class(encounter_bounding_box), intent(inout) :: self !! Swiftest encounter structure @@ -219,6 +268,16 @@ module subroutine encounter_util_final_list(self) type(encounter_list), intent(inout) :: self !! Swiftest encounter list object end subroutine encounter_util_final_list + module subroutine encounter_util_final_snapshot(self) + implicit none + type(encounter_snapshot), intent(inout) :: self !! Encounter snapshot object + end subroutine encounter_util_final_snapshot + + module subroutine encounter_util_final_storage(self) + implicit none + type(encounter_storage(*)), intent(inout) :: self !! SyMBA nbody system object + end subroutine encounter_util_final_storage + module subroutine encounter_util_resize_list(self, nnew) implicit none class(encounter_list), intent(inout) :: self !! Swiftest encounter list diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index d7c949a01..e4a458a3b 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -12,7 +12,8 @@ module fraggle_classes !! !! Definition of classes and methods specific to Fraggel: The Fragment Generation Model use swiftest_globals - use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system, swiftest_cb, swiftest_pl + use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system, swiftest_cb, swiftest_pl, swiftest_storage, netcdf_parameters + use encounter_classes, only : encounter_snapshot, encounter_io_parameters implicit none public @@ -35,7 +36,8 @@ module fraggle_classes real(DP), dimension(2) :: mass !! Two-body equivalent mass of the collider bodies prior to the collision real(DP), dimension(2) :: radius !! Two-body equivalent radii of the collider bodies prior to the collision contains - procedure :: regime => fraggle_regime_colliders !! Determine which fragmentation regime the set of colliders will be + procedure :: regime => fraggle_regime_colliders !! Determine which fragmentation regime the set of colliders will be + final :: fraggle_util_final_colliders !! Finalizer will deallocate all allocatables end type fraggle_colliders !******************************************************************************************************************************** @@ -103,8 +105,29 @@ module fraggle_classes procedure :: get_ang_mtm => fraggle_util_ang_mtm !! Calcualtes the current angular momentum of the fragments procedure :: get_energy_and_momentum => fraggle_util_get_energy_momentum !! Calculates total system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) procedure :: restructure => fraggle_util_restructure !! Restructure the inputs after a failed attempt failed to find a set of positions and velocities that satisfy the energy and momentum constraints + final :: fraggle_util_final_fragments !! Finalizer will deallocate all allocatables + end type fraggle_fragments + !! NetCDF dimension and variable names for the enounter save object + type, extends(encounter_io_parameters) :: fraggle_io_encounter_parameters + contains + procedure :: initialize => fraggle_io_encounter_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object + end type fraggle_io_encounter_parameters + + !> A class that that is used to store fragmentation data between file output + type, extends(swiftest_storage) :: fraggle_storage + contains + procedure :: dump => fraggle_io_encounter_dump !! Dumps contents of encounter history to file + final :: fraggle_util_final_storage + end type fraggle_storage + + type, extends(encounter_snapshot) :: fraggle_encounter_snapshot + contains + procedure :: write_frame => fraggle_io_encounter_write_frame !! Writes a frame of encounter data to file + final :: fraggle_util_final_snapshot + end type fraggle_encounter_snapshot + interface module subroutine fraggle_generate_fragments(self, colliders, system, param, lfailure) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters @@ -116,6 +139,25 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? end subroutine fraggle_generate_fragments + module subroutine fraggle_io_encounter_dump(self, param) + implicit none + class(fraggle_storage(*)), intent(inout) :: self !! Encounter storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine fraggle_io_encounter_dump + + module subroutine fraggle_io_encounter_initialize_output(self, param) + implicit none + class(fraggle_io_encounter_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(in) :: param + end subroutine fraggle_io_encounter_initialize_output + + module subroutine fraggle_io_encounter_write_frame(self, nc, param) + implicit none + class(fraggle_encounter_snapshot), intent(in) :: self !! Swiftest encounter structure + class(encounter_io_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine fraggle_io_encounter_write_frame + module subroutine fraggle_io_log_generate(frag) implicit none class(fraggle_fragments), intent(in) :: frag @@ -239,6 +281,26 @@ module subroutine fraggle_util_construct_temporary_system(frag, system, param, t class(swiftest_parameters), allocatable, intent(out) :: tmpparam !! Output temporary configuration run parameters end subroutine fraggle_util_construct_temporary_system + module subroutine fraggle_util_final_colliders(self) + implicit none + type(fraggle_colliders), intent(inout) :: self !! Fraggle encountar storage object + end subroutine fraggle_util_final_colliders + + module subroutine fraggle_util_final_fragments(self) + implicit none + type(fraggle_fragments), intent(inout) :: self !! Fraggle encountar storage object + end subroutine fraggle_util_final_fragments + + module subroutine fraggle_util_final_storage(self) + implicit none + type(fraggle_storage(*)), intent(inout) :: self !! Fraggle encountar storage object + end subroutine fraggle_util_final_storage + + module subroutine fraggle_util_final_snapshot(self) + implicit none + type(fraggle_encounter_snapshot), intent(inout) :: self !! Fraggle encountar storage object + end subroutine fraggle_util_final_snapshot + module subroutine fraggle_util_get_energy_momentum(self, colliders, system, param, lbefore) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none diff --git a/src/modules/helio_classes.f90 b/src/modules/helio_classes.f90 index 3f12628b7..bad042664 100644 --- a/src/modules/helio_classes.f90 +++ b/src/modules/helio_classes.f90 @@ -26,7 +26,7 @@ module helio_classes contains procedure :: step => helio_step_system !! Advance the Helio nbody system forward in time by one step procedure :: initialize => helio_setup_initialize_system !! Performs Helio-specific initilization steps, including converting to DH coordinates - final :: helio_util_final_system !! Finalizes the Helio system object - deallocates all allocatables + final :: helio_util_final_system !! Finalizes the Helio system object - deallocates all allocatables end type helio_nbody_system !******************************************************************************************************************************** @@ -53,7 +53,7 @@ module helio_classes 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_util_final_pl !! Finalizes the Helio massive body object - deallocates all allocatables + final :: helio_util_final_pl !! Finalizes the Helio massive body object - deallocates all allocatables end type helio_pl !******************************************************************************************************************************** @@ -70,7 +70,7 @@ module helio_classes 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_util_final_tp !! Finalizes the Helio test particle object - deallocates all allocatables + final :: helio_util_final_tp !! Finalizes the Helio test particle object - deallocates all allocatables end type helio_tp interface diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index 5c3dc7adc..ec7dfcf16 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -36,7 +36,7 @@ module rmvs_classes !> Replace the abstract procedures with concrete ones procedure :: initialize => rmvs_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_util_final_system !! Finalizes the RMVS nbody system object - deallocates all allocatables + final :: rmvs_util_final_system !! Finalizes the RMVS nbody system object - deallocates all allocatables end type rmvs_nbody_system type, private :: rmvs_interp @@ -46,7 +46,7 @@ module rmvs_classes real(DP), dimension(:, :), allocatable :: atide !! Encountering planet's tidal acceleration value contains procedure :: dealloc => rmvs_util_dealloc_interp !! Deallocates all allocatable arrays - final :: rmvs_util_final_interp !! Finalizes the RMVS interpolated system variables object - deallocates all allocatables + final :: rmvs_util_final_interp !! Finalizes the RMVS interpolated system variables object - deallocates all allocatables end type rmvs_interp !******************************************************************************************************************************** @@ -59,7 +59,7 @@ module rmvs_classes logical :: lplanetocentric = .false. !! Flag that indicates that the object is a planetocentric set of masive bodies used for close encounter calculations contains procedure :: dealloc => rmvs_util_dealloc_cb !! Deallocates all allocatable arrays - final :: rmvs_util_final_cb !! Finalizes the RMVS central body object - deallocates all allocatables + final :: rmvs_util_final_cb !! Finalizes the RMVS central body object - deallocates all allocatables end type rmvs_cb !******************************************************************************************************************************** @@ -94,7 +94,7 @@ module rmvs_classes procedure :: sort => rmvs_util_sort_tp !! Sorts body arrays by a sortable componen procedure :: rearrange => rmvs_util_sort_rearrange_tp !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods procedure :: spill => rmvs_util_spill_tp !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) - final :: rmvs_util_final_tp !! Finalizes the RMVS test particle object - deallocates all allocatables + final :: rmvs_util_final_tp !! Finalizes the RMVS test particle object - deallocates all allocatables end type rmvs_tp !******************************************************************************************************************************** @@ -119,7 +119,7 @@ module rmvs_classes procedure :: sort => rmvs_util_sort_pl !! Sorts body arrays by a sortable componen procedure :: rearrange => rmvs_util_sort_rearrange_pl !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods procedure :: spill => rmvs_util_spill_pl !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) - final :: rmvs_util_final_pl !! Finalizes the RMVS massive body object - deallocates all allocatables + final :: rmvs_util_final_pl !! Finalizes the RMVS massive body object - deallocates all allocatables end type rmvs_pl interface diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index a3e109bb5..abda5adc2 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -39,7 +39,7 @@ module swiftest_classes character(NAMELEN) :: space_dimname = "space" !! name of the space dimension integer(I4B) :: space_dimid !! ID for the space dimension integer(I4B) :: space_varid !! ID for the space variable - character(len=1),dimension(3) :: space_coords = ["x","y","z"] !! The space dimension coordinate labels + character(len=1), dimension(3) :: space_coords = ["x","y","z"] !! The space dimension coordinate labels ! Non-dimension ids and variable names character(NAMELEN) :: ptype_varname = "particle_type" !! name of the particle type variable @@ -133,7 +133,6 @@ module swiftest_classes character(NAMELEN) :: discard_body_id_varname = "discard_body_id" !! name of the id of the other body involved in the discard end type netcdf_variables - type, extends(netcdf_variables) :: netcdf_parameters contains procedure :: close => netcdf_close !! Closes an open NetCDF file @@ -143,6 +142,27 @@ module swiftest_classes procedure :: sync => netcdf_sync !! Syncrhonize the disk and memory buffer of the NetCDF file (e.g. commit the frame files stored in memory to disk) end type netcdf_parameters + type swiftest_storage_frame + class(*), allocatable :: item + contains + procedure :: store => 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 :: util_final_storage_frame + end type + + type :: swiftest_storage(nframes) + !! An class that establishes the pattern for various storage objects + integer(I4B), len :: nframes = 4096 !! Total number of frames that can be stored + type(swiftest_storage_frame), dimension(nframes) :: frame !! Array of stored frames + integer(I4B) :: iframe = 0 !! Index of the last frame stored in the system + integer(I4B), dimension(nframes) :: tslot !! The value of the time dimension index associated with each frame + real(DP), dimension(nframes) :: tvals !! Stored time values for snapshots + contains + procedure :: dump => io_dump_storage !! Dumps storage object contents to file + procedure :: reset => util_reset_storage !! Resets a storage object by deallocating all items and resetting the frame counter to 0 + final :: util_final_storage + end type swiftest_storage + !******************************************************************************************************************************** ! swiftest_parameters class definitions !******************************************************************************************************************************** @@ -527,7 +547,6 @@ module swiftest_classes procedure :: initialize => setup_initialize_system !! Initialize the system from input files procedure :: init_particle_info => setup_initialize_particle_info_system !! Initialize the system from input files ! procedure :: step_spin => tides_step_spin_system !! Steps the spins of the massive & central bodies due to tides. - procedure :: dealloc => util_dealloc_system !! Deallocates all allocatable components of the system procedure :: set_msys => util_set_msys !! Sets the value of msys from the masses of system bodies. procedure :: get_energy_and_momentum => util_get_energy_momentum_system !! Calculates the total system energy and momentum procedure :: rescale => util_rescale_system !! Rescales the system into a new set of units @@ -535,22 +554,6 @@ module swiftest_classes generic :: write_frame => write_frame_system, write_frame_netcdf !! Generic method call for reading a frame of output data end type swiftest_nbody_system - type swiftest_storage_frame - class(*), allocatable :: item - contains - procedure :: store => util_copy_store !! Stores a snapshot of the nbody system so that later it can be retrieved for saving to file. - generic :: assignment(=) => store - end type - - type :: swiftest_storage(nframes) - !! An class that establishes the pattern for various storage objects - integer(I4B), len :: nframes = 2048 !! Total number of frames that can be stored - type(swiftest_storage_frame), dimension(nframes) :: frame !! Array of stored frames - integer(I4B) :: iframe = 0 !! The current frame number - contains - procedure :: dump => io_dump_storage !! Dumps storage object contents to file - procedure :: reset => util_reset_storage !! Resets a storage object by deallocating all items and resetting the frame counter to 0 - end type swiftest_storage abstract interface subroutine abstract_accel(self, system, param, t, lbeg) @@ -580,14 +583,6 @@ subroutine abstract_kick_body(self, system, param, t, dt, lbeg) logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. end subroutine abstract_kick_body - function abstract_read_frame(self, iu, param) result(ierr) - import DP, I4B, swiftest_base, swiftest_parameters - class(swiftest_base), intent(inout) :: self !! Swiftest base object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - integer(I4B) :: ierr !! Error code: returns 0 if the read is successful - end function abstract_read_frame - subroutine abstract_set_mu(self, cb) import swiftest_body, swiftest_cb class(swiftest_body), intent(inout) :: self !! Swiftest body object @@ -1026,10 +1021,10 @@ end subroutine netcdf_sync module function netcdf_read_frame_system(self, nc, param) result(ierr) implicit none - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(netcdf_parameters), intent(inout) :: nc !! Parameters used to for reading a NetCDF dataset to file - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - integer(I4B) :: ierr !! Error code: returns 0 if the read is successful + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(netcdf_parameters), intent(inout) :: nc !! Parameters used to for reading a NetCDF dataset to file + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + integer(I4B) :: ierr !! Error code: returns 0 if the read is successful end function netcdf_read_frame_system module subroutine netcdf_read_hdr_system(self, nc, param) @@ -1390,10 +1385,10 @@ module subroutine util_dealloc_pl(self) class(swiftest_pl), intent(inout) :: self end subroutine util_dealloc_pl - module subroutine util_dealloc_system(self) + module subroutine util_final_system(self) implicit none class(swiftest_nbody_system), intent(inout) :: self - end subroutine util_dealloc_system + end subroutine util_final_system module subroutine util_dealloc_tp(self) implicit none @@ -1472,6 +1467,17 @@ end subroutine util_fill_arr_logical end interface interface + + module subroutine util_final_storage(self) + implicit none + type(swiftest_storage(*)) :: self + end subroutine util_final_storage + + module subroutine util_final_storage_frame(self) + implicit none + type(swiftest_storage_frame) :: self + end subroutine util_final_storage_frame + pure module subroutine util_flatten_eucl_ij_to_k(n, i, j, k) !$omp declare simd(util_flatten_eucl_ij_to_k) implicit none diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 6c3aaae12..faa1e2e80 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -16,7 +16,7 @@ module symba_classes use swiftest_classes, only : swiftest_parameters, swiftest_base, swiftest_particle_info, swiftest_storage, netcdf_parameters use helio_classes, only : helio_cb, helio_pl, helio_tp, helio_nbody_system use fraggle_classes, only : fraggle_colliders, fraggle_fragments - use encounter_classes, only : encounter_list + use encounter_classes, only : encounter_list, encounter_storage, encounter_snapshot implicit none public @@ -180,38 +180,16 @@ module symba_classes end type symba_plplenc - !! NetCDF dimension and variable names for the enounter save object - type, extends(netcdf_parameters) :: symba_io_encounter_parameters - integer(I4B) :: ienc_frame = 1 !! Current frame number for the encounter history - character(STRMAX) :: enc_file = "encounter.nc" !! Encounter output file name - character(NAMELEN) :: level_varname = "level" !! Recursion depth - integer(I4B) :: level_varid !! ID for the recursion level variable - integer(I4B) :: time_dimsize = 0 !! Number of time values in snapshot - integer(I4B) :: id_dimsize = 0 !! Number of potential id values in snapshot - contains - procedure :: initialize => symba_io_encounter_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object - end type symba_io_encounter_parameters - - type, extends(swiftest_storage) :: symba_encounter_storage - !! A class that that is used to store simulation history data between file output - type(symba_io_encounter_parameters) :: nc !! NetCDF parameter object - real(DP), dimension(nframes) :: tvals !! Stored time values for snapshots - contains - procedure :: dump => symba_io_encounter_dump !! Dumps contents of encounter history to file - final :: symba_util_final_encounter_storage - end type symba_encounter_storage - - !******************************************************************************************************************************** ! symba_nbody_system class definitions and method interfaces !******************************************************************************************************************************** type, extends(helio_nbody_system) :: symba_nbody_system - class(symba_merger), allocatable :: pl_adds !! List of added bodies in mergers or collisions - class(symba_pltpenc), allocatable :: pltpenc_list !! List of massive body-test particle encounters in a single step - class(symba_plplenc), allocatable :: plplenc_list !! List of massive body-massive body encounters in a single step - class(symba_plplenc), allocatable :: plplcollision_list !! List of massive body-massive body collisions in a single step - integer(I4B) :: irec !! System recursion level - type(symba_encounter_storage(nframes=:)), allocatable :: encounter_history !! Stores encounter history for later retrieval and saving to file + class(symba_merger), allocatable :: pl_adds !! List of added bodies in mergers or collisions + class(symba_pltpenc), allocatable :: pltpenc_list !! List of massive body-test particle encounters in a single step + class(symba_plplenc), allocatable :: plplenc_list !! List of massive body-massive body encounters in a single step + class(symba_plplenc), allocatable :: plplcollision_list !! List of massive body-massive body collisions in a single step + integer(I4B) :: irec !! System recursion level + type(encounter_storage(nframes=:)), allocatable :: encounter_history !! Stores encounter history for later retrieval and saving to file contains procedure :: write_discard => symba_io_write_discard !! Write out information about discarded and merged planets and test particles in SyMBA procedure :: initialize => symba_setup_initialize_system !! Performs SyMBA-specific initilization steps @@ -220,23 +198,13 @@ module symba_classes procedure :: set_recur_levels => symba_step_set_recur_levels_system !! Sets recursion levels of bodies and encounter lists to the current 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 - procedure :: dealloc => symba_util_dealloc_system !! Deallocates all allocatable arrays - procedure :: resize_storage => symba_util_resize_storage !! Resizes the encounter history storage object so that it contains enough spaces for the number of snapshots needed procedure :: snapshot => symba_util_take_encounter_snapshot !! Take a minimal snapshot of the system through an encounter procedure :: start_encounter => symba_io_start_encounter !! Initializes the new encounter history procedure :: stop_encounter => symba_io_stop_encounter !! Saves the encounter and/or fragmentation data to file(s) - final :: symba_util_final_system !! Finalizes the SyMBA nbody system object - deallocates all allocatables + final :: symba_util_final_system !! Finalizes the SyMBA nbody system object - deallocates all allocatables 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 - contains - procedure :: write_encounter_frame => symba_io_encounter_write_frame !! Writes a frame of encounter data to file - generic :: write_frame => write_encounter_frame - final :: symba_util_final_encounter_snapshot - end type symba_encounter_snapshot - interface module function symba_collision_check_encounter(self, system, param, t, dt, irec) result(lany_collision) @@ -419,25 +387,6 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) real(DP), intent(in) :: t !! current time end subroutine symba_util_take_encounter_snapshot - module subroutine symba_io_encounter_dump(self, param) - implicit none - class(symba_encounter_storage(*)), intent(inout) :: self !! Encounter storage object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine symba_io_encounter_dump - - module subroutine symba_io_encounter_initialize_output(self, param) - implicit none - class(symba_io_encounter_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset - class(swiftest_parameters), intent(in) :: param - end subroutine symba_io_encounter_initialize_output - - module subroutine symba_io_encounter_write_frame(self, nc, param) - implicit none - class(symba_encounter_snapshot), intent(in) :: self !! Swiftest encounter structure - class(symba_io_encounter_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine symba_io_encounter_write_frame - module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) implicit none class(symba_parameters), intent(inout) :: self !! Current run configuration parameters with SyMBA additionss @@ -657,11 +606,6 @@ 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_system(self) - implicit none - class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - end subroutine symba_util_dealloc_system - module subroutine symba_util_dealloc_pl(self) implicit none class(symba_pl), intent(inout) :: self !! SyMBA massive body object @@ -711,16 +655,6 @@ module subroutine symba_util_final_encounter_list(self) type(symba_encounter), intent(inout) :: self !! SyMBA encounter list object end subroutine symba_util_final_encounter_list - module subroutine symba_util_final_encounter_snapshot(self) - implicit none - type(symba_encounter_snapshot), intent(inout) :: self !! SyMBA nbody system object - end subroutine symba_util_final_encounter_snapshot - - module subroutine symba_util_final_encounter_storage(self) - implicit none - type(symba_encounter_storage(*)), intent(inout) :: self !! SyMBA nbody system object - end subroutine symba_util_final_encounter_storage - module subroutine symba_util_final_kin(self) implicit none type(symba_kinship), intent(inout) :: self !! SyMBA kinship object @@ -789,12 +723,6 @@ module subroutine symba_util_resize_pl(self, nnew) integer(I4B), intent(in) :: nnew !! New size neded end subroutine symba_util_resize_pl - module subroutine symba_util_resize_storage(self, nnew) - implicit none - class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - integer(I4B), intent(in) :: nnew !! New size of list needed - end subroutine symba_util_resize_storage - module subroutine symba_util_resize_tp(self, nnew) implicit none class(symba_tp), intent(inout) :: self !! SyMBA massive body object diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90 index 34922c29a..100e606c4 100644 --- a/src/modules/whm_classes.f90 +++ b/src/modules/whm_classes.f90 @@ -58,7 +58,7 @@ module whm_classes procedure :: spill => whm_util_spill_pl !!"Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) procedure :: setup => whm_setup_pl !! Constructor method - Allocates space for the input number of bodiess procedure :: step => whm_step_pl !! Steps the body forward one stepsize - final :: whm_util_final_pl !! Finalizes the WHM massive body object - deallocates all allocatables + final :: whm_util_final_pl !! Finalizes the WHM massive body object - deallocates all allocatables end type whm_pl !******************************************************************************************************************************** @@ -75,7 +75,7 @@ module whm_classes procedure :: accel => whm_kick_getacch_tp !! Compute heliocentric accelerations of test particles procedure :: kick => whm_kick_vh_tp !! Kick heliocentric velocities of test particles procedure :: step => whm_step_tp !! Steps the particle forward one stepsize - final :: whm_util_final_tp !! Finalizes the WHM test particle object - deallocates all allocatables + final :: whm_util_final_tp !! Finalizes the WHM test particle object - deallocates all allocatables end type whm_tp !******************************************************************************************************************************** @@ -87,7 +87,7 @@ module whm_classes !> Replace the abstract procedures with concrete ones procedure :: initialize => whm_setup_initialize_system !! Performs WHM-specific initilization steps, like calculating the Jacobi masses procedure :: step => whm_step_system !! Advance the WHM nbody system forward in time by one step - final :: whm_util_final_system !! Finalizes the WHM system object - deallocates all allocatables + final :: whm_util_final_system !! Finalizes the WHM system object - deallocates all allocatables end type whm_nbody_system interface diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index 31de061de..ad5d7bbeb 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -77,7 +77,7 @@ module function netcdf_get_old_t_final_system(self, param) result(old_t_final) integer(I4B) :: itmax, idmax real(DP), dimension(:), allocatable :: vals real(DP), dimension(1) :: rtemp - real(DP), dimension(NDIM) :: vectemp, rot0, Ip0, Lnow + real(DP), dimension(NDIM) :: rot0, Ip0, Lnow real(DP) :: KE_orb_orig, KE_spin_orig, PE_orig call param%nc%open(param) @@ -438,7 +438,7 @@ module function netcdf_read_frame_system(self, nc, param) result(ierr) implicit none ! Arguments class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset + class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Return integer(I4B) :: ierr !! Error code: returns 0 if the read is successful @@ -449,6 +449,8 @@ module function netcdf_read_frame_system(self, nc, param) result(ierr) integer(I4B), dimension(:), allocatable :: itemp logical, dimension(:), allocatable :: validmask, tpmask, plmask + tslot = param%ioutput + call nc%open(param, readonly=.true.) call self%read_hdr(nc, param) @@ -457,8 +459,6 @@ module function netcdf_read_frame_system(self, nc, param) result(ierr) call pl%setup(npl, param) call tp%setup(ntp, param) - tslot = param%ioutput - call check( nf90_inquire_dimension(nc%id, nc%id_dimid, len=idmax), "netcdf_read_frame_system nf90_inquire_dimension id_dimid" ) allocate(rtemp(idmax)) allocate(vectemp(NDIM,idmax)) @@ -1021,19 +1021,19 @@ module subroutine netcdf_write_frame_base(self, nc, param) !! Note: If outputting to orbital elements, but sure that the conversion is done prior to calling this method implicit none ! Arguments - class(swiftest_base), intent(in) :: self !! Swiftest particle object - class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_base), intent(in) :: self !! Swiftest particle object + class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i, j, tslot, idslot, old_mode integer(I4B), dimension(:), allocatable :: ind real(DP), dimension(NDIM) :: vh !! Temporary variable to store heliocentric velocity values when converting from pseudovelocity in GR-enabled runs real(DP) :: a, e, inc, omega, capom, capm, varpi, lam, f, cape, capf - call self%write_info(nc, param) - tslot = param%ioutput + call self%write_info(nc, param) + call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "netcdf_write_frame_base nf90_set_fill" ) select type(self) class is (swiftest_body) diff --git a/src/rmvs/rmvs_util.f90 b/src/rmvs/rmvs_util.f90 index 0fc1ed272..eb31db53d 100644 --- a/src/rmvs/rmvs_util.f90 +++ b/src/rmvs/rmvs_util.f90 @@ -232,7 +232,8 @@ module subroutine rmvs_util_final_system(self) ! Arguments type(rmvs_nbody_system), intent(inout) :: self !! RMVS nbody system object - call self%dealloc() + if (allocated(self%vbeg)) deallocate(self%vbeg) + call whm_util_final_system(self%whm_nbody_system) return end subroutine rmvs_util_final_system diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index ea6d3db23..e95505b9b 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -68,7 +68,6 @@ module subroutine setup_construct_system(system, param) allocate(symba_pltpenc :: system%pltpenc_list) allocate(symba_plplenc :: system%plplenc_list) allocate(symba_plplenc :: system%plplcollision_list) - allocate(symba_encounter_storage :: system%encounter_history) end select case (RINGMOONS) write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled' diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index bc0f1b985..e790e1685 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -11,194 +11,6 @@ use swiftest contains - module subroutine symba_io_encounter_dump(self, param) - !! author: David A. Minton - !! - !! Dumps the time history of an encounter to file. - implicit none - ! Arguments - class(symba_encounter_storage(*)), intent(inout) :: self !! Encounter storage object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: i - - ! Most of this is just temporary test code just to get something working. Eventually this should get cleaned up. - - do i = 1, self%nframes - if (allocated(self%frame(i)%item)) then - select type(snapshot => self%frame(i)%item) - class is (symba_encounter_snapshot) - self%nc%ienc_frame = i - call snapshot%write_frame(self%nc,param) - end select - else - exit - end if - end do - - - return - end subroutine symba_io_encounter_dump - - - module subroutine symba_io_encounter_initialize_output(self, param) - !! author: David A. Minton - !! - !! Initialize a NetCDF encounter file system. This is a simplified version of the main simulation output NetCDF file, but with fewer variables. - use, intrinsic :: ieee_arithmetic - use netcdf - implicit none - ! Arguments - class(symba_io_encounter_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: nvar, varid, vartype - real(DP) :: dfill - real(SP) :: sfill - logical :: fileExists - character(len=STRMAX) :: errmsg - integer(I4B) :: ndims - - associate(nc => self) - dfill = ieee_value(dfill, IEEE_QUIET_NAN) - sfill = ieee_value(sfill, IEEE_QUIET_NAN) - - select case (param%out_type) - case("NETCDF_FLOAT") - self%out_type = NF90_FLOAT - case("NETCDF_DOUBLE") - self%out_type = NF90_DOUBLE - end select - - ! Check if the file exists, and if it does, delete it - inquire(file=nc%enc_file, exist=fileExists) - if (fileExists) then - open(unit=LUN, file=nc%enc_file, status="old", err=667, iomsg=errmsg) - close(unit=LUN, status="delete") - end if - - call check( nf90_create(nc%enc_file, NF90_NETCDF4, nc%id), "symba_io_encounter_initialize_output nf90_create" ) - - ! Dimensions - call check( nf90_def_dim(nc%id, nc%time_dimname, nc%time_dimsize, nc%time_dimid), "symba_io_encounter_initialize_output nf90_def_dim time_dimid" ) ! Simulation time dimension - call check( nf90_def_dim(nc%id, nc%space_dimname, NDIM , nc%space_dimid), "symba_io_encounter_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension - call check( nf90_def_dim(nc%id, nc%id_dimname, param%maxid, nc%id_dimid), "symba_io_encounter_initialize_output nf90_def_dim id_dimid" ) ! dimension to store particle id numbers - call check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "symba_io_encounter_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) - - ! Dimension coordinates - call check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, nc%time_dimid, nc%time_varid), "symba_io_encounter_initialize_output nf90_def_var time_varid" ) - call check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "symba_io_encounter_initialize_output nf90_def_var space_varid" ) - call check( nf90_def_var(nc%id, nc%id_dimname, NF90_INT, nc%id_dimid, nc%id_varid), "symba_io_encounter_initialize_output nf90_def_var id_varid" ) - - ! Variables - call check( nf90_def_var(nc%id, nc%name_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], nc%name_varid), "symba_io_encounter_initialize_output nf90_def_var name_varid" ) - 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" ) - 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 - if (param%lrotation) then - call check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%Ip_varid), "symba_io_encounter_initialize_output nf90_def_var Ip_varid" ) - call check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%rot_varid), "symba_io_encounter_initialize_output nf90_def_var rot_varid" ) - end if - - call check( nf90_inquire(nc%id, nVariables=nvar), "symba_io_encounter_initialize_output nf90_inquire nVariables" ) - do varid = 1, nvar - call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "symba_io_encounter_initialize_output nf90_inquire_variable" ) - select case(vartype) - case(NF90_INT) - call check( nf90_def_var_fill(nc%id, varid, 0, NF90_FILL_INT), "symba_io_encounter_initialize_output nf90_def_var_fill NF90_INT" ) - case(NF90_FLOAT) - call check( nf90_def_var_fill(nc%id, varid, 0, sfill), "symba_io_encounter_initialize_output nf90_def_var_fill NF90_FLOAT" ) - case(NF90_DOUBLE) - call check( nf90_def_var_fill(nc%id, varid, 0, dfill), "symba_io_encounter_initialize_output nf90_def_var_fill NF90_DOUBLE" ) - case(NF90_CHAR) - call check( nf90_def_var_fill(nc%id, varid, 0, 0), "symba_io_encounter_initialize_output nf90_def_var_fill NF90_CHAR" ) - end select - end do - - ! Take the file out of define mode - call check( nf90_enddef(nc%id), "symba_io_encounter_initialize_output nf90_enddef" ) - - ! Add in the space dimension coordinates - call check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "symba_io_encounter_initialize_output nf90_put_var space" ) - end associate - - return - - 667 continue - write(*,*) "Error creating encounter output file. " // trim(adjustl(errmsg)) - call util_exit(FAILURE) - end subroutine symba_io_encounter_initialize_output - - - module subroutine symba_io_encounter_write_frame(self, nc, param) - !! author: David A. Minton - !! - !! Write a frame of output of an encounter trajectory. - use netcdf - implicit none - ! Arguments - class(symba_encounter_snapshot), intent(in) :: self !! Swiftest encounter structure - class(symba_io_encounter_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: i, tslot, idslot, old_mode, npl, ntp - character(len=NAMELEN) :: charstring - - call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "symba_io_encounter_write_frame nf90_set_fill" ) - - 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) - npl = pl%nbody - do i = 1, npl - idslot = pl%id(i) - call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[idslot]), "symba_io_encounter_write_frame nf90_put_var pl id_varid" ) - 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" ) - - 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" ) - - if (param%lrotation) then - call check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var pl Ip_varid" ) - call check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1,idslot, tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var pl rotx_varid" ) - end if - - charstring = trim(adjustl(pl%info(i)%name)) - call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "symba_io_encounter_write_frame nf90_put_var pl name_varid" ) - 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) - call check( nf90_put_var(nc%id, nc%id_varid, tp%id(i), start=[idslot]), "symba_io_encounter_write_frame nf90_put_var tp id_varid" ) - call check( nf90_put_var(nc%id, nc%rh_varid, tp%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var tp rh_varid" ) - call check( nf90_put_var(nc%id, nc%vh_varid, tp%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var tp vh_varid" ) - - charstring = trim(adjustl(tp%info(i)%name)) - call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "symba_io_encounter_write_frame nf90_put_var tp name_varid" ) - charstring = trim(adjustl(tp%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 tp particle_type_varid" ) - end do - end associate - - call check( nf90_set_fill(nc%id, old_mode, old_mode) ) - - return - end subroutine symba_io_encounter_write_frame - - module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! @@ -386,12 +198,11 @@ module subroutine symba_io_start_encounter(self, param, t) class(symba_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time - if (.not. allocated(self%encounter_history)) allocate(symba_encounter_storage :: self%encounter_history) + if (.not. allocated(self%encounter_history)) then + allocate(encounter_storage :: self%encounter_history) + end if call self%encounter_history%reset() - ! Empty out the time slot array for the next pass - self%encounter_history%tvals(:) = huge(1.0_DP) - ! Take the snapshot at the start of the encounter call self%snapshot(param, t) @@ -413,13 +224,7 @@ module subroutine symba_io_stop_encounter(self, param, t) ! Create and save the output file for this encounter - ! Figure out how many time slots we need - do i = 1, self%encounter_history%nframes - if (self%t + param%dt <= self%encounter_history%tvals(i)) then - self%encounter_history%nc%time_dimsize = i - exit - end if - end do + self%encounter_history%nc%time_dimsize = maxval(self%encounter_history%tslot(:)) write(self%encounter_history%nc%enc_file, '("encounter_",I0.6,".nc")') param%iloop diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index b96270fc9..d53bd442b 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -238,25 +238,6 @@ module subroutine symba_util_dealloc_merger(self) end subroutine symba_util_dealloc_merger - module subroutine symba_util_dealloc_system(self) - !! author: David A. Minton - !! - !! Deallocates all allocatabale arrays - implicit none - ! Arguments - class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - - if (allocated(self%pl_adds)) deallocate(self%pl_adds) - if (allocated(self%pltpenc_list)) deallocate(self%pltpenc_list) - if (allocated(self%plplenc_list)) deallocate(self%plplenc_list) - if (allocated(self%plplcollision_list)) deallocate(self%plplcollision_list) - - call util_dealloc_system(self) - - return - end subroutine symba_util_dealloc_system - - module subroutine symba_util_dealloc_pl(self) !! author: David A. Minton !! @@ -471,6 +452,7 @@ module subroutine symba_util_final_merger(self) return end subroutine symba_util_final_merger + module subroutine symba_util_final_pl(self) !! author: David A. Minton !! @@ -493,34 +475,15 @@ module subroutine symba_util_final_system(self) ! Argument type(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - call self%dealloc() - - return - end subroutine symba_util_final_system - - module subroutine symba_util_final_encounter_snapshot(self) - !! author: David A. Minton - !! - !! Finalize the SyMBA encounter system snapshot object - deallocates all allocatables - implicit none - type(symba_encounter_snapshot), intent(inout) :: self !! SyMBA nbody system object + if (allocated(self%pl_adds)) deallocate(self%pl_adds) + if (allocated(self%pltpenc_list)) deallocate(self%pltpenc_list) + if (allocated(self%plplenc_list)) deallocate(self%plplenc_list) + if (allocated(self%plplcollision_list)) deallocate(self%plplcollision_list) - call self%dealloc() + call helio_util_final_system(self%helio_nbody_system) return - end subroutine symba_util_final_encounter_snapshot - - - module subroutine symba_util_final_encounter_storage(self) - !! author: David A. Minton - !! - !! Finalize the SyMBA nbody system object - deallocates all allocatables - implicit none - ! Argument - type(symba_encounter_storage(*)), intent(inout) :: self !! SyMBA nbody system object - - return - end subroutine symba_util_final_encounter_storage + end subroutine symba_util_final_system module subroutine symba_util_final_tp(self) @@ -906,7 +869,7 @@ module subroutine symba_util_resize_pl(self, nnew) end subroutine symba_util_resize_pl - module subroutine symba_util_resize_storage(self, nnew) + subroutine symba_util_save_storage(system, snapshot, t) !! 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. @@ -915,40 +878,53 @@ module subroutine symba_util_resize_storage(self, nnew) !! Memory usage grows by a factor of 2 each time it fills up, but no more. implicit none ! Arguments - class(symba_nbody_system), intent(inout) :: self !! Swiftest encounter list - integer(I4B), intent(in) :: nnew !! New size of list needed + type(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(encounter_snapshot), intent(in) :: snapshot !! Encounter snapshot object + real(DP), intent(in) :: t !! The time of the snapshot ! Internals - type(symba_encounter_storage(nframes=:)), allocatable :: tmp - integer(I4B) :: i, nold, nbig, iframe_old = 0 - logical :: lmalloc + type(encounter_storage(nframes=:)), allocatable :: tmp + integer(I4B) :: i, nnew, nold, nbig + ! Advance the snapshot frame counter + system%encounter_history%iframe = system%encounter_history%iframe + 1 - lmalloc = allocated(self%encounter_history) - if (lmalloc) then - nold = self%encounter_history%nframes - iframe_old = self%encounter_history%iframe - else - nold = 0 - end if + ! Check to make sure the current encounter_history object is big enough. If not, grow it by a factor of 2 + nnew = system%encounter_history%iframe + nold = system%encounter_history%nframes if (nnew > nold) then nbig = nold do while (nbig < nnew) nbig = nbig * 2 end do - allocate(symba_encounter_storage(nbig) :: tmp) - 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) - end do - deallocate(self%encounter_history) - end if - call move_alloc(tmp,self%encounter_history) - self%encounter_history%iframe = iframe_old + allocate(encounter_storage(nbig) :: tmp) + tmp%tvals(1:nold) = system%encounter_history%tvals(1:nold) + tmp%tvals(nold+1:nbig) = huge(1.0_DP) + tmp%tslot(1:nold) = system%encounter_history%tslot(1:nold) + tmp%tslot(nold+1:nbig) = 0 + tmp%iframe = system%encounter_history%iframe + + do i = 1, nold + if (allocated(system%encounter_history%frame(i)%item)) call move_alloc(system%encounter_history%frame(i)%item, tmp%frame(i)%item) + end do + deallocate(system%encounter_history) + call move_alloc(tmp,system%encounter_history) + nnew = nbig 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 + do i = 1, nnew + if (t <= system%encounter_history%tvals(i)) then + system%encounter_history%tvals(i) = t + system%encounter_history%tslot(nnew) = i + system%encounter_history%frame(nnew) = snapshot + exit + end if + end do + return - end subroutine symba_util_resize_storage + end subroutine symba_util_save_storage module subroutine symba_util_resize_tp(self, nnew) @@ -1317,19 +1293,16 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) implicit none ! Internals class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! current time + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! current time ! Arguments - type(symba_encounter_snapshot), allocatable :: snapshot + type(encounter_snapshot) :: snapshot integer(I4B) :: i, npl_snap, ntp_snap associate(npl => self%pl%nbody, ntp => self%tp%nbody) - 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 @@ -1338,79 +1311,79 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) class is (symba_pl) select type (tp => self%tp) class is (symba_tp) - if (npl > 0) then - pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE .and. pl%levelg(1:npl) == self%irec - npl_snap = count(pl%lmask(1:npl)) - end if - if (ntp > 0) then - tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE .and. tp%levelg(1:ntp) == self%irec - ntp_snap = count(tp%lmask(1:ntp)) - end if - snapshot%pl%nbody = npl_snap - - ! Take snapshot of the currently encountering massive bodies - if (npl_snap > 0) then - allocate(snapshot%pl%id(npl_snap)) - allocate(snapshot%pl%info(npl_snap)) - allocate(snapshot%pl%Gmass(npl_snap)) - 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)) - allocate(snapshot%pl%rh(NDIM,npl_snap)) - allocate(snapshot%pl%vh(NDIM,npl_snap)) - do i = 1, NDIM - snapshot%pl%rh(i,:) = pack(pl%rh(i,1:npl), pl%lmask(1:npl)) - snapshot%pl%vh(i,:) = pack(pl%vb(i,1:npl), pl%lmask(1:npl)) - end do - if (param%lclose) then - allocate(snapshot%pl%radius(npl_snap)) - snapshot%pl%radius(:) = pack(pl%radius(1:npl), pl%lmask(1:npl)) - end if + allocate(symba_pl :: snapshot%pl) + allocate(symba_tp :: snapshot%tp) - if (param%lrotation) then - allocate(snapshot%pl%Ip(NDIM,npl_snap)) - allocate(snapshot%pl%rot(NDIM,npl_snap)) - do i = 1, NDIM - snapshot%pl%Ip(i,:) = pack(pl%Ip(i,1:npl), pl%lmask(1:npl)) - snapshot%pl%rot(i,:) = pack(pl%rot(i,1:npl), pl%lmask(1:npl)) - end do - end if - call snapshot%pl%sort("id", ascending=.true.) - end if + select type(pl_snap => snapshot%pl) + class is (symba_pl) + select type(tp_snap => snapshot%tp) + class is (symba_tp) - ! Take snapshot of the currently encountering test particles - snapshot%tp%nbody = ntp_snap - if (ntp_snap > 0) then - allocate(snapshot%tp%id(ntp_snap)) - allocate(snapshot%tp%info(ntp_snap)) - snapshot%tp%id(:) = pack(tp%id(1:ntp), tp%lmask(1:ntp)) - snapshot%tp%info(:) = pack(tp%info(1:ntp), tp%lmask(1:ntp)) - allocate(snapshot%tp%rh(NDIM,ntp_snap)) - allocate(snapshot%tp%vh(NDIM,ntp_snap)) - do i = 1, NDIM - snapshot%tp%rh(i,:) = pack(tp%rh(i,1:ntp), tp%lmask(1:ntp)) - snapshot%tp%vh(i,:) = pack(tp%vh(i,1:ntp), tp%lmask(1:ntp)) - end do - end if + if (npl > 0) then + pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE .and. pl%levelg(1:npl) == self%irec + npl_snap = count(pl%lmask(1:npl)) + end if + if (ntp > 0) then + tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE .and. tp%levelg(1:ntp) == self%irec + ntp_snap = count(tp%lmask(1:ntp)) + end if + pl_snap%nbody = npl_snap + + ! Take snapshot of the currently encountering massive bodies + if (npl_snap > 0) then + allocate(pl_snap%id(npl_snap)) + allocate(pl_snap%info(npl_snap)) + allocate(pl_snap%Gmass(npl_snap)) + + allocate(pl_snap%levelg(npl_snap)) + pl_snap%levelg(:) = pack(pl%levelg(1:npl), pl%lmask(1:npl)) + pl_snap%id(:) = pack(pl%id(1:npl), pl%lmask(1:npl)) + pl_snap%info(:) = pack(pl%info(1:npl), pl%lmask(1:npl)) + pl_snap%Gmass(:) = pack(pl%Gmass(1:npl), pl%lmask(1:npl)) + allocate(pl_snap%rh(NDIM,npl_snap)) + allocate(pl_snap%vh(NDIM,npl_snap)) + do i = 1, NDIM + pl_snap%rh(i,:) = pack(pl%rh(i,1:npl), pl%lmask(1:npl)) + pl_snap%vh(i,:) = pack(pl%vb(i,1:npl), pl%lmask(1:npl)) + end do + if (param%lclose) then + allocate(pl_snap%radius(npl_snap)) + pl_snap%radius(:) = pack(pl%radius(1:npl), pl%lmask(1:npl)) + end if + + if (param%lrotation) then + allocate(pl_snap%Ip(NDIM,npl_snap)) + allocate(pl_snap%rot(NDIM,npl_snap)) + do i = 1, NDIM + pl_snap%Ip(i,:) = pack(pl%Ip(i,1:npl), pl%lmask(1:npl)) + pl_snap%rot(i,:) = pack(pl%rot(i,1:npl), pl%lmask(1:npl)) + end do + end if + call pl_snap%sort("id", ascending=.true.) + end if + + ! Take snapshot of the currently encountering test particles + tp_snap%nbody = ntp_snap + if (ntp_snap > 0) then + allocate(tp_snap%id(ntp_snap)) + allocate(tp_snap%info(ntp_snap)) + tp_snap%id(:) = pack(tp%id(1:ntp), tp%lmask(1:ntp)) + tp_snap%info(:) = pack(tp%info(1:ntp), tp%lmask(1:ntp)) + allocate(tp_snap%rh(NDIM,ntp_snap)) + allocate(tp_snap%vh(NDIM,ntp_snap)) + do i = 1, NDIM + tp_snap%rh(i,:) = pack(tp%rh(i,1:ntp), tp%lmask(1:ntp)) + tp_snap%vh(i,:) = pack(tp%vh(i,1:ntp), tp%lmask(1:ntp)) + end do + end if + end select + end select - ! 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 - do i = 1, self%encounter_history%nframes - if (t <= self%encounter_history%tvals(i)) then - snapshot%tslot = i - self%encounter_history%tvals(i) = t - exit - end if - end do - - self%encounter_history%frame(self%encounter_history%iframe) = snapshot + ! Save the snapshot + call symba_util_save_storage(self,snapshot,t) end select - end select + end select end associate return diff --git a/src/util/util_dealloc.f90 b/src/util/util_dealloc.f90 index b3fb38fd9..54151f567 100644 --- a/src/util/util_dealloc.f90 +++ b/src/util/util_dealloc.f90 @@ -72,25 +72,6 @@ module subroutine util_dealloc_pl(self) return end subroutine util_dealloc_pl - - module subroutine util_dealloc_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 util_dealloc_system - - module subroutine util_dealloc_tp(self) !! author: David A. Minton !! diff --git a/src/util/util_final.f90 b/src/util/util_final.f90 new file mode 100644 index 000000000..4f4a6dd28 --- /dev/null +++ b/src/util/util_final.f90 @@ -0,0 +1,62 @@ +!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh +!! This file is part of Swiftest. +!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +!! You should have received a copy of the GNU General Public License along with Swiftest. +!! If not, see: https://www.gnu.org/licenses. + +submodule (swiftest_classes) s_util_final + use swiftest +contains + + module subroutine util_final_storage(self) + !! author: David A. Minton + !! + !! 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 + + return + end subroutine util_final_storage + + module subroutine util_final_storage_frame(self) + !! author: David A. Minton + !! + !! Finalizer for the storage frame data type + implicit none + type(swiftest_storage_frame) :: self + + if (allocated(self%item)) deallocate(self%item) + + return + end subroutine util_final_storage_frame + + + module subroutine util_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 util_final_system + + +end submodule s_util_final diff --git a/src/util/util_reset.f90 b/src/util/util_reset.f90 index 569846a68..7bb8d5ee3 100644 --- a/src/util/util_reset.f90 +++ b/src/util/util_reset.f90 @@ -24,6 +24,8 @@ module subroutine util_reset_storage(self) do i = 1, self%nframes if (allocated(self%frame(i)%item)) deallocate(self%frame(i)%item) end do + self%tslot(:) = 0 + self%tvals(:) = huge(1.0_DP) self%iframe = 0 return diff --git a/src/whm/whm_util.f90 b/src/whm/whm_util.f90 index 9c6efdd41..c58f5730f 100644 --- a/src/whm/whm_util.f90 +++ b/src/whm/whm_util.f90 @@ -117,7 +117,7 @@ module subroutine whm_util_final_system(self) ! Arguments type(whm_nbody_system), intent(inout) :: self !! WHM nbody system object - call self%dealloc() + call util_final_system(self) return end subroutine whm_util_final_system