From d06b7da85f3b9b32a1c2de6067b4f68ad5ed8c98 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 6 Dec 2022 17:56:58 -0500 Subject: [PATCH 01/13] More updates to movie. Now interpolate across missing data in recursion step data --- examples/Fragmentation/Fragmentation_Movie.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 80f75a86a..0bff9416a 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -83,9 +83,14 @@ class AnimatedScatter(object): """An animated scatter plot using matplotlib.animations.FuncAnimation.""" - def __init__(self, sim, animfile, title, nskip=1): - tgood = sim.enc.where(sim.enc['Gmass'] > 9e-8).time - self.ds = sim.enc.sel(time=tgood) + 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") + nframes = int(self.ds['time'].size) self.sim = sim self.title = title @@ -134,13 +139,12 @@ def center(Gmass, x, y): x = x[~np.isnan(x)] y = y[~np.isnan(y)] Gmass = Gmass[~np.isnan(Gmass)] - x = x[Gmass] x_com = np.sum(Gmass * x) / np.sum(Gmass) - y_com = #np.sum(Gmass * y) / np.sum(Gmass) + y_com = np.sum(Gmass * y) / np.sum(Gmass) return x_com, y_com Gmass, rh, point_rad = next(self.data_stream(frame)) - #x_com, y_com = center(Gmass, rh[:,0], rh[:,1]) + x_com, y_com = center(Gmass, rh[:,0], rh[:,1]) self.scatter_artist.set_offsets(np.c_[rh[:,0] - x_com, rh[:,1] - y_com]) self.scatter_artist.set_sizes(point_rad**2) return self.scatter_artist, @@ -183,4 +187,4 @@ def data_stream(self, frame=0): sim.set_parameter(fragmentation=True, fragmentation_save="TRAJECTORY", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) sim.run(dt=1e-4, tstop=2.0e-3, istep_out=1, dump_cadence=0) - anim = AnimatedScatter(sim,movie_filename,movie_titles[style],nskip=1) \ No newline at end of file + anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1) \ No newline at end of file From 3145e098043ce9dfc14d0617450d9612b07df758 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 7 Dec 2022 08:28:12 -0500 Subject: [PATCH 02/13] Fixed issues related to joining all the encounter data files --- python/swiftest/swiftest/io.py | 3 ++- python/swiftest/swiftest/simulation_class.py | 14 ++++++++------ 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index f03fb622f..a4f370773 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -817,11 +817,12 @@ def process_netcdf_input(ds, param): ds : xarray dataset """ + ds = ds.where(~np.isnan(ds.id) ,drop=True) if param['OUT_TYPE'] == "NETCDF_DOUBLE": ds = fix_types(ds,ftype=np.float64) elif param['OUT_TYPE'] == "NETCDF_FLOAT": ds = fix_types(ds,ftype=np.float32) - ds = ds.where(ds.id >=0 ,drop=True) + # Check if the name variable contains unique values. If so, make name the dimension instead of id if "id" in ds.dims: if len(np.unique(ds['name'])) == len(ds['name']): diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 4ac07386c..9e18a8ff8 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -323,12 +323,14 @@ def __init__(self,read_param: bool = False, read_old_output_file: bool = False, self.simdir = Path(simdir) if self.simdir.exists(): if not self.simdir.is_dir(): - msg = f"Cannot create the {self.simdir} directory: File exists." + msg = f"Cannot create the {self.simdir.resolve()} directory: File exists." msg += "\nDelete the file or change the location of param_file" - warnings.warn(msg,stacklevel=2) + raise NotADirectoryError(msg) else: - self.simdir.mkdir(parents=True, exist_ok=False) - + if read_old_output_file or read_param: + raise NotADirectoryError(f"Cannot find directory {self.simdir.resolve()} ") + else: + self.simdir.mkdir(parents=True, exist_ok=False) # Set the location of the parameter input file, choosing the default if it isn't specified. param_file = kwargs.pop("param_file",Path.cwd() / self.simdir / "param.in") @@ -376,7 +378,7 @@ def __init__(self,read_param: bool = False, read_old_output_file: bool = False, if os.path.exists(binpath): self.read_output_file() else: - warnings.warn(f"BIN_OUT file {binpath} not found.",stacklevel=2) + raise FileNotFoundError(f"BIN_OUT file {binpath} not found.") return def _run_swiftest_driver(self): @@ -2770,7 +2772,7 @@ def _preprocess(ds, param): return io.process_netcdf_input(ds,param) partial_func = partial(_preprocess, param=self.param) - self.enc = xr.open_mfdataset(enc_files,parallel=True,preprocess=partial_func,mask_and_scale=True) + self.enc = xr.open_mfdataset(enc_files,parallel=True,combine="nested",concat_dim="time",join="left",preprocess=partial_func,mask_and_scale=True) self.enc = io.process_netcdf_input(self.enc, self.param) return From f4fec1e02d447399afe8c2ca560ef3315bac0b3a Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 7 Dec 2022 09:06:48 -0500 Subject: [PATCH 03/13] Fixed more bugs in encounter data processing --- examples/.gitignore | 8 ++-- examples/Fragmentation/Fragmentation_Movie.py | 6 +-- python/swiftest/swiftest/simulation_class.py | 41 +++++++++++-------- 3 files changed, 30 insertions(+), 25 deletions(-) diff --git a/examples/.gitignore b/examples/.gitignore index 7c5c72692..ad990dfcb 100644 --- a/examples/.gitignore +++ b/examples/.gitignore @@ -1,6 +1,6 @@ * !.gitignore -!Basic_Simulation/* -!Fragmentation/* -!helio_gr_test/* -!whm_gr_test/* \ No newline at end of file +!Basic_Simulation +!Fragmentation +!helio_gr_test +!whm_gr_test \ No newline at end of file diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 0bff9416a..fa1203a8e 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -40,8 +40,8 @@ movie_titles = dict(zip(available_movie_styles, movie_title_list)) # These initial conditions were generated by trial and error -pos_vectors = {"disruption_headon" : [np.array([1.0, -2.807993e-05, 0.0]), - np.array([1.0, 2.807993e-05 ,0.0])], +pos_vectors = {"disruption_headon" : [np.array([1.0, -5.0e-05, 0.0]), + np.array([1.0, 5.0e-05 ,0.0])], "supercatastrophic_off_axis": [np.array([1.0, -4.2e-05, 0.0]), np.array([1.0, 4.2e-05, 0.0])], "hitandrun" : [np.array([1.0, -2.0e-05, 0.0]), @@ -185,6 +185,6 @@ 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=2.0e-3, istep_out=1, dump_cadence=0) + sim.run(dt=1e-4, tstop=4.0e-3, istep_out=1, dump_cadence=0) anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1) \ No newline at end of file diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 9e18a8ff8..9b9bd02dc 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -17,6 +17,7 @@ from swiftest import __file__ as _pyfile import json import os +from glob import glob from pathlib import Path import datetime import xarray as xr @@ -2667,23 +2668,23 @@ def convert(self, param_file, newcodename="Swiftest", plname="pl.swiftest.in", t Parameters ---------- - param_file : string - File name of the input parameter file - newcodename : string - Name of the desired format (Swift/Swifter/Swiftest) - plname : string - File name of the massive body input file - tpname : string - File name of the test particle input file - cbname : string - File name of the central body input file - conversion_questions : dictronary - Dictionary of additional parameters required to convert between formats + param_file : string + File name of the input parameter file + newcodename : string + Name of the desired format (Swift/Swifter/Swiftest) + plname : string + File name of the massive body input file + tpname : string + File name of the test particle input file + cbname : string + File name of the central body input file + conversion_questions : dictronary + Dictionary of additional parameters required to convert between formats Returns ------- - oldparam : xarray dataset - The old parameter configuration. + oldparam : xarray dataset + The old parameter configuration. """ oldparam = self.param if self.codename == newcodename: @@ -2720,12 +2721,12 @@ def read_output_file(self,read_init_cond : bool = True): Parameters ---------- - read_init_cond : bool - Read in an initial conditions file along with the output file. Default is True + read_init_cond : bool + Read in an initial conditions file along with the output file. Default is True Returns ------- - self.data : xarray dataset + self.data : xarray dataset """ # Make a temporary copy of the parameter dictionary so we can supply the absolute path of the binary file @@ -2765,7 +2766,8 @@ def read_output_file(self,read_init_cond : bool = True): def read_encounter(self): if self.verbose: print("Reading encounter history file as .enc") - enc_files = self.simdir.glob("**/encounter_*.nc") + enc_files = glob(f"{self.simdir}{os.path.sep}encounter_*.nc") + enc_files.sort() # This is needed in order to pass the param argument down to the io.process_netcdf_input function def _preprocess(ds, param): @@ -2774,6 +2776,9 @@ def _preprocess(ds, param): self.enc = xr.open_mfdataset(enc_files,parallel=True,combine="nested",concat_dim="time",join="left",preprocess=partial_func,mask_and_scale=True) self.enc = io.process_netcdf_input(self.enc, self.param) + # Remove any overlapping time values + tgood,tid = np.unique(self.enc.time,return_index=True) + self.enc = self.enc.isel(time=tid) return From d09d8c5e90b68a85e041fe72334a0292736ca076 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 7 Dec 2022 11:46:09 -0500 Subject: [PATCH 04/13] Made a number of improvements to the handling of encounters. time and id dimensions are now limited, which significantly reduces the file sizes of the individual encounter data files --- src/modules/symba_classes.f90 | 12 +++--- src/symba/symba_io.f90 | 79 ++++++++++++++++++++++++----------- src/symba/symba_util.f90 | 11 +++-- 3 files changed, 66 insertions(+), 36 deletions(-) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 68b94c373..6c3aaae12 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -182,12 +182,12 @@ module symba_classes !! NetCDF dimension and variable names for the enounter save object type, extends(netcdf_parameters) :: symba_io_encounter_parameters - integer(I4B) :: COLLIDER_DIM_SIZE = 2 !! Size of collider dimension - 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) :: 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 @@ -230,7 +230,7 @@ module symba_classes type, extends(symba_nbody_system) :: symba_encounter_snapshot - integer(I4B) :: tslot !! The index for the time array in the final NetCDF file + 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 diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 0fceaa724..bc0f1b985 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -31,6 +31,8 @@ module subroutine symba_io_encounter_dump(self, param) self%nc%ienc_frame = i call snapshot%write_frame(self%nc,param) end select + else + exit end if end do @@ -57,7 +59,6 @@ module subroutine symba_io_encounter_initialize_output(self, param) character(len=STRMAX) :: errmsg integer(I4B) :: ndims - associate(nc => self) dfill = ieee_value(dfill, IEEE_QUIET_NAN) sfill = ieee_value(sfill, IEEE_QUIET_NAN) @@ -69,7 +70,6 @@ module subroutine symba_io_encounter_initialize_output(self, param) 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 @@ -80,9 +80,9 @@ module subroutine symba_io_encounter_initialize_output(self, param) 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, NF90_UNLIMITED, 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, NF90_UNLIMITED, 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%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 @@ -141,37 +141,58 @@ module subroutine symba_io_encounter_write_frame(self, nc, param) 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 + 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, n - character(len=NAMELEN) :: charstring + integer(I4B) :: i, tslot, idslot, old_mode, npl, ntp + character(len=NAMELEN) :: charstring - tslot = self%tslot 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) - n = size(pl%id(:)) - do i = 1, n + npl = pl%nbody + do i = 1, npl idslot = pl%id(i) - call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "symba_io_encounter_write_frame nf90_put_var time_varid" ) - call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[idslot]), "symba_io_encounter_write_frame nf90_put_var 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 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 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 body 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 body radius_varid" ) + 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 body 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 body rotx_varid" ) + 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 name_varid" ) + 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 particle_type_varid" ) + 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 @@ -388,10 +409,20 @@ module subroutine symba_io_stop_encounter(self, param, t) class(symba_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time ! Internals - !character(STRMAX) + integer(I4B) :: i ! 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 + write(self%encounter_history%nc%enc_file, '("encounter_",I0.6,".nc")') param%iloop + call self%encounter_history%nc%initialize(param) call self%encounter_history%dump(param) call self%encounter_history%nc%close() diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 20d295b9c..b96270fc9 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -1326,9 +1326,10 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) associate(npl => self%pl%nbody, ntp => self%tp%nbody) allocate(symba_encounter_snapshot :: snapshot) + snapshot%t = t - if (npl > 0) allocate(symba_pl :: snapshot%pl) - if (ntp > 0) allocate(symba_tp :: snapshot%tp) + allocate(symba_pl :: snapshot%pl) + allocate(symba_tp :: snapshot%tp) if (npl + ntp == 0) return npl_snap = npl ntp_snap = ntp @@ -1345,6 +1346,7 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) 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 @@ -1377,6 +1379,7 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) end if ! 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)) @@ -1390,10 +1393,6 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) end do end if - if (npl_snap + ntp_snap == 0) return - - snapshot%t = t - ! Save the snapshot self%encounter_history%iframe = self%encounter_history%iframe + 1 call self%resize_storage(self%encounter_history%iframe) From be401c7f95ddb3b4bdef112c5e0f97566cdb9bde Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 7 Dec 2022 12:12:50 -0500 Subject: [PATCH 05/13] Removed weird duplicate file --- src/setup/symba_collision.f90 | 1090 --------------------------------- 1 file changed, 1090 deletions(-) delete mode 100644 src/setup/symba_collision.f90 diff --git a/src/setup/symba_collision.f90 b/src/setup/symba_collision.f90 deleted file mode 100644 index c4d04ee75..000000000 --- a/src/setup/symba_collision.f90 +++ /dev/null @@ -1,1090 +0,0 @@ -!! 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 (symba_classes) s_symba_collision - use swiftest - -contains - - module function symba_collision_casedisruption(system, param, colliders, frag) result(status) - !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Create the fragments resulting from a non-catastrophic disruption collision - !! - implicit none - ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object - class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object - ! Result - integer(I4B) :: status !! Status flag assigned to this outcome - ! Internals - integer(I4B) :: i, nfrag - logical :: lfailure - character(len=STRMAX) :: message - - select case(frag%regime) - case(COLLRESOLVE_REGIME_DISRUPTION) - message = "Disruption between" - case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - message = "Supercatastrophic disruption between" - end select - call symba_collision_collider_message(system%pl, colliders%idx, message) - call io_log_one_message(FRAGGLE_LOG_OUT, message) - - ! Collisional fragments will be uniformly distributed around the pre-impact barycenter - call frag%set_mass_dist(colliders, param) - - ! Generate the position and velocity distributions of the fragments - call frag%generate_fragments(colliders, system, param, lfailure) - - if (lfailure) then - call io_log_one_message(FRAGGLE_LOG_OUT, "No fragment solution found, so treat as a pure hit-and-run") - status = ACTIVE - nfrag = 0 - select type(pl => system%pl) - class is (symba_pl) - pl%status(colliders%idx(:)) = status - pl%ldiscard(colliders%idx(:)) = .false. - pl%lcollision(colliders%idx(:)) = .false. - end select - else - ! Populate the list of new bodies - nfrag = frag%nbody - write(message, *) nfrag - call io_log_one_message(FRAGGLE_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") - select case(frag%regime) - case(COLLRESOLVE_REGIME_DISRUPTION) - status = DISRUPTION - case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - status = SUPERCATASTROPHIC - end select - frag%id(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)] - param%maxid = frag%id(nfrag) - call symba_collision_mergeaddsub(system, param, colliders, frag, status) - end if - - return - end function symba_collision_casedisruption - - - module function symba_collision_casehitandrun(system, param, colliders, frag) result(status) - !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Create the fragments resulting from a non-catastrophic hit-and-run collision - !! - implicit none - ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object - class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object - ! Result - integer(I4B) :: status !! Status flag assigned to this outcom - ! Internals - integer(I4B) :: i, ibiggest, nfrag, jtarg, jproj - logical :: lpure - character(len=STRMAX) :: message - - message = "Hit and run between" - call symba_collision_collider_message(system%pl, colliders%idx, message) - call io_log_one_message(FRAGGLE_LOG_OUT, trim(adjustl(message))) - - if (colliders%mass(1) > colliders%mass(2)) then - jtarg = 1 - jproj = 2 - else - jtarg = 2 - jproj = 1 - end if - - if (frag%mass_dist(2) > 0.9_DP * colliders%mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched - call io_log_one_message(FRAGGLE_LOG_OUT, "Pure hit and run. No new fragments generated.") - nfrag = 0 - lpure = .true. - else ! Imperfect hit and run, so we'll keep the largest body and destroy the other - lpure = .false. - call frag%set_mass_dist(colliders, param) - - ! Generate the position and velocity distributions of the fragments - call frag%generate_fragments(colliders, system, param, lpure) - - if (lpure) then - call io_log_one_message(FRAGGLE_LOG_OUT, "Should have been a pure hit and run instead") - nfrag = 0 - else - nfrag = frag%nbody - write(message, *) nfrag - call io_log_one_message(FRAGGLE_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") - end if - end if - if (lpure) then ! Reset these bodies back to being active so that nothing further is done to them - status = HIT_AND_RUN_PURE - select type(pl => system%pl) - class is (symba_pl) - pl%status(colliders%idx(:)) = ACTIVE - pl%ldiscard(colliders%idx(:)) = .false. - pl%lcollision(colliders%idx(:)) = .false. - end select - else - ibiggest = colliders%idx(maxloc(system%pl%Gmass(colliders%idx(:)), dim=1)) - frag%id(1) = system%pl%id(ibiggest) - frag%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] - param%maxid = frag%id(nfrag) - status = HIT_AND_RUN_DISRUPT - call symba_collision_mergeaddsub(system, param, colliders, frag, status) - end if - - return - end function symba_collision_casehitandrun - - - module function symba_collision_casemerge(system, param, colliders, frag) result(status) - !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Merge massive bodies. - !! - !! Adapted from David E. Kaufmann's Swifter routines symba_merge_pl.f90 and symba_discard_merge_pl.f90 - !! - !! Adapted from Hal Levison's Swift routines symba5_merge.f and discard_mass_merge.f - implicit none - ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object - class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object - ! Result - integer(I4B) :: status !! Status flag assigned to this outcome - ! Internals - integer(I4B) :: i, j, k, ibiggest - real(DP) :: pe - real(DP), dimension(NDIM) :: L_spin_new - character(len=STRMAX) :: message - - message = "Merging" - call symba_collision_collider_message(system%pl, colliders%idx, message) - call io_log_one_message(FRAGGLE_LOG_OUT, message) - - select type(pl => system%pl) - class is (symba_pl) - - call frag%set_mass_dist(colliders, param) - ibiggest = colliders%idx(maxloc(pl%Gmass(colliders%idx(:)), dim=1)) - frag%id(1) = pl%id(ibiggest) - frag%xb(:,1) = frag%xbcom(:) - frag%vb(:,1) = frag%vbcom(:) - - if (param%lrotation) then - ! Conserve angular momentum by putting pre-impact orbital momentum into spin of the new body - L_spin_new(:) = colliders%L_orbit(:,1) + colliders%L_orbit(:,2) + colliders%L_spin(:,1) + colliders%L_spin(:,2) - - ! Assume prinicpal axis rotation on 3rd Ip axis - frag%rot(:,1) = L_spin_new(:) / (frag%Ip(3,1) * frag%mass(1) * frag%radius(1)**2) - else ! If spin is not enabled, we will consider the lost pre-collision angular momentum as "escaped" and add it to our bookkeeping variable - param%Lescape(:) = param%Lescape(:) + colliders%L_orbit(:,1) + colliders%L_orbit(:,2) - end if - - ! Keep track of the component of potential energy due to the pre-impact colliders%idx for book-keeping - pe = 0.0_DP - do j = 1, colliders%ncoll - do i = j + 1, colliders%ncoll - pe = pe - pl%Gmass(i) * pl%mass(j) / norm2(pl%xb(:, i) - pl%xb(:, j)) - end do - end do - system%Ecollisions = system%Ecollisions + pe - system%Euntracked = system%Euntracked - pe - - ! Update any encounter lists that have the removed bodies in them so that they instead point to the new - do k = 1, system%plplenc_list%nenc - do j = 1, colliders%ncoll - i = colliders%idx(j) - if (i == ibiggest) cycle - if (system%plplenc_list%id1(k) == pl%id(i)) then - system%plplenc_list%id1(k) = pl%id(ibiggest) - system%plplenc_list%index1(k) = i - end if - if (system%plplenc_list%id2(k) == pl%id(i)) then - system%plplenc_list%id2(k) = pl%id(ibiggest) - system%plplenc_list%index2(k) = i - end if - if (system%plplenc_list%id1(k) == system%plplenc_list%id2(k)) system%plplenc_list%status(k) = INACTIVE - end do - end do - - status = MERGED - - call symba_collision_mergeaddsub(system, param, colliders, frag, status) - - end select - - return - end function symba_collision_casemerge - - - subroutine symba_collision_collider_message(pl, collidx, collider_message) - !! author: David A. Minton - !! - !! Prints a nicely formatted message about which bodies collided, including their names and ids. - !! This subroutine appends the body names and ids to an input message. - implicit none - ! Arguments - class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object - integer(I4B), dimension(:), intent(in) :: collidx !! Index of collisional colliders%idx members - character(*), intent(inout) :: collider_message !! The message to print to the screen. - ! Internals - integer(I4B) :: i, n - character(len=STRMAX) :: idstr - - n = size(collidx) - if (n == 0) return - - do i = 1, n - if (i > 1) collider_message = trim(adjustl(collider_message)) // " and " - collider_message = " " // trim(adjustl(collider_message)) // " " // trim(adjustl(pl%info(collidx(i))%name)) - write(idstr, '(I10)') pl%id(collidx(i)) - collider_message = trim(adjustl(collider_message)) // " (" // trim(adjustl(idstr)) // ") " - end do - - return - end subroutine symba_collision_collider_message - - - module function symba_collision_check_encounter(self, system, param, t, dt, irec) result(lany_collision) - !! author: David A. Minton - !! - !! Check for merger between massive bodies and test particles in SyMBA - !! - !! Adapted from David E. Kaufmann's Swifter routine symba_merge.f90 and symba_merge_tp.f90 - !! - !! Adapted from Hal Levison's Swift routine symba5_merge.f - implicit none - ! Arguments - class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter list object - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! current time - real(DP), intent(in) :: dt !! step size - integer(I4B), intent(in) :: irec !! Current recursion level - ! Result - logical :: lany_collision !! Returns true if cany pair of encounters resulted in a collision - ! Internals - logical, dimension(:), allocatable :: lcollision, lmask - real(DP), dimension(NDIM) :: xr, vr - integer(I4B) :: i, j, k, nenc - real(DP) :: rlim, Gmtot - logical :: isplpl - character(len=STRMAX) :: timestr, idstri, idstrj, message - class(symba_encounter), allocatable :: tmp - - lany_collision = .false. - if (self%nenc == 0) return - - select type(self) - class is (symba_plplenc) - isplpl = .true. - class default - isplpl = .false. - end select - - select type(pl => system%pl) - class is (symba_pl) - select type(tp => system%tp) - class is (symba_tp) - nenc = self%nenc - allocate(lmask(nenc)) - lmask(:) = ((self%status(1:nenc) == ACTIVE) .and. (pl%levelg(self%index1(1:nenc)) >= irec)) - if (isplpl) then - lmask(:) = lmask(:) .and. (pl%levelg(self%index2(1:nenc)) >= irec) - else - lmask(:) = lmask(:) .and. (tp%levelg(self%index2(1:nenc)) >= irec) - end if - if (.not.any(lmask(:))) return - - allocate(lcollision(nenc)) - lcollision(:) = .false. - - if (isplpl) then - do concurrent(k = 1:nenc, lmask(k)) - i = self%index1(k) - j = self%index2(k) - xr(:) = pl%rh(:, i) - pl%rh(:, j) - vr(:) = pl%vb(:, i) - pl%vb(:, j) - rlim = pl%radius(i) + pl%radius(j) - Gmtot = pl%Gmass(i) + pl%Gmass(j) - lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), & - Gmtot, rlim, dt, self%lvdotr(k)) - end do - else - do concurrent(k = 1:nenc, lmask(k)) - i = self%index1(k) - j = self%index2(k) - xr(:) = pl%rh(:, i) - tp%rh(:, j) - vr(:) = pl%vb(:, i) - tp%vb(:, j) - lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), & - pl%Gmass(i), pl%radius(i), dt, self%lvdotr(k)) - end do - end if - - lany_collision = any(lcollision(:)) - if (lany_collision) then - call pl%xh2xb(system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary - do k = 1, nenc - i = self%index1(k) - j = self%index2(k) - if (lcollision(k)) self%status(k) = COLLISION - self%t(k) = t - self%x1(:,k) = pl%rh(:,i) + system%cb%xb(:) - self%v1(:,k) = pl%vb(:,i) - if (isplpl) then - self%x2(:,k) = pl%rh(:,j) + system%cb%xb(:) - self%v2(:,k) = pl%vb(:,j) - if (lcollision(k)) then - ! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional colliders%idx - if (pl%lcollision(i) .or. pl%lcollision(j)) call pl%make_colliders([i,j]) - - ! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step - pl%lcollision([i, j]) = .true. - pl%status([i, j]) = COLLISION - call pl%info(i)%set_value(status="COLLISION", discard_time=t, discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i)) - call pl%info(j)%set_value(status="COLLISION", discard_time=t, discard_rh=pl%rh(:,j), discard_vh=pl%vh(:,j)) - end if - else - self%x2(:,k) = tp%rh(:,j) + system%cb%xb(:) - self%v2(:,k) = tp%vb(:,j) - if (lcollision(k)) then - tp%status(j) = DISCARDED_PLR - tp%ldiscard(j) = .true. - write(idstri, *) pl%id(i) - write(idstrj, *) tp%id(j) - write(timestr, *) t - call tp%info(j)%set_value(status="DISCARDED_PLR", discard_time=t, discard_rh=tp%rh(:,j), discard_vh=tp%vh(:,j)) - write(message, *) "Particle " // trim(adjustl(tp%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & - // " collided with massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & - // " at t = " // trim(adjustl(timestr)) - call io_log_one_message(FRAGGLE_LOG_OUT, message) - end if - end if - end do - end if - end select - end select - - ! Extract the pl-pl or pl-tp encounter list and return the pl-pl or pl-tp collision_list - if (lany_collision) then - select type(self) - class is (symba_plplenc) - call self%extract_collisions(system, param) - class default - allocate(tmp, mold=self) - call self%spill(tmp, lcollision, ldestructive=.true.) ! Remove this encounter pair from the encounter list - end select - end if - - return - end function symba_collision_check_encounter - - - pure elemental function symba_collision_check_one(xr, yr, zr, vxr, vyr, vzr, Gmtot, rlim, dt, lvdotr) result(lcollision) - !! author: David A. Minton - !! - !! Check for a merger between a single pair of particles - !! - !! Adapted from David E. Kaufmann's Swifter routines symba_merge_tp.f90 and symba_merge_pl.f90 - !! - !! Adapted from Hal Levison's Swift routine symba5_merge.f - implicit none - ! Arguments - real(DP), intent(in) :: xr, yr, zr !! Relative position vector components - real(DP), intent(in) :: vxr, vyr, vzr !! Relative velocity vector components - real(DP), intent(in) :: Gmtot !! Sum of G*mass of colliding bodies - real(DP), intent(in) :: rlim !! Collision limit - Typically the sum of the radii of colliding bodies - real(DP), intent(in) :: dt !! Step size - logical, intent(in) :: lvdotr !! Logical flag indicating that these two bodies are approaching in the current substep - ! Result - logical :: lcollision !! Logical flag indicating whether these two bodies will collide or not - ! Internals - real(DP) :: r2, rlim2, a, e, q, vdotr, tcr2, dt2 - - r2 = xr**2 + yr**2 + zr**2 - rlim2 = rlim**2 - - if (r2 <= rlim2) then ! checks if bodies are actively colliding in this time step - lcollision = .true. - else ! if they are not actively colliding in this time step, checks if they are going to collide next time step based on velocities and q - lcollision = .false. - vdotr = xr * vxr + yr * vyr + zr * vzr - if (lvdotr .and. (vdotr > 0.0_DP)) then - tcr2 = r2 / (vxr**2 + vyr**2 + vzr**2) - dt2 = dt**2 - if (tcr2 <= dt2) then - call orbel_xv2aeq(Gmtot, xr, yr, zr, vxr, vyr, vzr, a, e, q) - lcollision = (q < rlim) - end if - end if - end if - - return - end function symba_collision_check_one - - - function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, colliders) result(lflag) - !! author: David A. Minton - !! - !! Loops through the pl-pl collision list and groups families together by index. Outputs the indices of all colliders%idx members, - !! and pairs of quantities (x and v vectors, mass, radius, L_spin, and Ip) that can be used to resolve the collisional outcome. - implicit none - ! Arguments - class(symba_pl), intent(inout) :: pl !! SyMBA massive body object - class(symba_cb), intent(inout) :: cb !! SyMBA central body object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B), dimension(2), intent(inout) :: idx_parent !! Index of the two bodies considered the "parents" of the collision - class(fraggle_colliders), intent(out) :: colliders - ! Result - logical :: lflag !! Logical flag indicating whether a colliders%idx was successfully created or not - ! Internals - type collidx_array - integer(I4B), dimension(:), allocatable :: id - integer(I4B), dimension(:), allocatable :: idx - end type collidx_array - type(collidx_array), dimension(2) :: parent_child_index_array - integer(I4B), dimension(2) :: nchild - integer(I4B) :: i, j, ncolliders, idx_child - real(DP), dimension(2) :: volume, density - real(DP) :: mchild, volchild - real(DP), dimension(NDIM) :: xc, vc, xcom, vcom, xchild, vchild, xcrossv - real(DP), dimension(NDIM,2) :: mxc, vcc - - nchild(:) = pl%kin(idx_parent(:))%nchild - ! If all of these bodies share a parent, but this is still a unique collision, move the last child - ! out of the parent's position and make it the secondary body - if (idx_parent(1) == idx_parent(2)) then - if (nchild(1) == 0) then ! There is only one valid body recorded in this pair (this could happen due to restructuring of the kinship relationships, though it should be rare) - lflag = .false. - call pl%reset_kinship([idx_parent(1)]) - return - end if - idx_parent(2) = pl%kin(idx_parent(1))%child(nchild(1)) - nchild(1) = nchild(1) - 1 - nchild(2) = 0 - pl%kin(idx_parent(:))%nchild = nchild(:) - pl%kin(idx_parent(2))%parent = idx_parent(1) - end if - - colliders%mass(:) = pl%mass(idx_parent(:)) ! Note: This is meant to mass, not G*mass, as the collisional regime determination uses mass values that will be converted to Si - colliders%radius(:) = pl%radius(idx_parent(:)) - volume(:) = (4.0_DP / 3.0_DP) * PI * colliders%radius(:)**3 - - ! Group together the ids and indexes of each collisional parent and its children - do j = 1, 2 - allocate(parent_child_index_array(j)%idx(nchild(j)+ 1)) - allocate(parent_child_index_array(j)%id(nchild(j)+ 1)) - associate(idx_arr => parent_child_index_array(j)%idx, & - id_arr => parent_child_index_array(j)%id, & - ncj => nchild(j), & - pl => pl, & - plkinj => pl%kin(idx_parent(j))) - idx_arr(1) = idx_parent(j) - if (ncj > 0) idx_arr(2:ncj + 1) = plkinj%child(1:ncj) - id_arr(:) = pl%id(idx_arr(:)) - end associate - end do - - ! Consolidate the groups of collsional parents with any children they may have into a single "colliders%idx" index array - ncolliders = 2 + sum(nchild(:)) - allocate(colliders%idx(ncolliders)) - colliders%idx = [parent_child_index_array(1)%idx(:),parent_child_index_array(2)%idx(:)] - - colliders%ncoll = count(pl%lcollision(colliders%idx(:))) - colliders%idx = pack(colliders%idx(:), pl%lcollision(colliders%idx(:))) - colliders%L_spin(:,:) = 0.0_DP - colliders%Ip(:,:) = 0.0_DP - - ! Find the barycenter of each body along with its children, if it has any - do j = 1, 2 - colliders%xb(:, j) = pl%rh(:, idx_parent(j)) + cb%xb(:) - colliders%vb(:, j) = pl%vb(:, idx_parent(j)) - ! Assume principal axis rotation about axis corresponding to highest moment of inertia (3rd Ip) - if (param%lrotation) then - colliders%Ip(:, j) = colliders%mass(j) * pl%Ip(:, idx_parent(j)) - colliders%L_spin(:, j) = colliders%Ip(3, j) * colliders%radius(j)**2 * pl%rot(:, idx_parent(j)) - end if - - if (nchild(j) > 0) then - do i = 1, nchild(j) ! Loop over all children and take the mass weighted mean of the properties - idx_child = parent_child_index_array(j)%idx(i + 1) - if (.not. pl%lcollision(idx_child)) cycle - mchild = pl%mass(idx_child) - xchild(:) = pl%rh(:, idx_child) + cb%xb(:) - vchild(:) = pl%vb(:, idx_child) - volchild = (4.0_DP / 3.0_DP) * PI * pl%radius(idx_child)**3 - volume(j) = volume(j) + volchild - ! Get angular momentum of the child-parent pair and add that to the spin - ! Add the child's spin - if (param%lrotation) then - xcom(:) = (colliders%mass(j) * colliders%xb(:,j) + mchild * xchild(:)) / (colliders%mass(j) + mchild) - vcom(:) = (colliders%mass(j) * colliders%vb(:,j) + mchild * vchild(:)) / (colliders%mass(j) + mchild) - xc(:) = colliders%xb(:, j) - xcom(:) - vc(:) = colliders%vb(:, j) - vcom(:) - xcrossv(:) = xc(:) .cross. vc(:) - colliders%L_spin(:, j) = colliders%L_spin(:, j) + colliders%mass(j) * xcrossv(:) - - xc(:) = xchild(:) - xcom(:) - vc(:) = vchild(:) - vcom(:) - xcrossv(:) = xc(:) .cross. vc(:) - colliders%L_spin(:, j) = colliders%L_spin(:, j) + mchild * xcrossv(:) - - colliders%L_spin(:, j) = colliders%L_spin(:, j) + mchild * pl%Ip(3, idx_child) & - * pl%radius(idx_child)**2 & - * pl%rot(:, idx_child) - colliders%Ip(:, j) = colliders%Ip(:, j) + mchild * pl%Ip(:, idx_child) - end if - - ! Merge the child and parent - colliders%mass(j) = colliders%mass(j) + mchild - colliders%xb(:, j) = xcom(:) - colliders%vb(:, j) = vcom(:) - end do - end if - density(j) = colliders%mass(j) / volume(j) - colliders%radius(j) = (3 * volume(j) / (4 * PI))**(1.0_DP / 3.0_DP) - if (param%lrotation) colliders%Ip(:, j) = colliders%Ip(:, j) / colliders%mass(j) - end do - lflag = .true. - - xcom(:) = (colliders%mass(1) * colliders%xb(:, 1) + colliders%mass(2) * colliders%xb(:, 2)) / sum(colliders%mass(:)) - vcom(:) = (colliders%mass(1) * colliders%vb(:, 1) + colliders%mass(2) * colliders%vb(:, 2)) / sum(colliders%mass(:)) - mxc(:, 1) = colliders%mass(1) * (colliders%xb(:, 1) - xcom(:)) - mxc(:, 2) = colliders%mass(2) * (colliders%xb(:, 2) - xcom(:)) - vcc(:, 1) = colliders%vb(:, 1) - vcom(:) - vcc(:, 2) = colliders%vb(:, 2) - vcom(:) - colliders%L_orbit(:,:) = mxc(:,:) .cross. vcc(:,:) - - ! Destroy the kinship relationships for all members of this colliders%idx - call pl%reset_kinship(colliders%idx(:)) - - return - end function symba_collision_consolidate_colliders - - - module subroutine symba_collision_encounter_extract_collisions(self, system, param) - !! author: David A. Minton - !! - !! Processes the pl-pl encounter list remove only those encounters that led to a collision - !! - implicit none - ! Arguments - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - ! Internals - logical, dimension(:), allocatable :: lplpl_collision - logical, dimension(:), allocatable :: lplpl_unique_parent - integer(I4B), dimension(:), pointer :: plparent - integer(I4B), dimension(:), allocatable :: collision_idx, unique_parent_idx - integer(I4B) :: i, index_coll, ncollisions, nunique_parent, nplplenc - - select type (pl => system%pl) - class is (symba_pl) - associate(plplenc_list => self, idx1 => self%index1, idx2 => self%index2, plparent => pl%kin%parent) - nplplenc = plplenc_list%nenc - allocate(lplpl_collision(nplplenc)) - lplpl_collision(:) = plplenc_list%status(1:nplplenc) == COLLISION - if (.not.any(lplpl_collision)) return - ! Collisions have been detected in this step. So we need to determine which of them are between unique bodies. - - ! Get the subset of pl-pl encounters that lead to a collision - ncollisions = count(lplpl_collision(:)) - allocate(collision_idx(ncollisions)) - collision_idx = pack([(i, i=1, nplplenc)], lplpl_collision) - - ! Get the subset of collisions that involve a unique pair of parents - allocate(lplpl_unique_parent(ncollisions)) - - lplpl_unique_parent(:) = plparent(idx1(collision_idx(:))) /= plparent(idx2(collision_idx(:))) - nunique_parent = count(lplpl_unique_parent(:)) - allocate(unique_parent_idx(nunique_parent)) - unique_parent_idx = pack(collision_idx(:), lplpl_unique_parent(:)) - - ! Scrub all pl-pl collisions involving unique pairs of parents, which will remove all duplicates and leave behind - ! all pairs that have themselves as parents but are not part of the unique parent list. This can hapepn in rare cases - ! due to restructuring of parent/child relationships when there are large numbers of multi-body collisions in a single - ! step - lplpl_unique_parent(:) = .true. - do index_coll = 1, ncollisions - associate(ip1 => plparent(idx1(collision_idx(index_coll))), ip2 => plparent(idx2(collision_idx(index_coll)))) - lplpl_unique_parent(:) = .not. ( any(plparent(idx1(unique_parent_idx(:))) == ip1) .or. & - any(plparent(idx2(unique_parent_idx(:))) == ip1) .or. & - any(plparent(idx1(unique_parent_idx(:))) == ip2) .or. & - any(plparent(idx2(unique_parent_idx(:))) == ip2) ) - end associate - end do - - ! Reassemble collision index list to include only those containing the unique pairs of parents, plus all the non-unique pairs that don't - ! contain a parent body on the unique parent list. - ncollisions = nunique_parent + count(lplpl_unique_parent) - collision_idx = [unique_parent_idx(:), pack(collision_idx(:), lplpl_unique_parent(:))] - - ! Create a mask that contains only the pl-pl encounters that did not result in a collision, and then discard them - lplpl_collision(:) = .false. - lplpl_collision(collision_idx(:)) = .true. - call plplenc_list%spill(system%plplcollision_list, lplpl_collision, ldestructive=.true.) ! Extract any encounters that are not collisions from the list. - end associate - end select - - return - end subroutine symba_collision_encounter_extract_collisions - - - module subroutine symba_collision_make_colliders_pl(self, idx) - !! author: Jennifer L.L. Pouplin, Carlisle A. wishard, and David A. Minton - !! - !! When a single body is involved in more than one collision in a single step, it becomes part of a colliders%idx. - !! The largest body involved in a multi-body collision is the "parent" and all bodies that collide with it are its "children," - !! including those that collide with the children. - !! - !! Adapted from David E. Kaufmann's Swifter routine symba_merge_pl.f90 - !! - !! Adapted from Hal Levison's Swift routine symba5_merge.f - implicit none - ! Arguments - class(symba_pl), intent(inout) :: self !! SyMBA massive body object - integer(I4B), dimension(2), intent(in) :: idx !! Array holding the indices of the two bodies involved in the collision - ! Internals - integer(I4B) :: i, j, index_parent, index_child, p1, p2 - integer(I4B) :: nchild_inherit, nchild_orig, nchild_new - integer(I4B), dimension(:), allocatable :: temp - - associate(pl => self) - p1 = pl%kin(idx(1))%parent - p2 = pl%kin(idx(2))%parent - if (p1 == p2) return ! This is a collision between to children of a shared parent. We will ignore it. - - if (pl%mass(p1) > pl%mass(p2)) then - index_parent = p1 - index_child = p2 - else - index_parent = p2 - index_child = p1 - end if - - ! Expand the child array (or create it if necessary) and copy over the previous lists of children - nchild_orig = pl%kin(index_parent)%nchild - nchild_inherit = pl%kin(index_child)%nchild - nchild_new = nchild_orig + nchild_inherit + 1 - allocate(temp(nchild_new)) - - if (nchild_orig > 0) temp(1:nchild_orig) = pl%kin(index_parent)%child(1:nchild_orig) - ! Find out if the child body has any children of its own. The new parent wil inherit these children - if (nchild_inherit > 0) then - temp(nchild_orig+1:nchild_orig+nchild_inherit) = pl%kin(index_child)%child(1:nchild_inherit) - do i = 1, nchild_inherit - j = pl%kin(index_child)%child(i) - ! Set the childrens' parent to the new parent - pl%kin(j)%parent = index_parent - end do - end if - call pl%reset_kinship([index_child]) - ! Add the new child to its parent - pl%kin(index_child)%parent = index_parent - temp(nchild_new) = index_child - ! Save the new child array to the parent - pl%kin(index_parent)%nchild = nchild_new - call move_alloc(from=temp, to=pl%kin(index_parent)%child) - end associate - - return - end subroutine symba_collision_make_colliders_pl - - - subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) - !! author: David A. Minton - !! - !! Fills the pl_discards and pl_adds with removed and added bodies - !! - implicit none - ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object - class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object - integer(I4B), intent(in) :: status !! Status flag to assign to adds - ! Internals - integer(I4B) :: i, ibiggest, ismallest, iother, nstart, nend, ncolliders, nfrag - logical, dimension(system%pl%nbody) :: lmask - class(symba_pl), allocatable :: plnew, plsub - character(*), parameter :: FRAGFMT = '("Newbody",I0.7)' - character(len=NAMELEN) :: newname - - select type(pl => system%pl) - class is (symba_pl) - select type(pl_discards => system%pl_discards) - class is (symba_merger) - associate(info => pl%info, pl_adds => system%pl_adds, cb => system%cb, npl => pl%nbody) - ! Add the colliders%idx bodies to the subtraction list - ncolliders = colliders%ncoll - nfrag = frag%nbody - - param%maxid_collision = max(param%maxid_collision, maxval(system%pl%info(:)%collision_id)) - param%maxid_collision = param%maxid_collision + 1 - - ! Setup new bodies - allocate(plnew, mold=pl) - call plnew%setup(nfrag, param) - ibiggest = colliders%idx(maxloc(pl%Gmass(colliders%idx(:)), dim=1)) - ismallest = colliders%idx(minloc(pl%Gmass(colliders%idx(:)), dim=1)) - - ! Copy over identification, information, and physical properties of the new bodies from the fragment list - plnew%id(1:nfrag) = frag%id(1:nfrag) - plnew%xb(:, 1:nfrag) = frag%xb(:, 1:nfrag) - plnew%vb(:, 1:nfrag) = frag%vb(:, 1:nfrag) - call pl%vb2vh(cb) - call pl%xh2xb(cb) - do i = 1, nfrag - plnew%rh(:,i) = frag%xb(:, i) - cb%xb(:) - plnew%vh(:,i) = frag%vb(:, i) - cb%vb(:) - end do - plnew%mass(1:nfrag) = frag%mass(1:nfrag) - plnew%Gmass(1:nfrag) = param%GU * frag%mass(1:nfrag) - plnew%radius(1:nfrag) = frag%radius(1:nfrag) - plnew%density(1:nfrag) = frag%mass(1:nfrag) / frag%radius(1:nfrag) - call plnew%set_rhill(cb) - - select case(status) - case(DISRUPTION) - plnew%status(1:nfrag) = NEW_PARTICLE - do i = 1, nfrag - write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Disruption", origin_time=system%t, name=newname, & - origin_rh=plnew%rh(:,i), & - origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision) - end do - do i = 1, ncolliders - if (colliders%idx(i) == ibiggest) then - iother = ismallest - else - iother = ibiggest - end if - call pl%info(colliders%idx(i))%set_value(status="Disruption", discard_time=system%t, & - discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) - end do - case(SUPERCATASTROPHIC) - plnew%status(1:nfrag) = NEW_PARTICLE - do i = 1, nfrag - write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=system%t, name=newname, & - origin_rh=plnew%rh(:,i), origin_vh=plnew%vh(:,i), & - collision_id=param%maxid_collision) - end do - do i = 1, ncolliders - if (colliders%idx(i) == ibiggest) then - iother = ismallest - else - iother = ibiggest - end if - call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=system%t, & - discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), & - discard_body_id=iother) - end do - case(HIT_AND_RUN_DISRUPT) - call plnew%info(1)%copy(pl%info(ibiggest)) - plnew%status(1) = OLD_PARTICLE - do i = 2, nfrag - write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Hit and run fragment", origin_time=system%t, name=newname, & - origin_rh=plnew%rh(:,i), origin_vh=plnew%vh(:,i), & - collision_id=param%maxid_collision) - end do - do i = 1, ncolliders - if (colliders%idx(i) == ibiggest) cycle - iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status="Hit and run fragmention", discard_time=system%t, & - discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), & - discard_body_id=iother) - end do - case(MERGED) - call plnew%info(1)%copy(pl%info(ibiggest)) - plnew%status(1) = OLD_PARTICLE - do i = 1, ncolliders - if (colliders%idx(i) == ibiggest) cycle - - iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=system%t, discard_rh=pl%rh(:,i), & - discard_vh=pl%vh(:,i), discard_body_id=iother) - end do - end select - - if (param%lrotation) then - plnew%Ip(:, 1:nfrag) = frag%Ip(:, 1:nfrag) - plnew%rot(:, 1:nfrag) = frag%rot(:, 1:nfrag) - end if - - ! if (param%ltides) then - ! plnew%Q = pl%Q(ibiggest) - ! plnew%k2 = pl%k2(ibiggest) - ! plnew%tlag = pl%tlag(ibiggest) - ! end if - - !Copy over or set integration parameters for new bodies - plnew%lcollision(1:nfrag) = .false. - plnew%ldiscard(1:nfrag) = .false. - plnew%levelg(1:nfrag) = pl%levelg(ibiggest) - plnew%levelm(1:nfrag) = pl%levelm(ibiggest) - - ! Log the properties of the new bodies - call fraggle_io_log_pl(plnew, param) - - ! Append the new merged body to the list - nstart = pl_adds%nbody + 1 - nend = pl_adds%nbody + nfrag - call pl_adds%append(plnew, lsource_mask=[(.true., i=1, nfrag)]) - ! Record how many bodies were added in this event - pl_adds%ncomp(nstart:nend) = plnew%nbody - - ! Add the discarded bodies to the discard list - pl%status(colliders%idx(:)) = MERGED - pl%ldiscard(colliders%idx(:)) = .true. - pl%lcollision(colliders%idx(:)) = .true. - lmask(:) = .false. - lmask(colliders%idx(:)) = .true. - - call plnew%setup(0, param) - deallocate(plnew) - - allocate(plsub, mold=pl) - call pl%spill(plsub, lmask, ldestructive=.false.) - - nstart = pl_discards%nbody + 1 - nend = pl_discards%nbody + ncolliders - call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, ncolliders)]) - - ! Record how many bodies were subtracted in this event - pl_discards%ncomp(nstart:nend) = ncolliders - - call plsub%setup(0, param) - deallocate(plsub) - end associate - end select - end select - - return - end subroutine symba_collision_mergeaddsub - - - module subroutine symba_collision_resolve_fragmentations(self, system, param) - !! author: David A. Minton - !! - !! Process list of collisions, determine the collisional regime, and then create fragments. - !! - implicit none - ! Arguments - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - ! Internals - ! Internals - integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision - logical :: lgoodcollision - integer(I4B) :: i - type(fraggle_colliders) :: colliders !! Fraggle colliders object - type(fraggle_fragments) :: frag !! Fraggle fragmentation system object - - associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2) - select type(pl => system%pl) - class is (symba_pl) - select type (cb => system%cb) - class is (symba_cb) - do i = 1, ncollisions - idx_parent(1) = pl%kin(idx1(i))%parent - idx_parent(2) = pl%kin(idx2(i))%parent - lgoodcollision = symba_collision_consolidate_colliders(pl, cb, param, idx_parent, colliders) - if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLISION)) cycle - - call colliders%regime(frag, system, param) - - select case (frag%regime) - case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - plplcollision_list%status(i) = symba_collision_casedisruption(system, param, colliders, frag) - case (COLLRESOLVE_REGIME_HIT_AND_RUN) - plplcollision_list%status(i) = symba_collision_casehitandrun(system, param, colliders, frag) - case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - plplcollision_list%status(i) = symba_collision_casemerge(system, param, colliders, frag) - case default - write(*,*) "Error in symba_collision, unrecognized collision regime" - call util_exit(FAILURE) - end select - end do - end select - end select - end associate - - return - end subroutine symba_collision_resolve_fragmentations - - - module subroutine symba_collision_resolve_mergers(self, system, param) - !! author: David A. Minton - !! - !! Process list of collisions and merge colliding bodies together. - !! - implicit none - ! Arguments - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - ! Internals - integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision - logical :: lgoodcollision - integer(I4B) :: i - type(fraggle_colliders) :: colliders !! Fraggle colliders object - type(fraggle_fragments) :: frag !! Fraggle fragmentation system object - - associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2) - select type(pl => system%pl) - class is (symba_pl) - select type(cb => system%cb) - class is (symba_cb) - do i = 1, ncollisions - idx_parent(1) = pl%kin(idx1(i))%parent - idx_parent(2) = pl%kin(idx2(i))%parent - lgoodcollision = symba_collision_consolidate_colliders(pl, cb, param, idx_parent, colliders) - if (.not. lgoodcollision) cycle - if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved - - frag%regime = COLLRESOLVE_REGIME_MERGE - frag%mtot = sum(colliders%mass(:)) - frag%mass_dist(1) = frag%mtot - frag%mass_dist(2) = 0.0_DP - frag%mass_dist(3) = 0.0_DP - frag%xbcom(:) = (colliders%mass(1) * colliders%xb(:,1) + colliders%mass(2) * colliders%xb(:,2)) / frag%mtot - frag%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / frag%mtot - plplcollision_list%status(i) = symba_collision_casemerge(system, param, colliders, frag) - end do - end select - end select - end associate - - return - end subroutine symba_collision_resolve_mergers - - - module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, irec) - !! author: David A. Minton - !! - !! Process the pl-pl collision list, then modifiy the massive bodies based on the outcome of the collision - !! - implicit none - ! Arguments - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - real(DP), intent(in) :: t !! Current simulation time - real(DP), intent(in) :: dt !! Current simulation step size - integer(I4B), intent(in) :: irec !! Current recursion level - ! Internals - real(DP) :: Eorbit_before, Eorbit_after - logical :: lplpl_collision - character(len=STRMAX) :: timestr - class(symba_parameters), allocatable :: tmp_param - - associate(plplenc_list => self, plplcollision_list => system%plplcollision_list) - select type(pl => system%pl) - class is (symba_pl) - select type(param) - class is (symba_parameters) - if (plplcollision_list%nenc == 0) return ! No collisions to resolve - ! Make sure that the heliocentric and barycentric coordinates are consistent with each other - call pl%vb2vh(system%cb) - call pl%xh2xb(system%cb) - - ! Get the energy before the collision is resolved - if (param%lenergy) then - call system%get_energy_and_momentum(param) - Eorbit_before = system%te - end if - - do - write(timestr,*) t - call io_log_one_message(FRAGGLE_LOG_OUT, "") - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & - "***********************************************************") - call io_log_one_message(FRAGGLE_LOG_OUT, "Collision between massive bodies detected at time t = " // & - trim(adjustl(timestr))) - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & - "***********************************************************") - allocate(tmp_param, source=param) - if (param%lfragmentation) then - call plplcollision_list%resolve_fragmentations(system, param) - else - call plplcollision_list%resolve_mergers(system, param) - end if - - ! Destroy the collision list now that the collisions are resolved - call plplcollision_list%setup(0_I8B) - - if ((system%pl_adds%nbody == 0) .and. (system%pl_discards%nbody == 0)) exit - - ! Save the add/discard information to file - call system%write_discard(tmp_param) - - ! Rearrange the arrays: Remove discarded bodies, add any new bodies, resort, and recompute all indices and encounter lists - call pl%rearray(system, tmp_param) - - ! Destroy the add/discard list so that we don't append the same body multiple times if another collision is detected - call system%pl_discards%setup(0, param) - call system%pl_adds%setup(0, param) - deallocate(tmp_param) - - ! Check whether or not any of the particles that were just added are themselves in a collision state. This will generate a new plplcollision_list - lplpl_collision = plplenc_list%collision_check(system, param, t, dt, irec) - - if (.not.lplpl_collision) exit - end do - - if (param%lenergy) then - call system%get_energy_and_momentum(param) - Eorbit_after = system%te - system%Ecollisions = system%Ecollisions + (Eorbit_after - Eorbit_before) - end if - - end select - end select - end associate - - return - end subroutine symba_collision_resolve_plplenc - - - module subroutine symba_collision_resolve_pltpenc(self, system, param, t, dt, irec) - !! author: David A. Minton - !! - !! Process the pl-tp collision list, then modifiy the massive bodies based on the outcome of the collision - !! - implicit none - ! Arguments - class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - real(DP), intent(in) :: t !! Current simulation tim - real(DP), intent(in) :: dt !! Current simulation step size - integer(I4B), intent(in) :: irec !! Current recursion level - - ! Make sure coordinate systems are all synced up due to being inside the recursion at this point - call system%pl%vb2vh(system%cb) - call system%tp%vb2vh(system%cb%vb) - call system%pl%b2h(system%cb) - call system%tp%b2h(system%cb) - - ! Discard the collider - call system%tp%discard(system, param) - - return - end subroutine symba_collision_resolve_pltpenc - -end submodule s_symba_collision \ No newline at end of file From ea48c5faada9e2c2522123b9bfe4e756329ecc03 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 7 Dec 2022 12:22:17 -0500 Subject: [PATCH 06/13] Changed the naming scheme behavior for disruptions. For non-supercat, the largest body keeps its id and name --- examples/Fragmentation/Fragmentation_Movie.py | 2 +- src/symba/symba_collision.f90 | 41 ++++++++----------- 2 files changed, 18 insertions(+), 25 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index fa1203a8e..4e2ada041 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -97,7 +97,7 @@ def __init__(self, sim, animfile, title, style, nskip=1): self.body_color_list = {'Initial conditions': 'xkcd:windows blue', 'Disruption': 'xkcd:baby poop', 'Supercatastrophic': 'xkcd:shocking pink', - 'Hit and run fragment': 'xkcd:blue with a hint of purple', + 'Hit and run fragmention': 'xkcd:blue with a hint of purple', 'Central body': 'xkcd:almost black'} # Set up the figure and axes... diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 06c474078..eb891eb23 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -26,7 +26,7 @@ module function symba_collision_casedisruption(system, param, colliders, frag) ! Result integer(I4B) :: status !! Status flag assigned to this outcome ! Internals - integer(I4B) :: i, nfrag + integer(I4B) :: i, ibiggest, nfrag logical :: lfailure character(len=STRMAX) :: message @@ -63,11 +63,16 @@ module function symba_collision_casedisruption(system, param, colliders, frag) select case(frag%regime) case(COLLRESOLVE_REGIME_DISRUPTION) status = DISRUPTION + ibiggest = colliders%idx(maxloc(system%pl%Gmass(colliders%idx(:)), dim=1)) + frag%id(1) = system%pl%id(ibiggest) + frag%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] + param%maxid = frag%id(nfrag) case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) status = SUPERCATASTROPHIC + frag%id(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)] + param%maxid = frag%id(nfrag) end select - frag%id(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)] - param%maxid = frag%id(nfrag) + call symba_collision_mergeaddsub(system, param, colliders, frag, status) end if @@ -720,7 +725,7 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) logical, dimension(system%pl%nbody) :: lmask class(symba_pl), allocatable :: plnew, plsub character(*), parameter :: FRAGFMT = '("Newbody",I0.7)' - character(len=NAMELEN) :: newname + character(len=NAMELEN) :: newname, origin_type select type(pl => system%pl) class is (symba_pl) @@ -757,23 +762,6 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) call plnew%set_rhill(cb) select case(status) - case(DISRUPTION) - plnew%status(1:nfrag) = NEW_PARTICLE - do i = 1, nfrag - write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Disruption", origin_time=system%t, name=newname, & - origin_rh=plnew%rh(:,i), & - origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision) - end do - do i = 1, ncolliders - if (colliders%idx(i) == ibiggest) then - iother = ismallest - else - iother = ibiggest - end if - call pl%info(colliders%idx(i))%set_value(status="Disruption", discard_time=system%t, & - discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) - end do case(SUPERCATASTROPHIC) plnew%status(1:nfrag) = NEW_PARTICLE do i = 1, nfrag @@ -792,19 +780,24 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), & discard_body_id=iother) end do - case(HIT_AND_RUN_DISRUPT) + case(DISRUPTION,HIT_AND_RUN_DISRUPT) + if (status == DISRUPTION) then + write(origin_type,*) "Disruption" + else if (status == HIT_AND_RUN_DISRUPT) then + write(origin_type,*) "Hit and run fragmention" + end if call plnew%info(1)%copy(pl%info(ibiggest)) plnew%status(1) = OLD_PARTICLE do i = 2, nfrag write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Hit and run fragment", origin_time=system%t, name=newname, & + call plnew%info(i)%set_value(origin_type=origin_type, origin_time=system%t, name=newname, & origin_rh=plnew%rh(:,i), origin_vh=plnew%vh(:,i), & collision_id=param%maxid_collision) end do do i = 1, ncolliders if (colliders%idx(i) == ibiggest) cycle iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status="Hit and run fragmention", discard_time=system%t, & + call pl%info(colliders%idx(i))%set_value(status=origin_type, discard_time=system%t, & discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), & discard_body_id=iother) end do From 4d4417cf1a81bd27ed2c42f90e0a70abc98fb96e Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 7 Dec 2022 16:46:13 -0500 Subject: [PATCH 07/13] Added in code that can merge together the encounter and simulation data --- examples/Fragmentation/Fragmentation_Movie.py | 39 ++++++++++++++----- 1 file changed, 29 insertions(+), 10 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 4e2ada041..329a6f343 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-5, tstop=4.0e-3, istep_out=10, dump_cadence=0) + print("Generating animation") anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1) \ No newline at end of file From f7afdf34fa032de0d28ac1e899ad086d329d8974 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 7 Dec 2022 19:32:17 -0500 Subject: [PATCH 08/13] Improved the stability of the encounter snapshot --- examples/Fragmentation/Fragmentation_Movie.py | 6 ++-- src/modules/symba_classes.f90 | 17 ++++++++--- src/symba/symba_io.f90 | 10 +++---- src/symba/symba_util.f90 | 30 +++++++++++++++---- 4 files changed, 45 insertions(+), 18 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 329a6f343..ab2e2b34d 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -41,8 +41,8 @@ movie_titles = dict(zip(available_movie_styles, movie_title_list)) # These initial conditions were generated by trial and error -pos_vectors = {"disruption_headon" : [np.array([1.0, -5.0e-05, 0.0]), - np.array([1.0, 5.0e-05 ,0.0])], +pos_vectors = {"disruption_headon" : [np.array([1.0, -5.0e-04, 0.0]), + np.array([1.0, 5.0e-04 ,0.0])], "supercatastrophic_off_axis": [np.array([1.0, -4.2e-05, 0.0]), np.array([1.0, 4.2e-05, 0.0])], "hitandrun" : [np.array([1.0, -2.0e-05, 0.0]), @@ -203,7 +203,7 @@ def data_stream(self, frame=0): minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades sim.set_parameter(fragmentation=True, fragmentation_save="TRAJECTORY", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) - sim.run(dt=1e-5, tstop=4.0e-3, istep_out=10, dump_cadence=0) + sim.run(dt=1e-4, tstop=1.0e-4, istep_out=1, dump_cadence=0) print("Generating animation") anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1) \ No newline at end of file diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 6c3aaae12..4c1b6b673 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -229,11 +229,15 @@ module symba_classes end type symba_nbody_system - type, extends(symba_nbody_system) :: symba_encounter_snapshot - integer(I4B) :: tslot !! The index for the time array in the final NetCDF file + type :: symba_encounter_snapshot + type(symba_pl) :: pl !! Massive body data structure + type(symba_tp) :: tp !! Test particle data structure + real(DP) :: t = 0.0_DP !! Time at the snapshot + integer(I4B) :: tslot = 0 !! The index for the time array in the final NetCDF file contains - procedure :: write_encounter_frame => symba_io_encounter_write_frame !! Writes a frame of encounter data to file - generic :: write_frame => write_encounter_frame + procedure :: write_encounter_frame => symba_io_encounter_write_frame !! Writes a frame of encounter data to file + procedure :: dealloc => symba_util_dealloc_snapshot !! Deallocates all allocatable arrays + generic :: write_frame => write_encounter_frame !! Writes a snaphot frame to file final :: symba_util_final_encounter_snapshot end type symba_encounter_snapshot @@ -657,6 +661,11 @@ module subroutine symba_util_dealloc_merger(self) class(symba_merger), intent(inout) :: self !! SyMBA body merger object end subroutine symba_util_dealloc_merger + module subroutine symba_util_dealloc_snapshot(self) + implicit none + class(symba_encounter_snapshot), intent(inout) :: self !! SyMBA nbody system object + end subroutine symba_util_dealloc_snapshot + module subroutine symba_util_dealloc_system(self) implicit none class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index bc0f1b985..a4e5a9922 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -95,7 +95,8 @@ module subroutine symba_io_encounter_initialize_output(self, param) call check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], nc%ptype_varid), "symba_io_encounter_initialize_output nf90_def_var ptype_varid" ) call check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%rh_varid), "symba_io_encounter_initialize_output nf90_def_var rh_varid" ) call check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%vh_varid), "symba_io_encounter_initialize_output nf90_def_var vh_varid" ) - call check( nf90_def_var(nc%id, nc%gmass_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%Gmass_varid), "symba_io_encounter_initialize_output nf90_def_var Gmass_varid" ) + call check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%Gmass_varid), "symba_io_encounter_initialize_output nf90_def_var Gmass_varid" ) + call check( nf90_def_var(nc%id, nc%level_varname, NF90_INT, [nc%id_dimid, nc%time_dimid], nc%level_varid), "symba_io_encounter_initialize_output nf90_def_var level_varid" ) if (param%lclose) then call check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%radius_varid), "symba_io_encounter_initialize_output nf90_def_var radius_varid" ) end if @@ -153,8 +154,7 @@ module subroutine symba_io_encounter_write_frame(self, nc, param) tslot = self%tslot call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "symba_io_encounter_write_frame nf90_put_var time_varid" ) - select type(pl => self%pl) - class is (symba_pl) + associate(pl => self%pl, tp => self%tp) npl = pl%nbody do i = 1, npl idslot = pl%id(i) @@ -162,6 +162,7 @@ module subroutine symba_io_encounter_write_frame(self, nc, param) call check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var pl rh_varid" ) call check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var pl vh_varid" ) call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[idslot, tslot]), "symba_io_encounter_write_frame nf90_put_var pl Gmass_varid" ) + call check( nf90_put_var(nc%id, nc%level_varid, pl%levelg(i), start=[idslot, tslot]), "symba_io_encounter_write_frame nf90_put_var pl level_varid" ) if (param%lclose) call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[idslot, tslot]), "symba_io_encounter_write_frame nf90_put_var pl radius_varid" ) @@ -175,10 +176,7 @@ module subroutine symba_io_encounter_write_frame(self, nc, param) charstring = trim(adjustl(pl%info(i)%particle_type)) call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "symba_io_encounter_write_frame nf90_put_var pl particle_type_varid" ) end do - - end select - associate(tp => self%tp) ntp = tp%nbody do i = 1, ntp idslot = tp%id(i) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index b96270fc9..b5b74d93e 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -238,10 +238,28 @@ module subroutine symba_util_dealloc_merger(self) end subroutine symba_util_dealloc_merger + module subroutine symba_util_dealloc_snapshot(self) + !! author: David A. Minton + !! + !! Deallocates allocatable arrays in an encounter snapshot + implicit none + ! Arguments + class(symba_encounter_snapshot), intent(inout) :: self !! SyMBA nbody system object + + + call self%pl%dealloc() + call self%tp%dealloc() + self%t = 0.0_DP + self%tslot = 0 + + return + end subroutine symba_util_dealloc_snapshot + + module subroutine symba_util_dealloc_system(self) !! author: David A. Minton !! - !! Deallocates all allocatabale arrays + !! Deallocates all allocatabale arrays in a SyMBA system object implicit none ! Arguments class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object @@ -937,6 +955,8 @@ module subroutine symba_util_resize_storage(self, nnew) nbig = nbig * 2 end do allocate(symba_encounter_storage(nbig) :: tmp) + tmp%tvals(1:nold) = self%encounter_history%tvals(1:nold) + tmp%tvals(nold+1:nbig) = huge(1.0_DP) if (lmalloc) then do i = 1, nold if (allocated(self%encounter_history%frame(i)%item)) call move_alloc(self%encounter_history%frame(i)%item, tmp%frame(i)%item) @@ -1328,8 +1348,6 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) allocate(symba_encounter_snapshot :: snapshot) snapshot%t = t - allocate(symba_pl :: snapshot%pl) - allocate(symba_tp :: snapshot%tp) if (npl + ntp == 0) return npl_snap = npl ntp_snap = ntp @@ -1353,6 +1371,8 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) allocate(snapshot%pl%id(npl_snap)) allocate(snapshot%pl%info(npl_snap)) allocate(snapshot%pl%Gmass(npl_snap)) + allocate(snapshot%pl%levelg(npl_snap)) + snapshot%pl%levelg(:) = pack(pl%levelg(1:npl), pl%lmask(1:npl)) snapshot%pl%id(:) = pack(pl%id(1:npl), pl%lmask(1:npl)) snapshot%pl%info(:) = pack(pl%info(1:npl), pl%lmask(1:npl)) snapshot%pl%Gmass(:) = pack(pl%Gmass(1:npl), pl%lmask(1:npl)) @@ -1394,8 +1414,6 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) end if ! Save the snapshot - self%encounter_history%iframe = self%encounter_history%iframe + 1 - call self%resize_storage(self%encounter_history%iframe) ! Find out which time slot this belongs in by searching for an existing slot ! with the same value of time or the first available one @@ -1403,6 +1421,8 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) if (t <= self%encounter_history%tvals(i)) then snapshot%tslot = i self%encounter_history%tvals(i) = t + self%encounter_history%iframe = i + call self%resize_storage(i) exit end if end do From 706d439c1e07c79ad2c4ec88e8814be2805618d2 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 7 Dec 2022 19:32:46 -0500 Subject: [PATCH 09/13] Simplified fraggle.log output in preparation for storing detailed state information into a NetCDF file --- src/fraggle/fraggle_io.f90 | 124 ++----------------------------------- 1 file changed, 5 insertions(+), 119 deletions(-) diff --git a/src/fraggle/fraggle_io.f90 b/src/fraggle/fraggle_io.f90 index 87d2d7d11..a8c772326 100644 --- a/src/fraggle/fraggle_io.f90 +++ b/src/fraggle/fraggle_io.f90 @@ -36,31 +36,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 +57,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 +79,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) From 954764f42da41dfeb245b5eeac4fe6dbe24da9e4 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 7 Dec 2022 19:38:47 -0500 Subject: [PATCH 10/13] Minor rearranging in preparation for building the infrastructure for fraggle outputs. --- src/modules/symba_classes.f90 | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 4c1b6b673..cf6872846 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -183,7 +183,7 @@ module symba_classes !! 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(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 @@ -202,6 +202,20 @@ module symba_classes end type symba_encounter_storage + type :: symba_encounter_snapshot + !! A simplified version of a SyMBA nbody system object for storing minimal snapshots of the system state during encounters + type(symba_pl) :: pl !! Massive body data structure + type(symba_tp) :: tp !! Test particle data structure + real(DP) :: t = 0.0_DP !! Time at the snapshot + integer(I4B) :: tslot = 0 !! The index for the time array in the final NetCDF file + contains + procedure :: write_encounter_frame => symba_io_encounter_write_frame !! Writes a frame of encounter data to file + procedure :: dealloc => symba_util_dealloc_snapshot !! Deallocates all allocatable arrays + generic :: write_frame => write_encounter_frame !! Writes a snaphot frame to file + final :: symba_util_final_encounter_snapshot + end type symba_encounter_snapshot + + !******************************************************************************************************************************** ! symba_nbody_system class definitions and method interfaces !******************************************************************************************************************************** @@ -229,18 +243,6 @@ module symba_classes end type symba_nbody_system - type :: symba_encounter_snapshot - type(symba_pl) :: pl !! Massive body data structure - type(symba_tp) :: tp !! Test particle data structure - real(DP) :: t = 0.0_DP !! Time at the snapshot - integer(I4B) :: tslot = 0 !! The index for the time array in the final NetCDF file - contains - procedure :: write_encounter_frame => symba_io_encounter_write_frame !! Writes a frame of encounter data to file - procedure :: dealloc => symba_util_dealloc_snapshot !! Deallocates all allocatable arrays - generic :: write_frame => write_encounter_frame !! Writes a snaphot frame to file - final :: symba_util_final_encounter_snapshot - end type symba_encounter_snapshot - interface module function symba_collision_check_encounter(self, system, param, t, dt, irec) result(lany_collision) From 51d0853e361acb65f68341395be14fdf366b04aa Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 7 Dec 2022 19:53:44 -0500 Subject: [PATCH 11/13] Set temporary snapshot to be non-allocatable when saving --- src/symba/symba_util.f90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index b5b74d93e..41fd746ad 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -1340,12 +1340,11 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) 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(symba_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 if (npl + ntp == 0) return From ba6744137aa5cf16e423abc5632d3a6af3e74a6b Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 7 Dec 2022 20:03:18 -0500 Subject: [PATCH 12/13] Updated structure of the symba storage class to make the netcdf parameter a class instead of a type for easier extending --- src/modules/symba_classes.f90 | 5 ++--- src/setup/setup.f90 | 1 - src/symba/symba_io.f90 | 5 ++++- src/symba/symba_util.f90 | 1 + 4 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index cf6872846..c16644b1d 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -194,14 +194,13 @@ module symba_classes 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 + class(symba_io_encounter_parameters), allocatable :: 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 - type :: symba_encounter_snapshot !! A simplified version of a SyMBA nbody system object for storing minimal snapshots of the system state during encounters type(symba_pl) :: pl !! Massive body data structure 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 a4e5a9922..c5ed78f2a 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -384,7 +384,10 @@ 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(symba_encounter_storage :: self%encounter_history) + allocate(symba_io_encounter_parameters :: self%encounter_history%nc) + end if call self%encounter_history%reset() ! Empty out the time slot array for the next pass diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 41fd746ad..b66fe7ed1 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -955,6 +955,7 @@ module subroutine symba_util_resize_storage(self, nnew) nbig = nbig * 2 end do allocate(symba_encounter_storage(nbig) :: tmp) + call move_alloc(self%encounter_history%nc, tmp%nc) tmp%tvals(1:nold) = self%encounter_history%tvals(1:nold) tmp%tvals(nold+1:nbig) = huge(1.0_DP) if (lmalloc) then From ca87dd1c0a555440189e75e5bc35949b307cd67d Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 8 Dec 2022 12:53:54 -0500 Subject: [PATCH 13/13] Major restructuring of encounter storage classes and methods --- examples/Fragmentation/Fragmentation_Movie.py | 6 +- src/CMakeLists.txt | 2 + src/encounter/encounter_io.f90 | 207 +++++++++++++ src/encounter/encounter_util.f90 | 30 ++ src/fraggle/fraggle_io.f90 | 20 ++ src/fraggle/fraggle_util.f90 | 64 +++++ src/helio/helio_util.f90 | 2 +- src/io/io.f90 | 2 - src/main/swiftest_driver.f90 | 4 +- src/modules/encounter_classes.f90 | 67 ++++- src/modules/fraggle_classes.f90 | 66 ++++- src/modules/helio_classes.f90 | 6 +- src/modules/rmvs_classes.f90 | 10 +- src/modules/swiftest_classes.f90 | 72 ++--- src/modules/symba_classes.f90 | 98 +------ src/modules/whm_classes.f90 | 6 +- src/netcdf/netcdf.f90 | 18 +- src/rmvs/rmvs_util.f90 | 3 +- src/symba/symba_io.f90 | 200 +------------ src/symba/symba_util.f90 | 271 ++++++++---------- src/util/util_dealloc.f90 | 19 -- src/util/util_final.f90 | 62 ++++ src/util/util_reset.f90 | 2 + src/whm/whm_util.f90 | 2 +- 24 files changed, 704 insertions(+), 535 deletions(-) create mode 100644 src/encounter/encounter_io.f90 create mode 100644 src/util/util_final.f90 diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index ab2e2b34d..a664ef49e 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -41,8 +41,8 @@ movie_titles = dict(zip(available_movie_styles, movie_title_list)) # These initial conditions were generated by trial and error -pos_vectors = {"disruption_headon" : [np.array([1.0, -5.0e-04, 0.0]), - np.array([1.0, 5.0e-04 ,0.0])], +pos_vectors = {"disruption_headon" : [np.array([1.0, -5.0e-05, 0.0]), + np.array([1.0, 5.0e-05 ,0.0])], "supercatastrophic_off_axis": [np.array([1.0, -4.2e-05, 0.0]), np.array([1.0, 4.2e-05, 0.0])], "hitandrun" : [np.array([1.0, -2.0e-05, 0.0]), @@ -203,7 +203,7 @@ def data_stream(self, frame=0): minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades sim.set_parameter(fragmentation=True, fragmentation_save="TRAJECTORY", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) - sim.run(dt=1e-4, tstop=1.0e-4, 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 a8c772326..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 !! 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 c16644b1d..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,51 +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 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 - class(symba_io_encounter_parameters), allocatable :: 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 - - type :: symba_encounter_snapshot - !! A simplified version of a SyMBA nbody system object for storing minimal snapshots of the system state during encounters - type(symba_pl) :: pl !! Massive body data structure - type(symba_tp) :: tp !! Test particle data structure - real(DP) :: t = 0.0_DP !! Time at the snapshot - integer(I4B) :: tslot = 0 !! The index for the time array in the final NetCDF file - contains - procedure :: write_encounter_frame => symba_io_encounter_write_frame !! Writes a frame of encounter data to file - procedure :: dealloc => symba_util_dealloc_snapshot !! Deallocates all allocatable arrays - generic :: write_frame => write_encounter_frame !! Writes a snaphot frame to file - final :: symba_util_final_encounter_snapshot - end type symba_encounter_snapshot - - !******************************************************************************************************************************** ! 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 @@ -233,12 +198,10 @@ 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 @@ -424,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 @@ -662,16 +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_snapshot(self) - implicit none - class(symba_encounter_snapshot), intent(inout) :: self !! SyMBA nbody system object - end subroutine symba_util_dealloc_snapshot - - module subroutine symba_util_dealloc_system(self) - implicit none - class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - 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 @@ -721,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 @@ -799,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/symba/symba_io.f90 b/src/symba/symba_io.f90 index c5ed78f2a..e790e1685 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -11,192 +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" ) - call check( nf90_def_var(nc%id, nc%level_varname, NF90_INT, [nc%id_dimid, nc%time_dimid], nc%level_varid), "symba_io_encounter_initialize_output nf90_def_var level_varid" ) - if (param%lclose) then - call check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%radius_varid), "symba_io_encounter_initialize_output nf90_def_var radius_varid" ) - end if - 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" ) - - 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]), "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" ) - call check( nf90_put_var(nc%id, nc%level_varid, pl%levelg(i), start=[idslot, tslot]), "symba_io_encounter_write_frame nf90_put_var pl level_varid" ) - - if (param%lclose) call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[idslot, tslot]), "symba_io_encounter_write_frame nf90_put_var pl radius_varid" ) - - 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 - - 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 !! @@ -385,14 +199,10 @@ module subroutine symba_io_start_encounter(self, param, t) real(DP), intent(in) :: t !! Current simulation time if (.not. allocated(self%encounter_history)) then - allocate(symba_encounter_storage :: self%encounter_history) - allocate(symba_io_encounter_parameters :: self%encounter_history%nc) + 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) @@ -414,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 b66fe7ed1..d53bd442b 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -238,43 +238,6 @@ module subroutine symba_util_dealloc_merger(self) end subroutine symba_util_dealloc_merger - module subroutine symba_util_dealloc_snapshot(self) - !! author: David A. Minton - !! - !! Deallocates allocatable arrays in an encounter snapshot - implicit none - ! Arguments - class(symba_encounter_snapshot), intent(inout) :: self !! SyMBA nbody system object - - - call self%pl%dealloc() - call self%tp%dealloc() - self%t = 0.0_DP - self%tslot = 0 - - return - end subroutine symba_util_dealloc_snapshot - - - module subroutine symba_util_dealloc_system(self) - !! author: David A. Minton - !! - !! Deallocates all allocatabale arrays in a SyMBA system object - 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 !! @@ -489,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 !! @@ -511,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 - - call self%dealloc() - - return - end subroutine symba_util_final_encounter_snapshot - + 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) - 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 + call helio_util_final_system(self%helio_nbody_system) return - end subroutine symba_util_final_encounter_storage + end subroutine symba_util_final_system module subroutine symba_util_final_tp(self) @@ -924,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. @@ -933,43 +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) - call move_alloc(self%encounter_history%nc, tmp%nc) - tmp%tvals(1:nold) = self%encounter_history%tvals(1:nold) + allocate(encounter_storage(nbig) :: tmp) + tmp%tvals(1:nold) = system%encounter_history%tvals(1:nold) tmp%tvals(nold+1:nbig) = huge(1.0_DP) - if (lmalloc) then - do i = 1, nold - if (allocated(self%encounter_history%frame(i)%item)) call move_alloc(self%encounter_history%frame(i)%item, tmp%frame(i)%item) - end do - deallocate(self%encounter_history) - end if - call move_alloc(tmp,self%encounter_history) - self%encounter_history%iframe = iframe_old + 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) @@ -1338,10 +1293,10 @@ 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) :: snapshot + type(encounter_snapshot) :: snapshot integer(I4B) :: i, npl_snap, ntp_snap associate(npl => self%pl%nbody, ntp => self%tp%nbody) @@ -1356,81 +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)) - allocate(snapshot%pl%levelg(npl_snap)) - snapshot%pl%levelg(:) = pack(pl%levelg(1:npl), pl%lmask(1:npl)) - snapshot%pl%id(:) = pack(pl%id(1:npl), pl%lmask(1:npl)) - snapshot%pl%info(:) = pack(pl%info(1:npl), pl%lmask(1:npl)) - snapshot%pl%Gmass(:) = pack(pl%Gmass(1:npl), pl%lmask(1:npl)) - 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 - ! Save the snapshot + 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 - ! 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 - self%encounter_history%iframe = i - call self%resize_storage(i) - 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