diff --git a/examples/.gitignore b/examples/.gitignore index ad990dfcb..2baf6fec8 100644 --- a/examples/.gitignore +++ b/examples/.gitignore @@ -3,4 +3,5 @@ !Basic_Simulation !Fragmentation !helio_gr_test -!whm_gr_test \ No newline at end of file +!whm_gr_test +!Chambers2013 \ No newline at end of file diff --git a/examples/Chambers2013/.gitignore b/examples/Chambers2013/.gitignore new file mode 100644 index 000000000..3e7954819 --- /dev/null +++ b/examples/Chambers2013/.gitignore @@ -0,0 +1,4 @@ +* +!.gitignore +!init_cond.py +!aescattermovie.py \ No newline at end of file diff --git a/examples/Chambers2013/aescattermovie.py b/examples/Chambers2013/aescattermovie.py new file mode 100755 index 000000000..17361a343 --- /dev/null +++ b/examples/Chambers2013/aescattermovie.py @@ -0,0 +1,121 @@ +#!/usr/bin/env python3 +import swiftest +import numpy as np +import matplotlib.pyplot as plt +from matplotlib import animation +import matplotlib.colors as mcolors + +titletext = "Chambers (2013)" +radscale = 2000 +xmin = 0.0 +xmax = 2.20 +ymin = 0.0 +ymax = 1.0 +framejump = 1 + +class AnimatedScatter(): + """An animated scatter plot using matplotlib.animations.FuncAnimation.""" + def __init__(self, ds, param): + + frame = 0 + nframes = int(ds['time'].size / framejump) + self.ds = ds + self.param = param + self.Rcb = self.ds['radius'].sel(name="Sun").isel(time=0).values[()] + + self.clist = {'Initial conditions' : 'xkcd:faded blue', + 'Disruption' : 'xkcd:marigold', + 'Supercatastrophic' : 'xkcd:shocking pink', + 'Hit and run fragmentation' : 'xkcd:baby poop green'} + + # Setup the figure and axes... + fig = plt.figure(figsize=(8,4.5), dpi=300) + plt.tight_layout(pad=0) + # set up the figure + self.ax = plt.Axes(fig, [0.1, 0.15, 0.8, 0.75]) + self.ax.set_xlim(xmin, xmax) + self.ax.set_ylim(ymin, ymax) + fig.add_axes(self.ax) + self.ani = animation.FuncAnimation(fig, self.update, interval=1, frames=nframes, init_func=self.setup_plot, blit=True) + self.ani.save('aescatter.mp4', fps=60, dpi=300, extra_args=['-vcodec', 'libx264']) + print('Finished writing aescattter.mp4') + + def scatters(self, pl, radmarker, origin): + scat = [] + for key, value in self.clist.items(): + idx = origin == key + s = self.ax.scatter(pl[idx, 0], pl[idx, 1], marker='o', s=radmarker[idx], c=value, alpha=0.75, label=key) + scat.append(s) + return scat + + def setup_plot(self): + # First frame + """Initial drawing of the scatter plot.""" + t, name, Gmass, radius, npl, pl, radmarker, origin = next(self.data_stream(0)) + + # set up the figure + self.ax.margins(x=10, y=1) + self.ax.set_xlabel("Semimajor Axis (AU)", fontsize='16', labelpad=1) + self.ax.set_ylabel("Eccentricity", fontsize='16', labelpad=1) + + self.title = self.ax.text(0.50, 1.05, "", bbox={'facecolor': 'w', 'alpha': 0.5, 'pad': 5}, transform=self.ax.transAxes, + ha="center") + + self.title.set_text(f"{titletext} - Time = ${t*1e-6:6.2f}$ My with ${npl:4.0f}$ particles") + slist = self.scatters(pl, radmarker, origin) + self.s0 = slist[0] + self.s1 = slist[1] + self.s2 = slist[2] + self.s3 = slist[3] + leg = plt.legend(loc="upper right", scatterpoints=1, fontsize=10) + for i,l in enumerate(leg.legendHandles): + leg.legendHandles[i]._sizes = [20] + return self.s0, self.s1, self.s2, self.s3, self.title + + def data_stream(self, frame=0): + while True: + d = self.ds.isel(time = frame) + name_good = d.name.where(d['status'] != 1, drop=True) + name_good = name_good.where(name_good != "Sun", drop=True) + d = d.sel(name=name_good) + d['radmarker'] = (d['radius'] / self.Rcb) * radscale + radius = d['radmarker'].values + + radius = d['radmarker'].values + Gmass = d['Gmass'].values + a = d['a'].values + e = d['e'].values + name = d['name'].values + npl = d['npl'].values + radmarker = d['radmarker'] + origin = d['origin_type'] + + t = self.ds.coords['time'].values[frame] + + yield t, name, Gmass, radius, npl, np.c_[a, e], radmarker, origin + + def update(self,frame): + """Update the scatter plot.""" + t, name, Gmass, radius, npl, pl, radmarker, origin = next(self.data_stream(framejump * frame)) + + self.title.set_text(f"{titletext} - Time = ${t*1e-6:6.3f}$ My with ${npl:4.0f}$ particles") + + # We need to return the updated artist for FuncAnimation to draw.. + # Note that it expects a sequence of artists, thus the trailing comma. + s = [self.s0, self.s1, self.s2, self.s3] + for i, (key, value) in enumerate(self.clist.items()): + idx = origin == key + s[i].set_sizes(radmarker[idx]) + s[i].set_offsets(pl[idx,:]) + s[i].set_facecolor(value) + + self.s0 = s[0] + self.s1 = s[1] + self.s2 = s[2] + self.s3 = s[3] + return self.s0, self.s1, self.s2, self.s3, self.title, + +sim = swiftest.Simulation(simdir="fragglesim",read_old_output=True) +print('Making animation') +anim = AnimatedScatter(sim.data,sim.param) +print('Animation finished') diff --git a/examples/Chambers2013/init_cond.py b/examples/Chambers2013/init_cond.py new file mode 100755 index 000000000..c0f4b9b25 --- /dev/null +++ b/examples/Chambers2013/init_cond.py @@ -0,0 +1,56 @@ +#!/usr/bin/env python3 +import swiftest +import numpy as np +from numpy.random import default_rng + +# Initialize simulation object +sim = swiftest.Simulation(simdir="fragglesim") + +sim.set_parameter(compute_conservation_values=True, rotation=True, init_cond_format="EL",collision_model="fraggle",encounter_save="none") + +# Add bodies described in Chambers (2013) Sec. 2.1, with the uniform spatial distribution and two bodies sizes (big and small) +Nb = 14 +Ns = 140 +Mb = 2.8e-7 * 14 / Nb +Ms = 2.8e-8 * 140 / Ns +dens = 3000.0 / (sim.param['MU2KG'] / sim.param['DU2M']**3) +Rb = (3 * Mb / (4 * np.pi * dens) )**(1.0 / 3.0) +Rs = (3 * Ms / (4 * np.pi * dens) )**(1.0 / 3.0) +mtiny = 1e-2 * Ms +mininum_fragment_mass = 1e-4 * Ms +rng = default_rng(seed=880334) + +# Define the initial orbital elements of the big and small bodies +avalb = rng.uniform(0.3, 2.0, Nb) +avals = rng.uniform(0.3, 2.0, Ns) +evalb = rng.uniform(0.0, 0.01, Nb) +evals = rng.uniform(0.0, 0.01, Ns) +incvalb = rng.uniform(0.0, 0.005 * 180 / np.pi, Nb) +incvals = rng.uniform(0.0, 0.005 * 180 / np.pi, Ns) +capomvalb = rng.uniform(0.0, 360.0, Nb) +capomvals = rng.uniform(0.0, 360.0, Ns) +omegavalb = rng.uniform(0.0, 360.0, Nb) +omegavals = rng.uniform(0.0, 360.0, Ns) +capmvalb = rng.uniform(0.0, 360.0, Nb) +capmvals = rng.uniform(0.0, 360.0, Ns) +Ipvalb = np.full((Nb,3), 0.4) +Ipvals = np.full((Ns,3), 0.4) +rotvalb = np.zeros_like(Ipvalb) +rotvals = np.zeros_like(Ipvals) +GMvalb = np.full(Nb, Mb * sim.GU) +GMvals = np.full(Ns, Ms * sim.GU) +Rvalb = np.full(Nb, Rb) +Rvals = np.full(Ns, Rs) + +# Give the bodies unique names +nameb = [f"Big{i:03}" for i in range(Nb)] +names = [f"Small{i:03}" for i in range(Ns)] + +# Add the modern planets and the Sun using the JPL Horizons Database. +sim.add_solar_system_body(["Sun","Jupiter","Saturn","Uranus","Neptune"]) +sim.add_body(name=nameb, a=avalb, e=evalb, inc=incvalb, capom=capomvalb, omega=omegavalb, capm=capmvalb, Gmass=GMvalb, radius=Rvalb, rot=rotvalb, Ip=Ipvalb) +sim.add_body(name=names, a=avals, e=evals, inc=incvals, capom=capomvals, omega=omegavals, capm=capmvals, Gmass=GMvals, radius=Rvals, rot=rotvals, Ip=Ipvals) +sim.set_parameter(mtiny=mtiny, minimum_fragment_mass=mininum_fragment_mass) + +sim.run(tstop=3e8, dt=6.0875/365.25, istep_out=60000, dump_cadence=1) + diff --git a/examples/Fragmentation/swiftest_fragmentation.py b/examples/Fragmentation/swiftest_fragmentation.py index d533bee0a..417c01577 100644 --- a/examples/Fragmentation/swiftest_fragmentation.py +++ b/examples/Fragmentation/swiftest_fragmentation.py @@ -29,8 +29,8 @@ Each subdirectory contains: data.nc : A NetCDF file containing the simulation output. init_cond.nc : A NetCDF file containing the initial conditions for the simulation. -collision_000001.nc : A NetCDF file containing the data for the collision. -encounter_000001.nc : A NetCDF file containing the data for the close encounter. +collision_000001.nc : A NetCDF file containing the data for the collisions. +encounter_000001.nc : A NetCDF file containing the data for the close encounters. dump_bin1.nc : A NetCDF file containing the necessary inputs to restart a simulation from t!=0. dump_bin2.nc : A NetCDF file containing the necessary inputs to restart a simulation from t!=0. dump_param1.in : An ASCII file containing the necessary parameters to restart a simulation. diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index ade7cac00..e76c772aa 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -63,7 +63,7 @@ # handles strings differently than Python's Xarray. string_varnames = ["name", "particle_type", "status", "origin_type", "stage", "regime"] char_varnames = ["space"] -int_varnames = ["id", "ntp", "npl", "nplm", "discard_body_id", "collision_id", "loopnum"] +int_varnames = ["id", "ntp", "npl", "nplm", "discard_body_id", "collision_id", "status"] def bool2yesno(boolval): """ @@ -804,8 +804,6 @@ def process_netcdf_input(ds, param): Performs several tasks to convert raw NetCDF files output by the Fortran program into a form that is used by the Python side. These include: - Ensuring all types are correct - - Removing any bad id values (empty id slots) - - Swapping the id and name dimension if the names are unique Parameters ---------- diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index cfe88f26f..30fe46108 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -1039,7 +1039,7 @@ def set_feature(self, in initial conditions. encounter_save : {"NONE","TRAJECTORY","CLOSEST","BOTH"}, default "NONE" Indicate if and how encounter data should be saved. If set to "TRAJECTORY", the position and velocity vectors - of all bodies undergoing close encounters are saved at each intermediate step to the encounter files. + of all bodies undergoing close encounter are saved at each intermediate step to the encounter files. If set to "CLOSEST", the position and velocities at the point of closest approach between pairs of bodies are computed and stored to the encounter files. If set to "BOTH", then this stores the values that would be computed in "TRAJECTORY" and "CLOSEST". If set to "NONE" no trajectory information is saved. @@ -2738,8 +2738,8 @@ def read_output_file(self,read_init_cond : bool = True): else: self.init_cond = self.data.isel(time=0) - self.read_encounters() - self.read_collisions() + self.read_encounter() + self.read_collision() if self.verbose: print("Finished reading Swiftest dataset files.") @@ -2752,23 +2752,14 @@ def read_output_file(self,read_init_cond : bool = True): warnings.warn('Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".',stacklevel=2) return - def read_encounters(self): - enc_files = glob(f"{self.simdir}{os.path.sep}encounter_*.nc") - if len(enc_files) == 0: - return - - if self.verbose: - print("Reading encounter history file as .encounters") - - enc_files.sort() + def read_encounter(self): + enc_file = self.simdir / "encounters.nc" + if not os.path.exists(enc_file): + return - # This is needed in order to pass the param argument down to the io.process_netcdf_input function - def _preprocess(ds, param): - return io.process_netcdf_input(ds,param) - partial_func = partial(_preprocess, param=self.param) - - self.encounters = xr.open_mfdataset(enc_files,parallel=True,combine="nested",concat_dim="time",join="left",preprocess=partial_func,mask_and_scale=True) + self.encounters = xr.open_dataset(enc_file) self.encounters = io.process_netcdf_input(self.encounters, self.param) + # Remove any overlapping time values tgood,tid = np.unique(self.encounters.time,return_index=True) self.encounters = self.encounters.isel(time=tid) @@ -2778,21 +2769,16 @@ def _preprocess(ds, param): return - def read_collisions(self): - col_files = glob(f"{self.simdir}{os.path.sep}collision_*.nc") - if len(col_files) == 0: - return + def read_collision(self): + + col_file = self.simdir / "collisions.nc" + if not os.path.exists(col_file): + return - col_files.sort() if self.verbose: - print("Reading collision history file as .collisions") - - # This is needed in order to pass the param argument down to the io.process_netcdf_input function - def _preprocess(ds, param): - return io.process_netcdf_input(ds,param) - partial_func = partial(_preprocess, param=self.param) + print("Reading collisions history file as .collisions") - self.collisions = xr.open_mfdataset(col_files,parallel=True, combine="nested", concat_dim="collision", preprocess=partial_func,mask_and_scale=True) + self.collisions = xr.open_dataset(col_file) self.collisions = io.process_netcdf_input(self.collisions, self.param) return @@ -2964,17 +2950,13 @@ def clean(self): old_files = [self.simdir / self.param['BIN_OUT'], self.simdir / "fraggle.log", self.simdir / "swiftest.log", + self.simdir / "collisions.log", + self.simdir / "collisions.nc", + self.simdir / "encounters.nc", + self.simdir / "param.restart.in", ] - glob_files = [self.simdir.glob("**/dump_param?.in")] \ - + [self.simdir.glob("**/dump_bin?.nc")] \ - + [self.simdir.glob("**/encounter_*.nc")] \ - + [self.simdir.glob("**/collision_*.nc")] for f in old_files: if f.exists(): os.remove(f) - for g in glob_files: - for f in g: - if f.exists(): - os.remove(f) return diff --git a/src/base/base_module.f90 b/src/base/base_module.f90 index 732045594..0610ece9e 100644 --- a/src/base/base_module.f90 +++ b/src/base/base_module.f90 @@ -20,7 +20,7 @@ module base !> User defined parameters that are read in from the parameters input file. !> Each paramter is initialized to a default values. type, abstract :: base_parameters - character(len=:), allocatable :: integrator !! Symbolic name of the nbody integrator used + character(len=:), allocatable :: integrator !! Name of the nbody integrator used character(len=:), allocatable :: param_file_name !! The name of the parameter file integer(I4B) :: maxid = -1 !! The current maximum particle id number integer(I4B) :: maxid_collision = 0 !! The current maximum collision id number @@ -29,7 +29,6 @@ module base real(DP) :: tstop = -1.0_DP !! Integration stop time real(DP) :: dt = -1.0_DP !! Time step integer(I8B) :: iloop = 0_I8B !! Main loop counter - integer(I4B) :: ioutput = 1 !! Output counter character(STRMAX) :: incbfile = CB_INFILE !! Name of input file for the central body character(STRMAX) :: inplfile = PL_INFILE !! Name of input file for massive bodies character(STRMAX) :: intpfile = TP_INFILE !! Name of input file for test particles @@ -85,15 +84,15 @@ module base logical :: ltides = .false. !! Include tidal dissipation ! Initial values to pass to the energy report subroutine (usually only used in the case of a restart, otherwise these will be updated with initial conditions values) - real(DP) :: Eorbit_orig = 0.0_DP !! Initial orbital energy + real(DP) :: E_orbit_orig = 0.0_DP !! Initial orbital energy real(DP) :: GMtot_orig = 0.0_DP !! Initial system mass - real(DP), dimension(NDIM) :: Ltot_orig = 0.0_DP !! Initial total angular momentum vector - real(DP), dimension(NDIM) :: Lorbit_orig = 0.0_DP !! Initial orbital angular momentum - real(DP), dimension(NDIM) :: Lspin_orig = 0.0_DP !! Initial spin angular momentum vector - real(DP), dimension(NDIM) :: Lescape = 0.0_DP !! Angular momentum of bodies that escaped the system (used for bookeeping) + real(DP), dimension(NDIM) :: L_total_orig = 0.0_DP !! Initial total angular momentum vector + real(DP), dimension(NDIM) :: L_orbit_orig = 0.0_DP !! Initial orbital angular momentum + real(DP), dimension(NDIM) :: L_spin_orig = 0.0_DP !! Initial spin angular momentum vector + real(DP), dimension(NDIM) :: L_escape = 0.0_DP !! Angular momentum of bodies that escaped the system (used for bookeeping) real(DP) :: GMescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping) - real(DP) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions - real(DP) :: Euntracked = 0.0_DP !! Energy gained from system due to escaped bodies + real(DP) :: E_collisions = 0.0_DP !! Energy lost from system due to collisions + real(DP) :: E_untracked = 0.0_DP !! Energy gained from system due to escaped bodies logical :: lfirstenergy = .true. !! This is the first time computing energe logical :: lfirstkick = .true. !! Initiate the first kick in a symplectic step logical :: lrestart = .false. !! Indicates whether or not this is a restarted run diff --git a/src/collision/collision_check.f90 b/src/collision/collision_check.f90 index 1f88c98aa..b67526977 100644 --- a/src/collision/collision_check.f90 +++ b/src/collision/collision_check.f90 @@ -244,7 +244,7 @@ module subroutine collision_check_pltp(self, nbody_system, param, t, dt, irec, l 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 swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) end if end do diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index b946a387e..9a23ee4f4 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -55,6 +55,13 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t) associate(impactors => nbody_system%collider%impactors, fragments => nbody_system%collider%fragments) select case (impactors%regime) case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + + ! Manually save the before/after snapshots because this case doesn't use the mergeaddsub procedure + select type(before => self%before) + class is (swiftest_nbody_system) + allocate(before%pl, source=pl) + end select + nfrag = size(impactors%id(:)) do i = 1, nfrag j = impactors%id(i) @@ -68,19 +75,18 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t) pl%ldiscard(j) = .false. pl%lcollision(j) = .false. end do - select type(before => self%before) - class is (swiftest_nbody_system) + select type(after => self%after) class is (swiftest_nbody_system) - allocate(after%pl, source=before%pl) ! Be sure to save the pl so that snapshots still work - end select + allocate(after%pl, source=pl) end select + case (COLLRESOLVE_REGIME_HIT_AND_RUN) call self%hitandrun(nbody_system, param, t) case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) call self%merge(nbody_system, param, t) ! Use the default collision model, which is merge case default - write(*,*) "Error in swiftest_collision, unrecognized collision regime" + call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Error in swiftest_collision, unrecognized collision regime") call util_exit(FAILURE) end select end associate @@ -107,7 +113,6 @@ module subroutine collision_generate_hitandrun(self, nbody_system, param, t) ! Internals character(len=STRMAX) :: message - select type(nbody_system) class is (swiftest_nbody_system) select type(pl => nbody_system%pl) @@ -117,18 +122,21 @@ module subroutine collision_generate_hitandrun(self, nbody_system, param, t) call collision_io_collider_message(nbody_system%pl, impactors%id, message) call swiftest_io_log_one_message(COLLISION_LOG_OUT, trim(adjustl(message))) + ! Manually save the before/after snapshots because this case doesn't use the mergeaddsub procedure + select type(before => self%before) + class is (swiftest_nbody_system) + allocate(before%pl, source=pl) + end select + status = HIT_AND_RUN_PURE pl%status(impactors%id(:)) = ACTIVE pl%ldiscard(impactors%id(:)) = .false. pl%lcollision(impactors%id(:)) = .false. - ! Be sure to save the pl so that snapshots still work - select type(before => self%before) - class is (swiftest_nbody_system) + select type(after => self%after) class is (swiftest_nbody_system) - allocate(after%pl, source=before%pl) + allocate(after%pl, source=pl) end select - end select end associate end select @@ -154,7 +162,7 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) real(DP), intent(in) :: t !! The time of the collision ! Internals integer(I4B) :: i, j, k, ibiggest - real(DP), dimension(NDIM) :: Lspin_new + real(DP), dimension(NDIM) :: L_spin_new real(DP) :: volume, G character(len=STRMAX) :: message @@ -193,12 +201,12 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) if (param%lrotation) then do concurrent(i = 1:NDIM) fragments%Ip(i,1) = sum(impactors%mass(:) * impactors%Ip(i,:)) - Lspin_new(i) = sum(impactors%Lorbit(i,:) + impactors%Lorbit(i,:)) + L_spin_new(i) = sum(impactors%L_orbit(i,:) + impactors%L_orbit(i,:)) end do fragments%Ip(:,1) = fragments%Ip(:,1) / fragments%mass(1) - fragments%rot(:,1) = Lspin_new(:) / (fragments%Ip(3,1) * fragments%mass(1) * fragments%radius(1)**2) + fragments%rot(:,1) = L_spin_new(:) / (fragments%Ip(3,1) * fragments%mass(1) * fragments%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 - nbody_system%Lescape(:) = nbody_system%Lescape(:) + impactors%Lorbit(:,1) + impactors%Lorbit(:,2) + nbody_system%L_escape(:) = nbody_system%L_escape(:) + impactors%L_orbit(:,1) + impactors%L_orbit(:,2) end if ! The fragment trajectory will be the barycentric trajectory diff --git a/src/collision/collision_io.f90 b/src/collision/collision_io.f90 index 820f69af1..c4270d4a8 100644 --- a/src/collision/collision_io.f90 +++ b/src/collision/collision_io.f90 @@ -98,21 +98,14 @@ module subroutine collision_io_netcdf_dump(self, param) select type(nc => self%nc) class is (collision_netcdf_parameters) select type(param) - class is (swiftest_parameters) + class is (base_parameters) if (self%iframe > 0) then - nc%file_number = nc%file_number + 1 call self%make_index_map() - nc%event_dimsize = self%nt - nc%name_dimsize = self%nid - - write(nc%file_name, '("collision_",I0.6,".nc")') nc%file_number - call nc%initialize(param) - - do i = 1, self%nframes + call nc%open(param) + do i = 1, self%iframe if (allocated(self%frame(i)%item)) then select type(snapshot => self%frame(i)%item) class is (collision_snapshot) - param%ioutput = i call snapshot%write_frame(self,param) end select else @@ -175,70 +168,72 @@ module subroutine collision_io_netcdf_initialize_output(self, param) nc%lfile_is_open = .true. ! Dimensions - call netcdf_io_check( nf90_def_dim(nc%id, nc%event_dimname, nc%event_dimsize, nc%event_dimid), "collision_io_netcdf_initialize_output nf90_def_dim event_dimid" ) ! Dimension to store individual collision events - call netcdf_io_check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "collision_io_netcdf_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension - call netcdf_io_check( nf90_def_dim(nc%id, nc%name_dimname, nc%name_dimsize, nc%name_dimid), "collision_io_netcdf_initialize_output nf90_def_dim name_dimid" ) ! Dimension to store particle id numbers - call netcdf_io_check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "collision_io_netcdf_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) - call netcdf_io_check( nf90_def_dim(nc%id, nc%stage_dimname, 2, nc%stage_dimid), "collision_io_netcdf_initialize_output nf90_def_dim stage_dimid" ) ! Dimension for stage variables (aka "before" vs. "after" + call netcdf_io_check( nf90_def_dim(nc%id, nc%collision_id_varname, NF90_UNLIMITED, nc%collision_id_dimid), "collision_io_netcdf_initialize_output nf90_def_dim collision_id_dimid" ) ! Dimension to store individual collision events + call netcdf_io_check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "collision_io_netcdf_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension + call netcdf_io_check( nf90_def_dim(nc%id, nc%name_dimname, NF90_UNLIMITED, nc%name_dimid), "collision_io_netcdf_initialize_output nf90_def_dim name_dimid" ) ! Dimension to store particle id numbers + call netcdf_io_check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "collision_io_netcdf_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) + call netcdf_io_check( nf90_def_dim(nc%id, nc%stage_dimname, 2, nc%stage_dimid), "collision_io_netcdf_initialize_output nf90_def_dim stage_dimid" ) ! Dimension for stage variables (aka "before" vs. "after" ! Dimension coordinates - call netcdf_io_check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "collision_io_netcdf_initialize_output nf90_def_var space_varid" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%name_dimname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%name_varid), "collision_io_netcdf_initialize_output nf90_def_var name_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%stage_dimname, NF90_CHAR, [nc%str_dimid, nc%stage_dimid], nc%stage_varid), "collision_io_netcdf_initialize_output nf90_def_var stage_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%collision_id_varname, NF90_INT, [nc%collision_id_dimid], nc%collision_id_varid), "collision_io_netcdf_initialize_output nf90_def_var collision_id_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "collision_io_netcdf_initialize_output nf90_def_var space_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%name_dimname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%name_varid), "collision_io_netcdf_initialize_output nf90_def_var name_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%stage_dimname, NF90_CHAR, [nc%str_dimid, nc%stage_dimid], nc%stage_varid), "collision_io_netcdf_initialize_output nf90_def_var stage_varid" ) ! Variables call netcdf_io_check( nf90_def_var(nc%id, nc%id_varname, NF90_INT, nc%name_dimid, nc%id_varid), "collision_io_netcdf_initialize_output nf90_def_var id_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, & - nc%event_dimid, nc%time_varid), "collision_io_netcdf_initialize_output nf90_def_var time_varid" ) + nc%collision_id_dimid, nc%time_varid), "collision_io_netcdf_initialize_output nf90_def_var time_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%regime_varname, NF90_CHAR, & - [nc%str_dimid, nc%event_dimid], nc%regime_varid), "collision_io_netcdf_initialize_output nf90_def_var regime_varid") + [nc%str_dimid, nc%collision_id_dimid], nc%regime_varid), "collision_io_netcdf_initialize_output nf90_def_var regime_varid") call netcdf_io_check( nf90_def_var(nc%id, nc%Qloss_varname, nc%out_type, & - [ nc%event_dimid], nc%Qloss_varid), "collision_io_netcdf_initialize_output nf90_def_var Qloss_varid") + [ nc%collision_id_dimid], nc%Qloss_varid), "collision_io_netcdf_initialize_output nf90_def_var Qloss_varid") call netcdf_io_check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, & - [nc%str_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%ptype_varid), "collision_io_netcdf_initialize_output nf90_def_var ptype_varid") + [nc%str_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%ptype_varid), "collision_io_netcdf_initialize_output nf90_def_var ptype_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%loop_varname, NF90_INT, & - [ nc%event_dimid], nc%loop_varid), "collision_io_netcdf_initialize_output nf90_def_var loop_varid") call netcdf_io_check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type,& - [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%rh_varid), "collision_io_netcdf_initialize_output nf90_def_var rh_varid") + [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%rh_varid), "collision_io_netcdf_initialize_output nf90_def_var rh_varid") call netcdf_io_check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type,& - [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%vh_varid), "collision_io_netcdf_initialize_output nf90_def_var vh_varid") + [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%vh_varid), "collision_io_netcdf_initialize_output nf90_def_var vh_varid") call netcdf_io_check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type,& - [ nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%Gmass_varid), "collision_io_netcdf_initialize_output nf90_def_var Gmass_varid") + [ nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%Gmass_varid), "collision_io_netcdf_initialize_output nf90_def_var Gmass_varid") call netcdf_io_check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type,& - [ nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%radius_varid), "collision_io_netcdf_initialize_output nf90_def_var radius_varid") + [ nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%radius_varid), "collision_io_netcdf_initialize_output nf90_def_var radius_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type,& - [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%Ip_varid), "collision_io_netcdf_initialize_output nf90_def_var Ip_varid") + if (param%lrotation) then + call netcdf_io_check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type,& + [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%Ip_varid), "collision_io_netcdf_initialize_output nf90_def_var Ip_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type,& - [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%rot_varid), "collision_io_netcdf_initialize_output nf90_def_var rot_varid") - - if (param%lenergy) then + call netcdf_io_check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type,& + [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%rot_varid), "collision_io_netcdf_initialize_output nf90_def_var rot_varid") + end if + if (param%lenergy) then call netcdf_io_check( nf90_def_var(nc%id, nc%ke_orb_varname, nc%out_type,& - [ nc%stage_dimid, nc%event_dimid], nc%KE_orb_varid), "collision_io_netcdf_initialize_output nf90_def_var KE_orb_varid") + [ nc%stage_dimid, nc%collision_id_dimid], nc%KE_orb_varid), "collision_io_netcdf_initialize_output nf90_def_var KE_orb_varid") call netcdf_io_check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type,& - [ nc%stage_dimid, nc%event_dimid], nc%KE_spin_varid), "collision_io_netcdf_initialize_output nf90_def_var KE_spin_varid" ) + [ nc%stage_dimid, nc%collision_id_dimid], nc%KE_spin_varid), "collision_io_netcdf_initialize_output nf90_def_var KE_spin_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type,& - [ nc%stage_dimid, nc%event_dimid], nc%PE_varid), "collision_io_netcdf_initialize_output nf90_def_var PE_varid" ) + [ nc%stage_dimid, nc%collision_id_dimid], nc%PE_varid), "collision_io_netcdf_initialize_output nf90_def_var PE_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%be_varname, nc%out_type,& - [ nc%stage_dimid, nc%event_dimid], nc%BE_varid), "collision_io_netcdf_initialize_output nf90_def_var BE_varid" ) + [ nc%stage_dimid, nc%collision_id_dimid], nc%BE_varid), "collision_io_netcdf_initialize_output nf90_def_var BE_varid" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%L_orb_varname, nc%out_type, & - [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%L_orb_varid), "collision_io_netcdf_initialize_output nf90_def_var L_orb_varid" ) + if (param%lrotation) then + call netcdf_io_check( nf90_def_var(nc%id, nc%L_orbit_varname, nc%out_type, & + [ nc%space_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%L_orbit_varid), "collision_io_netcdf_initialize_output nf90_def_var L_orbit_varid" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%Lspin_varname, nc%out_type,& - [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%Lspin_varid), "collision_io_netcdf_initialize_output nf90_def_var Lspin_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%L_spin_varname, nc%out_type,& + [ nc%space_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%L_spin_varid), "collision_io_netcdf_initialize_output nf90_def_var L_spin_varid" ) + end if end if call netcdf_io_check( nf90_inquire(nc%id, nVariables=nvar), "collision_io_netcdf_initialize_output nf90_inquire nVariables" ) @@ -255,6 +250,7 @@ module subroutine collision_io_netcdf_initialize_output(self, param) call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "collision_io_netcdf_initialize_output nf90_def_var_fill NF90_CHAR" ) end select end do + ! Take the file out of define mode call netcdf_io_check( nf90_enddef(nc%id), "collision_io_netcdf_initialize_output nf90_enddef" ) @@ -274,6 +270,87 @@ module subroutine collision_io_netcdf_initialize_output(self, param) end subroutine collision_io_netcdf_initialize_output + module subroutine collision_io_netcdf_open(self, param, readonly) + !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton + !! + !! Opens a NetCDF file and does the variable inquiries to activate variable ids + use netcdf + implicit none + ! Arguments + class(collision_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(base_parameters), intent(in) :: param !! Current run configuration parameters + logical, optional, intent(in) :: readonly !! Logical flag indicating that this should be open read only + ! Internals + integer(I4B) :: mode + character(len=STRMAX) :: errmsg + logical fileExists + + mode = NF90_WRITE + if (present(readonly)) then + if (readonly) mode = NF90_NOWRITE + end if + + select type(param) + class is (base_parameters) + associate(nc => self) + + inquire(file=nc%file_name, exist=fileExists) + if (.not.fileExists) then + call nc%initialize(param) + return + end if + + write(errmsg,*) "collision_io_netcdf_open nf90_open ",trim(adjustl(nc%file_name)) + call netcdf_io_check( nf90_open(nc%file_name, mode, nc%id), errmsg) + self%lfile_is_open = .true. + + ! Dimensions + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%collision_id_varname, nc%time_dimid), "collision_io_netcdf_open nf90_inq_dimid time_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%space_dimname, nc%space_dimid), "collision_io_netcdf_open nf90_inq_dimid space_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%name_dimname, nc%name_dimid), "collision_io_netcdf_open nf90_inq_dimid name_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%str_dimname, nc%str_dimid), "collision_io_netcdf_open nf90_inq_dimid str_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%stage_dimname, nc%stage_dimid), "collision_io_netcdf_open nf90_inq_dimid stage_dimid" ) + + ! Dimension coordinates + call netcdf_io_check( nf90_inq_varid(nc%id, nc%collision_id_varname, nc%time_varid), "collision_io_netcdf_open nf90_inq_varid time_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%space_dimname, nc%space_varid), "collision_io_netcdf_open nf90_inq_varid space_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%name_dimname, nc%name_varid), "collision_io_netcdf_open nf90_inq_varid name_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%stage_dimname, nc%stage_varid), "collision_io_netcdf_open nf90_inq_varid stage_varid" ) + + ! Required Variables + call netcdf_io_check( nf90_inq_varid(nc%id, nc%id_varname, nc%id_varid), "collision_io_netcdf_open nf90_inq_varid name_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%time_dimname, nc%time_varid), "collision_io_netcdf_open nf90_inq_varid time_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%regime_varname, nc%regime_varid), "collision_io_netcdf_open nf90_inq_varid regime_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%Qloss_varname, nc%Qloss_varid), "collision_io_netcdf_open nf90_inq_varid Qloss_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%ptype_varname, nc%ptype_varid), "collision_io_netcdf_open nf90_inq_varid ptype_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%rh_varname, nc%rh_varid), "collision_io_netcdf_open nf90_inq_varid rh_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%vh_varname, nc%vh_varid), "collision_io_netcdf_open nf90_inq_varid vh_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%Gmass_varname, nc%Gmass_varid), "collision_io_netcdf_open nf90_inq_varid Gmass_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%radius_varname, nc%radius_varid), "collision_io_netcdf_open nf90_inq_varid radius_varid" ) + + if (param%lrotation) then + call netcdf_io_check( nf90_inq_varid(nc%id, nc%Ip_varname, nc%Ip_varid), "collision_io_netcdf_open nf90_inq_varid Ip_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%rot_varname, nc%rot_varid), "collision_io_netcdf_open nf90_inq_varid rot_varid" ) + end if + + if (param%lenergy) then + call netcdf_io_check( nf90_inq_varid(nc%id, nc%ke_orb_varname, nc%ke_orb_varid), "collision_io_netcdf_open nf90_inq_varid ke_orb_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%pe_varname, nc%pe_varid), "collision_io_netcdf_open nf90_inq_varid pe_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%be_varname, nc%be_varid), "collision_io_netcdf_open nf90_inq_varid be_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%L_orbit_varname, nc%L_orbit_varid), "collision_io_netcdf_open nf90_inq_varid L_orbit_varid" ) + if (param%lrotation) then + call netcdf_io_check( nf90_inq_varid(nc%id, nc%ke_spin_varname, nc%ke_spin_varid), "collision_io_netcdf_open nf90_inq_varid ke_spin_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%L_spin_varname, nc%L_spin_varid), "collision_io_netcdf_open nf90_inq_varid L_spin_varid" ) + end if + end if + + end associate + end select + + return + end subroutine collision_io_netcdf_open + + module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) !! author: David A. Minton !! @@ -291,11 +368,11 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) select type(nc => history%nc) class is (collision_netcdf_parameters) - associate(collider => self%collider, impactors => self%collider%impactors, fragments => self%collider%fragments, eslot => param%ioutput) - call netcdf_io_check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "collision_io_netcdf_write_frame_snapshot nf90_set_fill" ) + associate(collider => self%collider, impactors => self%collider%impactors, fragments => self%collider%fragments, eslot => self%collider%collision_id) + call netcdf_io_check( nf90_set_fill(nc%id, NF90_NOFILL, old_mode), "collision_io_netcdf_write_frame_snapshot nf90_set_fill" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%collision_id_varid, eslot, start=[eslot]), "collision_io_netcdf_write_frame_snapshot nf90_put_var collision_id_varid" ) call netcdf_io_check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[eslot]), "collision_io_netcdf_write_frame_snapshot nf90_put_var time_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%loop_varid, int(self%iloop,kind=I4B), start=[eslot]), "collision_io_netcdf_write_frame_snapshot nf90_put_varloop_varid" ) charstring = trim(adjustl(REGIME_NAMES(impactors%regime))) call netcdf_io_check( nf90_put_var(nc%id, nc%regime_varid, charstring, start=[1, eslot], count=[len(charstring), 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var regime_varid" ) @@ -325,19 +402,22 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) call netcdf_io_check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var vh_varid" ) call netcdf_io_check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[ idslot, stage, eslot]), "collision_io_netcdf_write_frame_snapshot nf90_put_var Gmass_varid" ) call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[ idslot, stage, eslot]), "collision_io_netcdf_write_frame_snapshot nf90_put_var radius_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var Ip_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var rotx_varid" ) + if (param%lrotation) then + call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var Ip_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var rotx_varid" ) + end if end do end do end select end select + if (param%lenergy) then call netcdf_io_check( nf90_put_var(nc%id, nc%ke_orb_varid, collider%ke_orbit(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var ke_orb_varid before" ) call netcdf_io_check( nf90_put_var(nc%id, nc%ke_spin_varid, collider%ke_spin(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var ke_spin_varid before" ) call netcdf_io_check( nf90_put_var(nc%id, nc%pe_varid, collider%pe(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var pe_varid before" ) call netcdf_io_check( nf90_put_var(nc%id, nc%be_varid, collider%be(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var pe_varid before" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%L_orb_varid, collider%Lorbit(:,:), start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var L_orb_varid before" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%Lspin_varid, collider%Lspin(:,:), start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var Lspin_varid before" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%L_orbit_varid, collider%L_orbit(:,:), start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var L_orbit_varid before" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%L_spin_varid, collider%L_spin(:,:), start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var L_spin_varid before" ) end if call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode) ) diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 3fa91a220..50e2bb75a 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -18,7 +18,8 @@ module collision implicit none public - character(len=*), parameter :: COLLISION_LOG_OUT = "collision.log" !! Name of log file for collision diagnostic information + character(len=*), parameter :: COLLISION_OUTFILE = 'collisions.nc' !! Name of NetCDF output file for collision information + character(len=*), parameter :: COLLISION_LOG_OUT = "collisions.log" !! Name of log file for collision diagnostic information !>Symbolic names for collisional outcomes from collresolve_resolve: integer(I4B), parameter :: COLLRESOLVE_REGIME_MERGE = 1 @@ -57,8 +58,8 @@ module collision real(DP), dimension(NDIM,2) :: rc !! Two-body equivalent position vectors of the collider bodies prior to collision in collision center of mass coordinates real(DP), dimension(NDIM,2) :: vc !! Two-body equivalent velocity vectors of the collider bodies prior to collision in collision center of mass coordinates real(DP), dimension(NDIM,2) :: rot !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision - real(DP), dimension(NDIM,2) :: Lspin !! Two-body equivalent spin angular momentum vectors of the collider bodies prior to collision - real(DP), dimension(NDIM,2) :: Lorbit !! Two-body equivalent orbital angular momentum vectors of the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: L_spin !! Two-body equivalent spin angular momentum vectors of the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: L_orbit !! Two-body equivalent orbital angular momentum vectors of the collider bodies prior to collision real(DP), dimension(NDIM,2) :: Ip !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision real(DP), dimension(2) :: Gmass !! Two-body equivalent G*mass of the collider bodies prior to the collision real(DP), dimension(2) :: mass !! Two-body equivalent mass of the collider bodies prior to the collision @@ -112,10 +113,10 @@ module collision real(DP), dimension(NDIM,nbody) :: t_unit !! Array of tangential direction unit vectors of individual fragments in the collisional coordinate frame real(DP), dimension(NDIM,nbody) :: n_unit !! Array of normal direction unit vectors of individual fragments in the collisional coordinate frame integer(I1B), dimension(nbody) :: origin_body !! Array of indices indicating which impactor body (1 or 2) the fragment originates from - real(DP), dimension(NDIM) :: Lorbit_tot !! Orbital angular momentum vector of all fragments - real(DP), dimension(NDIM) :: Lspin_tot !! Spin angular momentum vector of all fragments - real(DP), dimension(NDIM,nbody) :: Lorbit !! Orbital angular momentum vector of each individual fragment - real(DP), dimension(NDIM,nbody) :: Lspin !! Spin angular momentum vector of each individual fragment + real(DP), dimension(NDIM) :: L_orbit_tot !! Orbital angular momentum vector of all fragments + real(DP), dimension(NDIM) :: L_spin_tot !! Spin angular momentum vector of all fragments + real(DP), dimension(NDIM,nbody) :: L_orbit !! Orbital angular momentum vector of each individual fragment + real(DP), dimension(NDIM,nbody) :: L_spin !! Spin angular momentum vector of each individual fragment real(DP) :: ke_orbit_tot !! Orbital kinetic energy of all fragments real(DP) :: ke_spin_tot !! Spin kinetic energy of all fragments real(DP) :: pe !! Potential energy of all fragments @@ -138,27 +139,27 @@ module collision !! to resolve collision by defining extended types of encounters_impactors and/or encounetr_fragments !! !! The generate method for this class is the merge model. This allows any extended type to have access to the merge procedure by selecting the collision_basic parent class - class(collision_fragments(:)), allocatable :: fragments !! Object containing information on the pre-collision system - class(collision_impactors), allocatable :: impactors !! Object containing information on the post-collision system - class(base_nbody_system), allocatable :: before !! A snapshot of the subset of the nbody_system involved in the collision - class(base_nbody_system), allocatable :: after !! A snapshot of the subset of the nbody_system containing products of the collision - integer(I4B) :: status !! Status flag to pass to the collision list once the collision has been resolved + class(collision_fragments(:)), allocatable :: fragments !! Object containing information on the pre-collision system + class(collision_impactors), allocatable :: impactors !! Object containing information on the post-collision system + class(base_nbody_system), allocatable :: before !! A snapshot of the subset of the nbody_system involved in the collision + class(base_nbody_system), allocatable :: after !! A snapshot of the subset of the nbody_system containing products of the collision + integer(I4B) :: status !! Status flag to pass to the collision list once the collision has been resolved + integer(I4B) :: collision_id !! ID number of this collision event ! For the following variables, index 1 refers to the *entire* n-body nbody_system in its pre-collisional state and index 2 refers to the nbody_system in its post-collisional state - real(DP), dimension(NDIM,2) :: Lorbit !! Before/after orbital angular momentum - real(DP), dimension(NDIM,2) :: Lspin !! Before/after spin angular momentum - real(DP), dimension(NDIM,2) :: Ltot !! Before/after total nbody_system angular momentum + real(DP), dimension(NDIM,2) :: L_orbit !! Before/after orbital angular momentum + real(DP), dimension(NDIM,2) :: L_spin !! Before/after spin angular momentum + real(DP), dimension(NDIM,2) :: L_total !! Before/after total nbody_system angular momentum real(DP), dimension(2) :: ke_orbit !! Before/after orbital kinetic energy real(DP), dimension(2) :: ke_spin !! Before/after spin kinetic energy real(DP), dimension(2) :: pe !! Before/after potential energy real(DP), dimension(2) :: be !! Before/after binding energy - real(DP), dimension(2) :: Etot !! Before/after total nbody_system energy + real(DP), dimension(2) :: te !! Before/after total system energy contains procedure :: setup => collision_util_setup_collider !! Initializer for the encounter collision system and the before/after snapshots procedure :: setup_impactors => collision_util_setup_impactors_collider !! Initializer for the impactors for the encounter collision system. Deallocates old impactors before creating new ones procedure :: setup_fragments => collision_util_setup_fragments_collider !! Initializer for the fragments of the collision system. procedure :: add_fragments => collision_util_add_fragments_to_collider !! Add fragments to nbody_system - procedure :: construct_temporary_system => collision_util_construct_temporary_system !! Constructs temporary n-body nbody_system in order to compute pre- and post-impact energy and momentum procedure :: get_energy_and_momentum => collision_util_get_energy_momentum !! Calculates total nbody_system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) procedure :: reset => collision_util_reset_system !! Deallocates all allocatables procedure :: set_budgets => collision_util_set_budgets !! Sets the energy and momentum budgets of the fragments based on the collider value @@ -168,6 +169,7 @@ module collision procedure :: merge => collision_generate_merge !! Merges the impactors to make a single final body end type collision_basic + type, extends(collision_basic) :: collision_bounce contains procedure :: generate => collision_generate_bounce !! If a collision would result in a disruption, "bounce" the bodies instead. @@ -175,8 +177,6 @@ module collision end type collision_bounce - - !! NetCDF dimension and variable names for the enounter save object type, extends(encounter_netcdf_parameters) :: collision_netcdf_parameters integer(I4B) :: stage_dimid !! ID for the stage dimension @@ -184,10 +184,7 @@ module collision character(NAMELEN) :: stage_dimname = "stage" !! name of the stage dimension (before/after) character(len=6), dimension(2) :: stage_coords = ["before", "after"] !! The stage coordinate labels - character(NAMELEN) :: event_dimname = "collision" !! Name of collision event dimension - integer(I4B) :: event_dimid !! ID for the collision event dimension - integer(I4B) :: event_varid !! ID for the collision event variable - integer(I4B) :: event_dimsize = 0 !! Number of events + integer(I4B) :: collision_id_dimid !! ID for the collision event dimension character(NAMELEN) :: Qloss_varname = "Qloss" !! name of the energy loss variable integer(I4B) :: Qloss_varid !! ID for the energy loss variable @@ -195,6 +192,7 @@ module collision integer(I4B) :: regime_varid !! ID for the collision regime variable contains procedure :: initialize => collision_io_netcdf_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object + procedure :: open => collision_io_netcdf_open !! Opens an old file final :: collision_final_netcdf_parameters !! Finalizer closes the NetCDF file end type collision_netcdf_parameters @@ -276,6 +274,13 @@ module subroutine collision_io_netcdf_initialize_output(self, param) class(base_parameters), intent(in) :: param !! Current run configuration parameters end subroutine collision_io_netcdf_initialize_output + module subroutine collision_io_netcdf_open(self, param, readonly) + implicit none + class(collision_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(base_parameters), intent(in) :: param !! Current run configuration parameters + logical, optional, intent(in) :: readonly !! Logical flag indicating that this should be open read only + end subroutine collision_io_netcdf_open + module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) implicit none class(collision_snapshot), intent(in) :: self !! Swiftest encounter structure @@ -375,15 +380,6 @@ module subroutine collision_util_add_fragments_to_collider(self, nbody_system, p class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters end subroutine collision_util_add_fragments_to_collider - module subroutine collision_util_construct_temporary_system(self, nbody_system, param, tmpsys, tmpparam) - implicit none - class(collision_basic), intent(inout) :: self !! Collision system object - class(base_nbody_system), intent(in) :: nbody_system !! Original swiftest nbody system object - class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters - class(base_nbody_system), allocatable, intent(out) :: tmpsys !! Output temporary swiftest nbody system object - class(base_parameters), allocatable, intent(out) :: tmpparam !! Output temporary configuration run parameters - end subroutine collision_util_construct_temporary_system - module subroutine collision_util_get_angular_momentum(self) implicit none class(collision_fragments(*)), intent(inout) :: self !! Fragment system object diff --git a/src/collision/collision_regime.f90 b/src/collision/collision_regime.f90 index fb0f2f7da..90a5402e2 100644 --- a/src/collision/collision_regime.f90 +++ b/src/collision/collision_regime.f90 @@ -264,7 +264,7 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, regime = COLLRESOLVE_REGIME_SUPERCATASTROPHIC ! supercatastrophic Qloss = (c_star + 1.0_DP) * U_binding ! Qr - Qr_supercat else - write(*,*) "Error no regime found in symba_regime" + call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Error no regime found in symba_regime") end if end if diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index 21dacc918..a7236faf3 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -16,7 +16,7 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa !! author: David A. Minton !! !! Loops through the pl-pl collision list and groups families together by index. Outputs the indices of all impactors%id members, - !! and pairs of quantities (x and v vectors, mass, radius, Lspin, and Ip) that can be used to resolve the collisional outcome. + !! 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(collision_impactors), intent(out) :: self !! Collision impactors object @@ -83,7 +83,7 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa impactors%ncoll = count(pl%lcollision(impactors%id(:))) impactors%id = pack(impactors%id(:), pl%lcollision(impactors%id(:))) - impactors%Lspin(:,:) = 0.0_DP + impactors%L_spin(:,:) = 0.0_DP impactors%Ip(:,:) = 0.0_DP ! Find the barycenter of each body along with its children, if it has any @@ -93,7 +93,7 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa ! Assume principal axis rotation about axis corresponding to highest moment of inertia (3rd Ip) if (param%lrotation) then impactors%Ip(:, j) = impactors%mass(j) * pl%Ip(:, idx_parent(j)) - impactors%Lspin(:, j) = impactors%Ip(3, j) * impactors%radius(j)**2 * pl%rot(:, idx_parent(j)) + impactors%L_spin(:, j) = impactors%Ip(3, j) * impactors%radius(j)**2 * pl%rot(:, idx_parent(j)) end if if (nchild(j) > 0) then @@ -113,14 +113,14 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa xc(:) = impactors%rb(:, j) - xcom(:) vc(:) = impactors%vb(:, j) - vcom(:) xcrossv(:) = xc(:) .cross. vc(:) - impactors%Lspin(:, j) = impactors%Lspin(:, j) + impactors%mass(j) * xcrossv(:) + impactors%L_spin(:, j) = impactors%L_spin(:, j) + impactors%mass(j) * xcrossv(:) xc(:) = xchild(:) - xcom(:) vc(:) = vchild(:) - vcom(:) xcrossv(:) = xc(:) .cross. vc(:) - impactors%Lspin(:, j) = impactors%Lspin(:, j) + mchild * xcrossv(:) + impactors%L_spin(:, j) = impactors%L_spin(:, j) + mchild * xcrossv(:) - impactors%Lspin(:, j) = impactors%Lspin(:, j) + mchild * pl%Ip(3, idx_child) & + impactors%L_spin(:, j) = impactors%L_spin(:, j) + mchild * pl%Ip(3, idx_child) & * pl%radius(idx_child)**2 & * pl%rot(:, idx_child) impactors%Ip(:, j) = impactors%Ip(:, j) + mchild * pl%Ip(:, idx_child) @@ -144,7 +144,7 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa mxc(:, 2) = impactors%mass(2) * (impactors%rb(:, 2) - xcom(:)) vcc(:, 1) = impactors%vb(:, 1) - vcom(:) vcc(:, 2) = impactors%vb(:, 2) - vcom(:) - impactors%Lorbit(:,:) = mxc(:,:) .cross. vcc(:,:) + impactors%L_orbit(:,:) = mxc(:,:) .cross. vcc(:,:) ! Destroy the kinship relationships for all members of this impactors%id call pl%reset_kinship(impactors%id(:)) @@ -334,9 +334,6 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status) nimpactors = impactors%ncoll nfrag = fragments%nbody - param%maxid_collision = max(param%maxid_collision, maxval(nbody_system%pl%info(:)%collision_id)) - param%maxid_collision = param%maxid_collision + 1 - ! Setup new bodies allocate(plnew, mold=pl) call plnew%setup(nfrag, param) @@ -440,12 +437,6 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status) plnew%info(1:nfrag)%particle_type = PL_TYPE_NAME end where - ! Log the properties of the new bodies - select type(after => collider%after) - class is (swiftest_nbody_system) - allocate(after%pl, source=plnew) - end select - ! Append the new merged body to the list call pl_adds%append(plnew, lsource_mask=[(.true., i=1, nfrag)]) @@ -457,16 +448,22 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status) lmask(:) = .false. lmask(impactors%id(:)) = .true. - call plnew%setup(0, param) - deallocate(plnew) - allocate(plsub, mold=pl) call pl%spill(plsub, lmask, ldestructive=.false.) call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, nimpactors)]) - call plsub%setup(0, param) - deallocate(plsub) + ! Log the properties of the old and new bodies + select type(before => collider%before) + class is (swiftest_nbody_system) + call move_alloc(plsub, before%pl) + end select + + select type(after => collider%after) + class is (swiftest_nbody_system) + call move_alloc(plnew, after%pl) + end select + end associate end select @@ -490,14 +487,13 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) real(DP), intent(in) :: dt !! Current simulation step size integer(I4B), intent(in) :: irec !! Current recursion level ! Internals - real(DP) :: Eorbit_before, Eorbit_after + real(DP) :: E_orbit_before, E_orbit_after logical :: lplpl_collision character(len=STRMAX) :: timestr integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision logical :: lgoodcollision integer(I4B) :: i, loop, ncollisions integer(I4B), parameter :: MAXCASCADE = 1000 - real(DP) :: dpe select type (nbody_system) class is (swiftest_nbody_system) @@ -518,7 +514,7 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) ! Get the energy before the collision is resolved if (param%lenergy) then call nbody_system%get_energy_and_momentum(param) - Eorbit_before = nbody_system%te + E_orbit_before = nbody_system%te end if do loop = 1, MAXCASCADE @@ -538,16 +534,20 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) idx_parent(2) = pl%kin(idx2(i))%parent call impactors%consolidate(nbody_system, param, idx_parent, lgoodcollision) if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLIDED)) cycle + + ! Advance the collision id number and save it + param%maxid_collision = max(param%maxid_collision, maxval(nbody_system%pl%info(:)%collision_id)) + param%maxid_collision = param%maxid_collision + 1 + collider%collision_id = param%maxid_collision + + ! Get the collision regime call impactors%get_regime(nbody_system, param) call collision_history%take_snapshot(param,nbody_system, t, "before") - call nbody_system%get_energy_and_momentum(param) - + ! Generate the new bodies resulting from the collision call collider%generate(nbody_system, param, t) - call nbody_system%get_energy_and_momentum(param) - call collision_history%take_snapshot(param,nbody_system, t, "after") plpl_collision%status(i) = collider%status @@ -574,9 +574,8 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) if (.not.lplpl_collision) exit if (loop == MAXCASCADE) then - write(*,*) - write(*,*) "A runaway collisional cascade has been detected in collision_resolve_plpl." - write(*,*) "Consider reducing the step size or changing the parameters in the collisional model to reduce the number of fragments." + call swiftest_io_log_one_message(COLLISION_LOG_OUT,"An runaway collisional cascade has been detected in collision_resolve_plpl.") + call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Consider reducing the step size or changing the parameters in the collisional model to reduce the number of fragments.") call util_exit(FAILURE) end if end associate @@ -584,8 +583,8 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) if (param%lenergy) then call nbody_system%get_energy_and_momentum(param) - Eorbit_after = nbody_system%te - nbody_system%Ecollisions = nbody_system%Ecollisions + (Eorbit_after - Eorbit_before) + E_orbit_after = nbody_system%te + nbody_system%E_collisions = nbody_system%E_collisions + (E_orbit_after - E_orbit_before) end if end associate diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 7f073eecf..fda192a5b 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -61,67 +61,6 @@ module subroutine collision_util_add_fragments_to_collider(self, nbody_system, p end subroutine collision_util_add_fragments_to_collider - module subroutine collision_util_construct_temporary_system(self, nbody_system, param, tmpsys, tmpparam) - !! Author: David A. Minton - !! - !! Constructs a temporary internal system consisting of active bodies and additional fragments. This internal temporary system is used to calculate system energy with and without fragments - implicit none - ! Arguments - class(collision_basic), intent(inout) :: self !! Fraggle collision system object - class(base_nbody_system), intent(in) :: nbody_system !! Original swiftest nbody system object - class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters - class(base_nbody_system), allocatable, intent(out) :: tmpsys !! Output temporary swiftest nbody system object - class(base_parameters), allocatable, intent(out) :: tmpparam !! Output temporary configuration run parameters - ! Internals - logical, dimension(:), allocatable :: linclude - integer(I4B) :: npl_tot - ! The following are needed in order to deal with typing requirements - class(swiftest_nbody_system), allocatable :: tmpsys_local - class(swiftest_parameters), allocatable :: tmpparam_local - - select type(nbody_system) - class is (swiftest_nbody_system) - select type(param) - class is (swiftest_parameters) - associate(fragments => self%fragments, nfrag => self%fragments%nbody, pl => nbody_system%pl, npl => nbody_system%pl%nbody, cb => nbody_system%cb) - ! Set up a new system based on the original - if (allocated(tmpparam)) deallocate(tmpparam) - if (allocated(tmpsys)) deallocate(tmpsys) - allocate(tmpparam_local, source=param) - select type(tmpparam_local) - class is (swiftest_parameters) - tmpparam_local%system_history%nc%lfile_is_open = .false. - end select - call swiftest_util_setup_construct_system(tmpsys_local, tmpparam_local) - - ! No test particles necessary for energy/momentum calcs - call tmpsys_local%tp%setup(0, tmpparam_local) - - ! Replace the empty central body object with a copy of the original - deallocate(tmpsys_local%cb) - allocate(tmpsys_local%cb, source=cb) - - ! Make space for the fragments - npl_tot = npl + nfrag - call tmpsys_local%pl%setup(npl_tot, tmpparam_local) - allocate(linclude(npl_tot)) - - ! Fill up the temporary system with all of the original bodies, leaving the spaces for fragments empty until we add them in later - linclude(1:npl) = .true. - linclude(npl+1:npl_tot) = .false. - call tmpsys_local%pl%fill(pl, linclude) - - call move_alloc(tmpsys_local, tmpsys) - call move_alloc(tmpparam_local, tmpparam) - - end associate - end select - end select - - return - end subroutine collision_util_construct_temporary_system - - module subroutine collision_util_get_idvalues_snapshot(self, idvals) !! author: David A. Minton !! @@ -188,12 +127,12 @@ module subroutine collision_util_get_angular_momentum(self) associate(fragments => self, nfrag => self%nbody) do i = 1, nfrag - fragments%Lorbit(:,i) = fragments%mass(i) * (fragments%rc(:,i) .cross. fragments%vc(:, i)) - fragments%Lspin(:,i) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i) * fragments%rot(:,i) + fragments%L_orbit(:,i) = fragments%mass(i) * (fragments%rc(:,i) .cross. fragments%vc(:, i)) + fragments%L_spin(:,i) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i) * fragments%rot(:,i) end do - fragments%Lorbit_tot(:) = sum(fragments%Lorbit, dim=2) - fragments%Lspin_tot(:) = sum(fragments%Lspin, dim=2) + fragments%L_orbit_tot(:) = sum(fragments%L_orbit, dim=2) + fragments%L_spin_tot(:) = sum(fragments%L_spin, dim=2) end associate return @@ -253,9 +192,8 @@ module subroutine collision_util_get_energy_momentum(self, nbody_system, param, logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the nbody_system, with impactors included and fragments excluded or vice versa ! Internals integer(I4B) :: stage,i,j - real(DP), dimension(NDIM) :: Lorbit, Lspin + real(DP), dimension(NDIM) :: L_orbit, L_spin real(DP) :: ke_orbit, ke_spin, pe, be - real(DP), dimension(self%fragments%nbody) :: pepl select type(nbody_system) class is (swiftest_nbody_system) @@ -265,8 +203,8 @@ module subroutine collision_util_get_energy_momentum(self, nbody_system, param, if (lbefore) then - Lorbit(:) = sum(impactors%Lorbit(:,:), dim=2) - Lspin(:) = sum(impactors%Lspin(:,:), dim=2) + L_orbit(:) = sum(impactors%L_orbit(:,:), dim=2) + L_spin(:) = sum(impactors%L_spin(:,:), dim=2) ke_orbit = 0.0_DP ke_spin = 0.0_DP do concurrent(i = 1:2) @@ -281,8 +219,8 @@ module subroutine collision_util_get_energy_momentum(self, nbody_system, param, else call fragments%get_angular_momentum() - Lorbit(:) = fragments%Lorbit_tot(:) - Lspin(:) = fragments%Lspin_tot(:) + L_orbit(:) = fragments%L_orbit_tot(:) + L_spin(:) = fragments%L_spin_tot(:) call fragments%get_energy() ke_orbit = fragments%ke_orbit_tot @@ -297,14 +235,14 @@ module subroutine collision_util_get_energy_momentum(self, nbody_system, param, else stage = 2 end if - self%Lorbit(:,stage) = Lorbit(:) - self%Lspin(:,stage) = Lspin(:) - self%Ltot(:,stage) = Lorbit(:) + Lspin(:) + self%L_orbit(:,stage) = L_orbit(:) + self%L_spin(:,stage) = L_spin(:) + self%L_total(:,stage) = L_orbit(:) + L_spin(:) self%ke_orbit(stage) = ke_orbit self%ke_spin(stage) = ke_spin self%pe(stage) = pe self%be(stage) = be - self%Etot(stage) = ke_orbit + ke_spin + pe + be + self%te(stage) = ke_orbit + ke_spin + pe + be end associate end select end select @@ -351,8 +289,8 @@ module subroutine collision_util_reset_impactors(self) self%rb(:,:) = 0.0_DP self%vb(:,:) = 0.0_DP self%rot(:,:) = 0.0_DP - self%Lspin(:,:) = 0.0_DP - self%Lorbit(:,:) = 0.0_DP + self%L_spin(:,:) = 0.0_DP + self%L_orbit(:,:) = 0.0_DP self%Ip(:,:) = 0.0_DP self%mass(:) = 0.0_DP self%radius(:) = 0.0_DP @@ -420,13 +358,13 @@ module subroutine collision_util_reset_system(self) if (allocated(after%tp)) deallocate(after%tp) end select - self%Lorbit(:,:) = 0.0_DP - self%Lspin(:,:) = 0.0_DP - self%Ltot(:,:) = 0.0_DP + self%L_orbit(:,:) = 0.0_DP + self%L_spin(:,:) = 0.0_DP + self%L_total(:,:) = 0.0_DP self%ke_orbit(:) = 0.0_DP self%ke_spin(:) = 0.0_DP self%pe(:) = 0.0_DP - self%Etot(:) = 0.0_DP + self%te(:) = 0.0_DP if (allocated(self%impactors)) call self%impactors%reset() if (allocated(self%fragments)) deallocate(self%fragments) @@ -445,8 +383,8 @@ module subroutine collision_util_set_budgets(self) associate(impactors => self%impactors, fragments => self%fragments) - fragments%L_budget(:) = self%Ltot(:,1) - fragments%E_budget = self%Etot(1) - impactors%Qloss + fragments%L_budget(:) = self%L_total(:,1) + fragments%E_budget = self%te(1) - impactors%Qloss end associate @@ -509,7 +447,7 @@ module subroutine collision_util_set_coordinate_impactors(self) ! Arguments class(collision_impactors), intent(inout) :: self !! Collisional nbody_system ! Internals - real(DP), dimension(NDIM) :: delta_r, delta_v, Ltot + real(DP), dimension(NDIM) :: delta_r, delta_v, L_total real(DP) :: L_mag, mtot associate(impactors => self) @@ -521,11 +459,11 @@ module subroutine collision_util_set_coordinate_impactors(self) ! y-axis is the separation distance impactors%y_unit(:) = .unit.delta_r(:) - Ltot = impactors%Lorbit(:,1) + impactors%Lorbit(:,2) + impactors%Lspin(:,1) + impactors%Lspin(:,2) + L_total = impactors%L_orbit(:,1) + impactors%L_orbit(:,2) + impactors%L_spin(:,1) + impactors%L_spin(:,2) - L_mag = .mag.Ltot(:) + L_mag = .mag.L_total(:) if (L_mag > sqrt(tiny(L_mag))) then - impactors%z_unit(:) = .unit.Ltot(:) + impactors%z_unit(:) = .unit.L_total(:) else ! Not enough angular momentum to determine a z-axis direction. We'll just pick a random direction call random_number(impactors%z_unit(:)) impactors%z_unit(:) = .unit.impactors%z_unit(:) @@ -628,8 +566,8 @@ module subroutine collision_util_setup_fragments_collider(self, nfrag) self%fragments%density(:) = 0.0_DP self%fragments%rmag(:) = 0.0_DP self%fragments%vmag(:) = 0.0_DP - self%fragments%Lorbit_tot(:) = 0.0_DP - self%fragments%Lspin_tot(:) = 0.0_DP + self%fragments%L_orbit_tot(:) = 0.0_DP + self%fragments%L_spin_tot(:) = 0.0_DP self%fragments%L_budget(:) = 0.0_DP self%fragments%ke_orbit_tot = 0.0_DP self%fragments%ke_spin_tot = 0.0_DP @@ -727,8 +665,7 @@ module subroutine collision_util_snapshot(self, param, nbody_system, t, arg) real(DP), intent(in), optional :: t !! Time of snapshot if different from nbody_system time character(*), intent(in), optional :: arg !! "before": takes a snapshot just before the collision. "after" takes the snapshot just after the collision. ! Arguments - class(collision_snapshot), allocatable :: snapshot - class(swiftest_pl), allocatable :: pl + class(collision_snapshot), allocatable, save :: snapshot character(len=:), allocatable :: stage if (present(arg)) then @@ -744,27 +681,52 @@ module subroutine collision_util_snapshot(self, param, nbody_system, t, arg) select case(stage) case("before") ! Saves the states of the bodies involved in the collision before the collision is resolved - associate (idx => nbody_system%collider%impactors%id, ncoll => nbody_system%collider%impactors%ncoll) - allocate(pl, mold=nbody_system%pl) - call pl%setup(ncoll, param) - pl%id(:) = nbody_system%pl%id(idx(:)) - pl%Gmass(:) = nbody_system%pl%Gmass(idx(:)) - pl%radius(:) = nbody_system%pl%radius(idx(:)) - pl%rot(:,:) = nbody_system%pl%rot(:,idx(:)) - pl%Ip(:,:) = nbody_system%pl%Ip(:,idx(:)) - pl%rh(:,:) = nbody_system%pl%rh(:,idx(:)) - pl%vh(:,:) = nbody_system%pl%vh(:,idx(:)) - pl%info(:) = nbody_system%pl%info(idx(:)) - select type (before => nbody_system%collider%before) - class is (swiftest_nbody_system) - call move_alloc(pl, before%pl) - end select - end associate - case("after") allocate(collision_snapshot :: snapshot) allocate(snapshot%collider, source=nbody_system%collider) snapshot%t = t + + ! Get and record the energy of the system before the collision + call nbody_system%get_energy_and_momentum(param) + snapshot%collider%L_orbit(:,1) = nbody_system%L_orbit(:) + snapshot%collider%L_spin(:,1) = nbody_system%L_spin(:) + snapshot%collider%L_total(:,1) = nbody_system%L_total(:) + snapshot%collider%ke_orbit(1) = nbody_system%ke_orbit + snapshot%collider%ke_spin(1) = nbody_system%ke_spin + snapshot%collider%pe(1) = nbody_system%pe + snapshot%collider%be(1) = nbody_system%be + snapshot%collider%te(1) = nbody_system%te + + case("after") + ! Get record the energy of the sytem after the collision + call nbody_system%get_energy_and_momentum(param) + snapshot%collider%L_orbit(:,2) = nbody_system%L_orbit(:) + snapshot%collider%L_spin(:,2) = nbody_system%L_spin(:) + snapshot%collider%L_total(:,2) = nbody_system%L_total(:) + snapshot%collider%ke_orbit(2) = nbody_system%ke_orbit + snapshot%collider%ke_spin(2) = nbody_system%ke_spin + snapshot%collider%pe(2) = nbody_system%pe + snapshot%collider%be(2) = nbody_system%be + snapshot%collider%te(2) = nbody_system%te + + select type(before_snap => snapshot%collider%before ) + class is (swiftest_nbody_system) + select type(before_orig => nbody_system%collider%before) + class is (swiftest_nbody_system) + call move_alloc(before_orig%pl, before_snap%pl) + end select + end select + + select type(after_snap => snapshot%collider%after ) + class is (swiftest_nbody_system) + select type(after_orig => nbody_system%collider%after) + class is (swiftest_nbody_system) + call move_alloc(after_orig%pl, after_snap%pl) + end select + end select + + ! Save the snapshot for posterity call collision_util_save_snapshot(nbody_system%collision_history,snapshot) + deallocate(snapshot) case default write(*,*) "collision_util_snapshot requies either 'before' or 'after' passed to 'arg'" end select diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index cba31cfb8..7c76cd4ca 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -26,18 +26,13 @@ module subroutine encounter_io_netcdf_dump(self, param) class is (encounter_netcdf_parameters) if (self%iframe > 0) then ! Create and save the output files for this encounter and fragmentation - nc%file_number = nc%file_number + 1 call self%make_index_map() - nc%time_dimsize = self%nt - nc%name_dimsize = self%nid - write(nc%file_name, '("encounter_",I0.6,".nc")') nc%file_number - call nc%initialize(param) - + call nc%open(param) do i = 1, self%iframe if (allocated(self%frame(i)%item)) then select type(snapshot => self%frame(i)%item) class is (encounter_snapshot) - param%ioutput = self%tmap(i) + nc%tslot = nc%max_tslot + self%tmap(i) call snapshot%write_frame(self,param) end select else @@ -47,6 +42,8 @@ module subroutine encounter_io_netcdf_dump(self, param) call nc%close() call self%reset() + ! Update the time slot tracker + nc%max_tslot = nc%max_tslot + maxval(self%tmap(1:self%iframe)) end if end select @@ -97,9 +94,9 @@ module subroutine encounter_io_netcdf_initialize_output(self, param) nc%lfile_is_open = .true. ! Dimensions - call netcdf_io_check( nf90_def_dim(nc%id, nc%time_dimname, nc%time_dimsize, nc%time_dimid), "encounter_io_netcdf_initialize_output nf90_def_dim time_dimid" ) ! Simulation time dimension + call netcdf_io_check( nf90_def_dim(nc%id, nc%time_dimname, NF90_UNLIMITED, nc%time_dimid), "encounter_io_netcdf_initialize_output nf90_def_dim time_dimid" ) ! Simulation time dimension call netcdf_io_check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "encounter_io_netcdf_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension - call netcdf_io_check( nf90_def_dim(nc%id, nc%name_dimname, nc%name_dimsize, nc%name_dimid), "encounter_io_netcdf_initialize_output nf90_def_dim name_dimid" ) ! dimension to store particle id numbers + call netcdf_io_check( nf90_def_dim(nc%id, nc%name_dimname, NF90_UNLIMITED, nc%name_dimid), "encounter_io_netcdf_initialize_output nf90_def_dim name_dimid" ) ! dimension to store particle id numbers call netcdf_io_check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "encounter_io_netcdf_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) ! Dimension coordinates @@ -113,7 +110,6 @@ module subroutine encounter_io_netcdf_initialize_output(self, param) call netcdf_io_check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%rh_varid), "encounter_io_netcdf_initialize_output nf90_def_var rh_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%vh_varid), "encounter_io_netcdf_initialize_output nf90_def_var vh_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%Gmass_varid), "encounter_io_netcdf_initialize_output nf90_def_var Gmass_varid" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%loop_varname, NF90_INT, [nc%time_dimid], nc%loop_varid), "encounter_io_netcdf_initialize_output nf90_def_var loop_varid" ) if (param%lclose) then call netcdf_io_check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%radius_varid), "encounter_io_netcdf_initialize_output nf90_def_var radius_varid" ) end if @@ -153,6 +149,76 @@ module subroutine encounter_io_netcdf_initialize_output(self, param) end subroutine encounter_io_netcdf_initialize_output + module subroutine encounter_io_netcdf_open(self, param, readonly) + !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton + !! + !! Opens a NetCDF file and does the variable inquiries to activate variable ids + use netcdf + implicit none + ! Arguments + class(encounter_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(base_parameters), intent(in) :: param !! Current run configuration parameters + logical, optional, intent(in) :: readonly !! Logical flag indicating that this should be open read only + ! Internals + integer(I4B) :: mode + character(len=STRMAX) :: errmsg + logical fileExists + + mode = NF90_WRITE + if (present(readonly)) then + if (readonly) mode = NF90_NOWRITE + end if + + select type(param) + class is (base_parameters) + associate(nc => self) + + inquire(file=nc%file_name, exist=fileExists) + if (.not.fileExists) then + call nc%initialize(param) + return + end if + + write(errmsg,*) "encounter_io_netcdf_open nf90_open ",trim(adjustl(nc%file_name)) + call netcdf_io_check( nf90_open(nc%file_name, mode, nc%id), errmsg) + self%lfile_is_open = .true. + + ! Dimensions + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%time_dimname, nc%time_dimid), "encounter_io_netcdf_open nf90_inq_dimid time_dimid" ) + call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%time_dimid, nc%time_dimname, len=nc%max_tslot), "encounter_io_netcdf_open nf90_inquire_dimension max_tslot" ) + + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%space_dimname, nc%space_dimid), "encounter_io_netcdf_open nf90_inq_dimid space_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%name_dimname, nc%name_dimid), "encounter_io_netcdf_open nf90_inq_dimid name_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%str_dimname, nc%str_dimid), "encounter_io_netcdf_open nf90_inq_dimid str_dimid" ) + + ! Dimension coordinates + call netcdf_io_check( nf90_inq_varid(nc%id, nc%time_dimname, nc%time_varid), "encounter_io_netcdf_open nf90_inq_varid time_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%space_dimname, nc%space_varid), "encounter_io_netcdf_open nf90_inq_varid space_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%name_dimname, nc%name_varid), "encounter_io_netcdf_open nf90_inq_varid name_varid" ) + + ! Required Variables + call netcdf_io_check( nf90_inq_varid(nc%id, nc%id_varname, nc%id_varid), "encounter_io_netcdf_open nf90_inq_varid name_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%ptype_varname, nc%ptype_varid), "encounter_io_netcdf_open nf90_inq_varid ptype_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%time_dimname, nc%time_varid), "encounter_io_netcdf_open nf90_inq_varid time_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%rh_varname, nc%rh_varid), "encounter_io_netcdf_open nf90_inq_varid rh_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%vh_varname, nc%vh_varid), "encounter_io_netcdf_open nf90_inq_varid vh_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%Gmass_varname, nc%Gmass_varid), "encounter_io_netcdf_open nf90_inq_varid Gmass_varid" ) + if (param%lclose) then + call netcdf_io_check( nf90_inq_varid(nc%id, nc%radius_varname, nc%radius_varid), "encounter_io_netcdf_open nf90_inq_varid radius_varid" ) + end if + + if (param%lrotation) then + call netcdf_io_check( nf90_inq_varid(nc%id, nc%Ip_varname, nc%Ip_varid), "encounter_io_netcdf_open nf90_inq_varid Ip_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%rot_varname, nc%rot_varid), "encounter_io_netcdf_open nf90_inq_varid rot_varid" ) + end if + + end associate + end select + + return + end subroutine encounter_io_netcdf_open + + module subroutine encounter_io_netcdf_write_frame_snapshot(self, history, param) !! author: David A. Minton !! @@ -176,11 +242,10 @@ module subroutine encounter_io_netcdf_write_frame_snapshot(self, history, param) class is (swiftest_tp) select type (nc => history%nc) class is (encounter_netcdf_parameters) - associate(tslot => param%ioutput) - call netcdf_io_check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "encounter_io_netcdf_write_frame_snapshot nf90_set_fill" ) + associate(tslot => nc%tslot) + call netcdf_io_check( nf90_set_fill(nc%id, NF90_NOFILL, old_mode), "encounter_io_netcdf_write_frame_snapshot nf90_set_fill" ) call netcdf_io_check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "encounter_io_netcdf_write_frame_snapshot nf90_put_var time_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%loop_varid, int(self%iloop,kind=I4B), start=[tslot]), "encounter_io_netcdf_write_frame_snapshot nf90_put_var pl loop_varid" ) npl = pl%nbody do i = 1, npl diff --git a/src/encounter/encounter_module.f90 b/src/encounter/encounter_module.f90 index 32ff67d25..ed7ae4dad 100644 --- a/src/encounter/encounter_module.f90 +++ b/src/encounter/encounter_module.f90 @@ -18,6 +18,7 @@ module encounter public integer(I4B), parameter :: SWEEPDIM = 3 + character(len=*), parameter :: ENCOUNTER_OUTFILE = 'encounters.nc' !! Name of NetCDF output file for encounter information type, abstract :: encounter_list integer(I8B) :: nenc = 0 !! Total number of encounters @@ -60,13 +61,9 @@ module encounter !> NetCDF dimension and variable names for the enounter save object type, extends(netcdf_parameters) :: encounter_netcdf_parameters - character(NAMELEN) :: loop_varname = "loopnum" !! Loop number for encounter - integer(I4B) :: loop_varid !! ID for the recursion level variable - integer(I4B) :: time_dimsize = 0 !! Number of time values in snapshot - integer(I4B) :: name_dimsize = 0 !! Number of potential id values in snapshot - integer(I4B) :: file_number = 1 !! The number to append on the output file contains procedure :: initialize => encounter_io_netcdf_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object + procedure :: open => encounter_io_netcdf_open !! Open an encounter NetCDF file final :: encounter_final_netcdf_parameters !! Finalizer will close the NetCDF file end type encounter_netcdf_parameters @@ -229,6 +226,13 @@ module subroutine encounter_io_netcdf_initialize_output(self, param) class(base_parameters), intent(in) :: param end subroutine encounter_io_netcdf_initialize_output + module subroutine encounter_io_netcdf_open(self, param, readonly) + implicit none + class(encounter_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(base_parameters), intent(in) :: param !! Current run configuration parameters + logical, optional, intent(in) :: readonly !! Logical flag indicating that this should be open read only + end subroutine encounter_io_netcdf_open + module subroutine encounter_io_netcdf_write_frame_snapshot(self, history, param) implicit none class(encounter_snapshot), intent(in) :: self !! Swiftest encounter structure @@ -272,7 +276,6 @@ module subroutine encounter_util_dealloc_list(self) class(encounter_list), intent(inout) :: self !! Swiftest encounter list object end subroutine encounter_util_dealloc_list - module subroutine encounter_util_get_idvalues_snapshot(self, idvals) implicit none class(encounter_snapshot), intent(in) :: self !! Encounter snapshot object diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index d8dbeb0de..c72254ca9 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -136,8 +136,8 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) call self%set_original_scale() - dE = self%Etot(2) - self%Etot(1) - dL = .mag.(self%Ltot(:,2) - self%Ltot(:,1)) + dE = self%te(2) - self%te(1) + dL = .mag.(self%L_total(:,2) - self%L_total(:,1)) write(message,*) dE call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Estimated energy change: " // trim(adjustl(message))) @@ -326,7 +326,7 @@ module subroutine fraggle_generate_rot_vec(collider) ! Arguments class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object ! Internals - real(DP), dimension(NDIM) :: Lbefore, Lafter, Lspin, rotdir + real(DP), dimension(NDIM) :: Lbefore, Lafter, L_spin, rotdir real(DP) :: v_init, v_final, mass_init, mass_final, rotmag integer(I4B) :: i @@ -341,8 +341,8 @@ module subroutine fraggle_generate_rot_vec(collider) Lbefore(:) = mass_init * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (impactors%vb(:,2) - impactors%vb(:,1)) Lafter(:) = mass_final * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (v_final * impactors%bounce_unit(:)) - Lspin(:) = impactors%Lspin(:,1) + (Lbefore(:) - Lafter(:)) - fragments%rot(:,1) = Lspin(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1)) + L_spin(:) = impactors%L_spin(:,1) + (Lbefore(:) - Lafter(:)) + fragments%rot(:,1) = L_spin(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1)) ! Add in some random spin noise. The magnitude will be scaled by the before-after amount and the direction will be random do concurrent(i = 2:nfrag) @@ -488,7 +488,7 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure) if (lsupercat) then ! Put some of the residual angular momentum into velocity shear. Not too much, or we get some weird trajectories call fragments%get_angular_momentum() - Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:)) + Lresidual(:) = fragments%L_budget(:) - (fragments%L_orbit_tot(:) + fragments%L_spin_tot(:)) do concurrent(i = istart:nfrag) vunit(:) = .unit. (Lresidual(:) .cross. fragments%r_unit(:,i)) vshear(:) = vunit(:) * (.mag.Lresidual(:) / ((nfrag-istart+1)*fragments%mass(i) * fragments%rmag(i))) @@ -497,14 +497,14 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure) end if ! Check for any residual angular momentum, and if there is any, put it into spin call fragments%get_angular_momentum() - Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:)) + Lresidual(:) = fragments%L_budget(:) - (fragments%L_orbit_tot(:) + fragments%L_spin_tot(:)) do concurrent(i = 1:nfrag) - fragments%Lspin(:,i) = fragments%Lspin(:,i) + Lresidual(:) / nfrag - fragments%rot(:,i) = fragments%Lspin(:,i) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i)) + fragments%L_spin(:,i) = fragments%L_spin(:,i) + Lresidual(:) / nfrag + fragments%rot(:,i) = fragments%L_spin(:,i) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i)) end do call fragments%get_angular_momentum() - Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:)) + Lresidual(:) = fragments%L_budget(:) - (fragments%L_orbit_tot(:) + fragments%L_spin_tot(:)) end do ! We didn't converge. Try another configuration and see if we get a better result diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index b69a7aea3..caf72fbd3 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -260,11 +260,11 @@ module subroutine fraggle_util_set_natural_scale_factors(self) impactors%Gmass(:) = impactors%Gmass(:) / (collider%dscale**3/collider%tscale**2) impactors%Mcb = impactors%Mcb / collider%mscale impactors%radius(:) = impactors%radius(:) / collider%dscale - impactors%Lspin(:,:) = impactors%Lspin(:,:) / collider%Lscale - impactors%Lorbit(:,:) = impactors%Lorbit(:,:) / collider%Lscale + impactors%L_spin(:,:) = impactors%L_spin(:,:) / collider%Lscale + impactors%L_orbit(:,:) = impactors%L_orbit(:,:) / collider%Lscale do concurrent(i = 1:2) - impactors%rot(:,i) = impactors%Lspin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i)) + impactors%rot(:,i) = impactors%L_spin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i)) end do fragments%mtot = fragments%mtot / collider%mscale @@ -309,10 +309,10 @@ module subroutine fraggle_util_set_original_scale_factors(self) impactors%vb = impactors%vb * collider%vscale impactors%rc = impactors%rc * collider%dscale impactors%vc = impactors%vc * collider%vscale - impactors%Lspin = impactors%Lspin * collider%Lscale - impactors%Lorbit = impactors%Lorbit * collider%Lscale + impactors%L_spin = impactors%L_spin * collider%Lscale + impactors%L_orbit = impactors%L_orbit * collider%Lscale do concurrent(i = 1:2) - impactors%rot(:,i) = impactors%Lspin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i)) + impactors%rot(:,i) = impactors%L_spin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i)) end do fragments%mtot = fragments%mtot * collider%mscale @@ -327,14 +327,14 @@ module subroutine fraggle_util_set_original_scale_factors(self) impactors%Qloss = impactors%Qloss * collider%Escale - collider%Lorbit(:,:) = collider%Lorbit(:,:) * collider%Lscale - collider%Lspin(:,:) = collider%Lspin(:,:) * collider%Lscale - collider%Ltot(:,:) = collider%Ltot(:,:) * collider%Lscale + collider%L_orbit(:,:) = collider%L_orbit(:,:) * collider%Lscale + collider%L_spin(:,:) = collider%L_spin(:,:) * collider%Lscale + collider%L_total(:,:) = collider%L_total(:,:) * collider%Lscale collider%ke_orbit(:) = collider%ke_orbit(:) * collider%Escale collider%ke_spin(:) = collider%ke_spin(:) * collider%Escale collider%pe(:) = collider%pe(:) * collider%Escale collider%be(:) = collider%be(:) * collider%Escale - collider%Etot(:) = collider%Etot(:) * collider%Escale + collider%te(:) = collider%te(:) * collider%Escale collider%mscale = 1.0_DP collider%dscale = 1.0_DP diff --git a/src/helio/helio_drift.f90 b/src/helio/helio_drift.f90 index ad476a2b2..134421d2b 100644 --- a/src/helio/helio_drift.f90 +++ b/src/helio/helio_drift.f90 @@ -28,6 +28,7 @@ module subroutine helio_drift_body(self, nbody_system, param, dt) integer(I4B) :: i !! Loop counter integer(I4B), dimension(:),allocatable :: iflag !! Vectorized error code flag real(DP), dimension(:), allocatable :: mu + character(len=STRMAX) message if (self%nbody == 0) return @@ -40,7 +41,10 @@ module subroutine helio_drift_body(self, nbody_system, param, dt) if (any(iflag(1:n) /= 0)) then where(iflag(1:n) /= 0) self%status(1:n) = DISCARDED_DRIFTERR do i = 1, n - if (iflag(i) /= 0) write(*, *) " Body ", self%id(i), " lost due to error in Danby drift" + if (iflag(i) /= 0) then + write(message, *) " Body ", self%id(i), " lost due to error in Danby drift" + call swiftest_io_log_one_message(COLLISION_LOG_OUT,message) + end if end do end if end associate diff --git a/src/netcdf_io/netcdf_io_implementations.f90 b/src/netcdf_io/netcdf_io_implementations.f90 index 9614e2a32..17f957a98 100644 --- a/src/netcdf_io/netcdf_io_implementations.f90 +++ b/src/netcdf_io/netcdf_io_implementations.f90 @@ -21,7 +21,7 @@ module subroutine netcdf_io_check(status, call_identifier) integer, intent (in) :: status !! The status code returned by a NetCDF function character(len=*), intent(in), optional :: call_identifier !! String that indicates which calling function caused the error for diagnostic purposes - if(status /= nf90_noerr) then + if(status /= NF90_NOERR) then if (present(call_identifier)) write(*,*) "NetCDF error in ",trim(call_identifier) write(*,*) trim(nf90_strerror(status)) call util_exit(FAILURE) diff --git a/src/netcdf_io/netcdf_io_module.f90 b/src/netcdf_io/netcdf_io_module.f90 index 3926d2eea..3adb7b1a5 100644 --- a/src/netcdf_io/netcdf_io_module.f90 +++ b/src/netcdf_io/netcdf_io_module.f90 @@ -24,8 +24,8 @@ module netcdf_io logical :: lfile_is_open = .false. !! Flag indicating that the linked file is currently open integer(I4B) :: out_type !! output type (will be assigned either NF90_DOUBLE or NF90_FLOAT, depending on the user parameter) integer(I4B) :: id !! ID for the output file - integer(I4B) :: name_chunk !! Chunk size for the id dimension variables - integer(I4B) :: time_chunk !! Chunk size for the time dimension variables + integer(I4B) :: max_tslot = 0 !! Records the last index value of time in the NetCDF file + integer(I4B) :: tslot = 1 !! The current time slot that gets passed to the NetCDF reader/writer ! Dimension ids and variable names character(NAMELEN) :: str_dimname = "string32" !! name of the character string dimension @@ -44,6 +44,8 @@ module netcdf_io ! Non-dimension ids and variable names character(NAMELEN) :: id_varname = "id" !! name of the particle id variable integer(I4B) :: id_varid !! ID for the id variable + character(NAMELEN) :: status_varname = "status" !! name of the particle status variable + integer(I4B) :: status_varid !! ID for the status variable character(NAMELEN) :: ptype_varname = "particle_type" !! name of the particle type variable integer(I4B) :: ptype_varid !! ID for the particle type variable character(NAMELEN) :: npl_varname = "npl" !! name of the number of active massive bodies variable @@ -78,7 +80,7 @@ module netcdf_io integer(I4B) :: vh_varid !! ID for the heliocentric velocity vector variable character(NAMELEN) :: gr_pseudo_vh_varname = "gr_pseudo_vh" !! name of the heliocentric pseudovelocity vector variable (used in GR only) integer(I4B) :: gr_pseudo_vh_varid !! ID for the heliocentric pseudovelocity vector variable (used in GR) - character(NAMELEN) :: gmass_varname = "Gmass" !! name of the mass variable + character(NAMELEN) :: Gmass_varname = "Gmass" !! name of the mass variable integer(I4B) :: Gmass_varid !! ID for the mass variable character(NAMELEN) :: rhill_varname = "rhill" !! name of the hill radius variable integer(I4B) :: rhill_varid !! ID for the hill radius variable @@ -104,16 +106,16 @@ module netcdf_io integer(I4B) :: PE_varid !! ID for the system potential energy variable character(NAMELEN) :: be_varname = "BE" !! name of the system binding energy variable integer(I4B) :: BE_varid !! ID for the system binding energy variable - character(NAMELEN) :: L_orb_varname = "L_orb" !! name of the orbital angular momentum vector variable - integer(I4B) :: L_orb_varid !! ID for the system orbital angular momentum vector variable - character(NAMELEN) :: Lspin_varname = "Lspin" !! name of the spin angular momentum vector variable - integer(I4B) :: Lspin_varid !! ID for the system spin angular momentum vector variable + character(NAMELEN) :: L_orbit_varname = "L_orbit" !! name of the orbital angular momentum vector variable + integer(I4B) :: L_orbit_varid !! ID for the system orbital angular momentum vector variable + character(NAMELEN) :: L_spin_varname = "L_spin" !! name of the spin angular momentum vector variable + integer(I4B) :: L_spin_varid !! ID for the system spin angular momentum vector variable character(NAMELEN) :: L_escape_varname = "L_escape" !! name of the escaped angular momentum vector variable integer(I4B) :: L_escape_varid !! ID for the escaped angular momentum vector variable - character(NAMELEN) :: Ecollisions_varname = "Ecollisions" !! name of the escaped angular momentum y variable - integer(I4B) :: Ecollisions_varid !! ID for the energy lost in collisions variable - character(NAMELEN) :: Euntracked_varname = "Euntracked" !! name of the energy that is untracked due to loss (untracked potential energy due to mergers and body energy for escaped bodies) - integer(I4B) :: Euntracked_varid !! ID for the energy that is untracked due to loss (untracked potential energy due to mergers and body energy for escaped bodies) + character(NAMELEN) :: E_collisions_varname = "E_collisions" !! name of the escaped angular momentum y variable + integer(I4B) :: E_collisions_varid !! ID for the energy lost in collisions variable + character(NAMELEN) :: E_untracked_varname = "E_untracked" !! name of the energy that is untracked due to loss (untracked potential energy due to mergers and body energy for escaped bodies) + integer(I4B) :: E_untracked_varid !! ID for the energy that is untracked due to loss (untracked potential energy due to mergers and body energy for escaped bodies) character(NAMELEN) :: GMescape_varname = "GMescape" !! name of the G*Mass of bodies that escape the system integer(I4B) :: GMescape_varid !! ID for the G*Mass of bodies that escape the system character(NAMELEN) :: origin_type_varname = "origin_type" !! name of the origin type variable (Initial Conditions, Disruption, etc.) diff --git a/src/rmvs/rmvs_discard.f90 b/src/rmvs/rmvs_discard.f90 index 5a8fd94f5..7e0cb9905 100644 --- a/src/rmvs/rmvs_discard.f90 +++ b/src/rmvs/rmvs_discard.f90 @@ -25,7 +25,7 @@ module subroutine rmvs_discard_tp(self, nbody_system, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i - character(len=STRMAX) :: timestr, idstri, idstrj + character(len=STRMAX) :: timestr, idstri, idstrj, message if (self%nbody == 0) return @@ -38,9 +38,10 @@ module subroutine rmvs_discard_tp(self, nbody_system, param) write(idstri, *) tp%id(i) write(idstrj, *) pl%id(iplperP) write(timestr, *) t - write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) & + write(message, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) & // ") q with respect to massive body " // trim(adjustl(pl%info(iplperP)%name)) & // " (" // trim(adjustl(idstrj)) // ") is too small at t = " // trim(adjustl(timestr)) + call swiftest_io_log_one_message(COLLISION_LOG_OUT,message) tp%ldiscard(i) = .true. tp%lmask(i) = .false. call tp%info(i)%set_value(status="DISCARDED_PLQ", discard_time=t, discard_rh=tp%rh(:,i), & diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 1bc313336..19b8779f6 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -93,6 +93,7 @@ subroutine rmvs_interp_out(cb, pl, dt) real(DP), dimension(:,:), allocatable :: xtmp, vtmp real(DP), dimension(:), allocatable :: GMcb, dto integer(I4B), dimension(:), allocatable :: iflag + character(len=STRMAX) :: message dntenc = real(NTENC, kind=DP) associate (npl => pl%nbody) @@ -112,11 +113,12 @@ subroutine rmvs_interp_out(cb, pl, dt) if (any(iflag(1:npl) /= 0)) then do i = 1, npl if (iflag(i) /= 0) then - write(*, *) " Planet ", pl%id(i), " is lost!!!!!!!!!!" - write(*, *) GMcb(i), dto(i) - write(*, *) xtmp(:,i) - write(*, *) vtmp(:,i) - write(*, *) " STOPPING " + write(message, *) " Planet ", pl%id(i), " is lost!!!!!!!!!!",new_line('a'), & + GMcb(i), dto(i),new_line('a'), & + xtmp(:,i),new_line('a'), & + vtmp(:,i),new_line('a'), & + " STOPPING " + call swiftest_io_log_one_message(COLLISION_LOG_OUT,message) call util_exit(FAILURE) end if end do @@ -134,11 +136,12 @@ subroutine rmvs_interp_out(cb, pl, dt) if (any(iflag(1:npl) /= 0)) then do i = 1, npl if (iflag(i) /= 0) then - write(*, *) " Planet ", pl%id(i), " is lost!!!!!!!!!!" - write(*, *) GMcb(i), -dto(i) - write(*, *) xtmp(:,i) - write(*, *) vtmp(:,i) - write(*, *) " STOPPING " + write(message, *) " Planet ", pl%id(i), " is lost!!!!!!!!!!",new_line('a'), & + GMcb(i), -dto(i), new_line('a'), & + xtmp(:,i), new_line('a'), & + vtmp(:,i), new_line('a'), & + " STOPPING " + call swiftest_io_log_one_message(COLLISION_LOG_OUT,message) call util_exit(FAILURE) end if end do @@ -237,6 +240,7 @@ subroutine rmvs_interp_in(cb, pl, nbody_system, param, dt, outer_index) real(DP), dimension(:,:), allocatable :: xtmp, vtmp, rh_original, ah_original real(DP), dimension(:), allocatable :: GMcb, dti integer(I4B), dimension(:), allocatable :: iflag + character(len=STRMAX) :: message associate (npl => nbody_system%pl%nbody) dntphenc = real(NTPHENC, kind=DP) @@ -279,11 +283,12 @@ subroutine rmvs_interp_in(cb, pl, nbody_system, param, dt, outer_index) if (any(iflag(1:npl) /= 0)) then do i = 1, npl if (iflag(i) /=0) then - write(*, *) " Planet ", pl%id(i), " is lost!!!!!!!!!!" - write(*, *) GMcb(i), dti(i) - write(*, *) xtmp(:,i) - write(*, *) vtmp(:,i) - write(*, *) " STOPPING " + write(message, *) " Planet ", pl%id(i), " is lost!!!!!!!!!!", new_line('a'), & + GMcb(i), dti(i), new_line('a'), & + xtmp(:,i), new_line('a'), & + vtmp(:,i), new_line('a'), & + " STOPPING " + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) call util_exit(failure) end if end do diff --git a/src/swiftest/swiftest_discard.f90 b/src/swiftest/swiftest_discard.f90 index eafb53452..ad5be95ae 100644 --- a/src/swiftest/swiftest_discard.f90 +++ b/src/swiftest/swiftest_discard.f90 @@ -200,7 +200,7 @@ subroutine swiftest_discard_peri_tp(tp, nbody_system, param) integer(I4B) :: i, j, ih real(DP) :: r2 real(DP), dimension(NDIM) :: dx - character(len=STRMAX) :: idstr, timestr + character(len=STRMAX) :: idstr, timestr, message associate(cb => nbody_system%cb, ntp => tp%nbody, pl => nbody_system%pl, npl => nbody_system%pl%nbody, t => nbody_system%t) call tp%get_peri(nbody_system, param) @@ -220,8 +220,10 @@ subroutine swiftest_discard_peri_tp(tp, nbody_system, param) tp%status(i) = DISCARDED_PERI write(idstr, *) tp%id(i) write(timestr, *) nbody_system%t - write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & + write(message, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & " perihelion distance too small at t = " // trim(adjustl(timestr)) + + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) tp%ldiscard(i) = .true. call tp%info(i)%set_value(status="DISCARDED_PERI", discard_time=nbody_system%t, discard_rh=tp%rh(:,i), & discard_vh=tp%vh(:,i), discard_body_id=pl%id(j)) diff --git a/src/swiftest/swiftest_drift.f90 b/src/swiftest/swiftest_drift.f90 index 89f8afa16..b7811f88c 100644 --- a/src/swiftest/swiftest_drift.f90 +++ b/src/swiftest/swiftest_drift.f90 @@ -34,6 +34,7 @@ module subroutine swiftest_drift_body(self, nbody_system, param, dt) ! Internals integer(I4B) :: i integer(I4B), dimension(:), allocatable :: iflag + character(len=STRMAX) :: message associate(n => self%nbody) allocate(iflag(n)) @@ -42,7 +43,10 @@ module subroutine swiftest_drift_body(self, nbody_system, param, dt) if (any(iflag(1:n) /= 0)) then where(iflag(1:n) /= 0) self%status(1:n) = DISCARDED_DRIFTERR do i = 1, n - if (iflag(i) /= 0) write(*, *) " Body ", self%id(i), " lost due to error in Danby drift" + if (iflag(i) /= 0) then + write(message, *) " Body ", self%id(i), " lost due to error in Danby drift" + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + end if end do end if end associate diff --git a/src/swiftest/swiftest_driver.f90 b/src/swiftest/swiftest_driver.f90 index 4d8279326..e878bdd72 100644 --- a/src/swiftest/swiftest_driver.f90 +++ b/src/swiftest/swiftest_driver.f90 @@ -46,7 +46,7 @@ program swiftest_driver !> Read in the user-defined parameters file and the initial conditions of the nbody_system allocate(swiftest_parameters :: param) param%integrator = trim(adjustl(integrator)) - call param%set_display(display_style) + param%display_style = trim(adjustl(display_style)) call param%read_in(param_file_name) associate(t0 => param%t0, & @@ -56,16 +56,13 @@ program swiftest_driver iloop => param%iloop, & istep_out => param%istep_out, & dump_cadence => param%dump_cadence, & - ioutput => param%ioutput, & display_style => param%display_style, & display_unit => param%display_unit) - ! Set up loop and output cadence variables nloops = ceiling((tstop - t0) / dt, kind=I8B) istart = ceiling((tstart - t0) / dt + 1.0_DP, kind=I8B) iloop = istart - 1 - ioutput = max(int(istart / istep_out, kind=I4B),1) ! Set up nbody_system storage for intermittent file dumps if (dump_cadence == 0) dump_cadence = ceiling(nloops / (1.0_DP * istep_out), kind=I8B) diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index d88402be2..4fe297197 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -41,11 +41,11 @@ module subroutine swiftest_io_compact_output(self, param, timer) formatted_output = formatted_output // fmt("NPLM",pl%nplm) end select if (param%lenergy) then - formatted_output = formatted_output // fmt("LTOTERR",self%Ltot_error) // fmt("ETOTERR",self%Etot_error) // fmt("MTOTERR",self%Mtot_error) & - // fmt("KEOERR",self%ke_orbit_error) // fmt("PEERR",self%pe_error) // fmt("EORBERR",self%Eorbit_error) & - // fmt("EUNTRERR",self%Euntracked_error) // fmt("LESCERR",self%Lescape_error) // fmt("MESCERR",self%Mescape_error) + formatted_output = formatted_output // fmt("LTOTERR",self%L_total_error) // fmt("ETOTERR",self%te_error) // fmt("MTOTERR",self%Mtot_error) & + // fmt("KEOERR",self%ke_orbit_error) // fmt("PEERR",self%pe_error) // fmt("EORBERR",self%E_orbit_error) & + // fmt("EUNTRERR",self%E_untracked_error) // fmt("LESCERR",self%L_escape_error) // fmt("MESCERR",self%Mescape_error) if (param%lclose) formatted_output = formatted_output // fmt("ECOLLERR",self%Ecoll_error) - if (param%lrotation) formatted_output = formatted_output // fmt("KESPINERR",self%ke_spin_error) // fmt("LSPINERR",self%Lspin_error) + if (param%lrotation) formatted_output = formatted_output // fmt("KESPINERR",self%ke_spin_error) // fmt("LSPINERR",self%L_spin_error) end if if (.not. timer%main_is_started) then ! This is the start of a new run @@ -120,14 +120,14 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal) class(swiftest_parameters), intent(inout) :: param !! Input colleciton of user-defined parameters logical, intent(in) :: lterminal !! Indicates whether to output information to the terminal screen ! Internals - real(DP), dimension(NDIM) :: Ltot_now, Lorbit_now, Lspin_now - real(DP) :: ke_orbit_now, ke_spin_now, pe_now, Eorbit_now, be_now + real(DP), dimension(NDIM) :: L_total_now, L_orbit_now, L_spin_now + real(DP) :: ke_orbit_now, ke_spin_now, pe_now, E_orbit_now, be_now real(DP) :: GMtot_now character(len=STRMAX) :: errmsg integer(I4B), parameter :: EGYIU = 72 character(len=*), parameter :: EGYTERMFMT = '(" DL/L0 = ", ES12.5 & - "; DEcollisions/|E0| = ", ES12.5, & - "; D(Eorbit+Ecollisions)/|E0| = ", ES12.5, & + "; DE_collisions/|E0| = ", ES12.5, & + "; D(E_orbit+E_collisions)/|E0| = ", ES12.5, & "; DM/M0 = ", ES12.5)' associate(nbody_system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody, display_unit => param%display_unit, nc => param%system_history%nc) @@ -140,10 +140,10 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal) ke_spin_now = nbody_system%ke_spin pe_now = nbody_system%pe be_now = nbody_system%be - Lorbit_now(:) = nbody_system%Lorbit(:) - Lspin_now(:) = nbody_system%Lspin(:) - Eorbit_now = ke_orbit_now + ke_spin_now + pe_now + be_now - Ltot_now(:) = nbody_system%Ltot(:) + nbody_system%Lescape(:) + L_orbit_now(:) = nbody_system%L_orbit(:) + L_spin_now(:) = nbody_system%L_spin(:) + E_orbit_now = ke_orbit_now + ke_spin_now + pe_now + be_now + L_total_now(:) = nbody_system%L_total(:) + nbody_system%L_escape(:) GMtot_now = nbody_system%GMtot + nbody_system%GMescape if (param%lfirstenergy) then @@ -151,35 +151,34 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal) nbody_system%ke_spin_orig = ke_spin_now nbody_system%pe_orig = pe_now nbody_system%be_orig = be_now - nbody_system%Eorbit_orig = Eorbit_now + nbody_system%E_orbit_orig = E_orbit_now nbody_system%GMtot_orig = GMtot_now - nbody_system%Lorbit_orig(:) = Lorbit_now(:) - nbody_system%Lspin_orig(:) = Lspin_now(:) - nbody_system%Ltot_orig(:) = Ltot_now(:) + nbody_system%L_orbit_orig(:) = L_orbit_now(:) + nbody_system%L_spin_orig(:) = L_spin_now(:) + nbody_system%L_total_orig(:) = L_total_now(:) param%lfirstenergy = .false. end if if (.not.param%lfirstenergy) then - nbody_system%ke_orbit_error = (ke_orbit_now - nbody_system%ke_orbit_orig) / abs(nbody_system%Eorbit_orig) - nbody_system%ke_spin_error = (ke_spin_now - nbody_system%ke_spin_orig) / abs(nbody_system%Eorbit_orig) - nbody_system%pe_error = (pe_now - nbody_system%pe_orig) / abs(nbody_system%Eorbit_orig) - nbody_system%be_error = (be_now - nbody_system%be_orig) / abs(nbody_system%Eorbit_orig) - nbody_system%Eorbit_error = (Eorbit_now - nbody_system%Eorbit_orig) / abs(nbody_system%Eorbit_orig) - nbody_system%Ecoll_error = nbody_system%Ecollisions / abs(nbody_system%Eorbit_orig) - nbody_system%Euntracked_error = nbody_system%Euntracked / abs(nbody_system%Eorbit_orig) - nbody_system%Etot_error = (Eorbit_now - nbody_system%Ecollisions - nbody_system%Eorbit_orig - nbody_system%Euntracked) / abs(nbody_system%Eorbit_orig) - - nbody_system%Lorbit_error = norm2(Lorbit_now(:) - nbody_system%Lorbit_orig(:)) / norm2(nbody_system%Ltot_orig(:)) - nbody_system%Lspin_error = norm2(Lspin_now(:) - nbody_system%Lspin_orig(:)) / norm2(nbody_system%Ltot_orig(:)) - nbody_system%Lescape_error = norm2(nbody_system%Lescape(:)) / norm2(nbody_system%Ltot_orig(:)) - nbody_system%Ltot_error = norm2(Ltot_now(:) - nbody_system%Ltot_orig(:)) / norm2(nbody_system%Ltot_orig(:)) + nbody_system%ke_orbit_error = (ke_orbit_now - nbody_system%ke_orbit_orig) / abs(nbody_system%E_orbit_orig) + nbody_system%ke_spin_error = (ke_spin_now - nbody_system%ke_spin_orig) / abs(nbody_system%E_orbit_orig) + nbody_system%pe_error = (pe_now - nbody_system%pe_orig) / abs(nbody_system%E_orbit_orig) + nbody_system%be_error = (be_now - nbody_system%be_orig) / abs(nbody_system%E_orbit_orig) + nbody_system%E_orbit_error = (E_orbit_now - nbody_system%E_orbit_orig) / abs(nbody_system%E_orbit_orig) + nbody_system%Ecoll_error = nbody_system%E_collisions / abs(nbody_system%E_orbit_orig) + nbody_system%E_untracked_error = nbody_system%E_untracked / abs(nbody_system%E_orbit_orig) + nbody_system%te_error = (E_orbit_now - nbody_system%E_collisions - nbody_system%E_orbit_orig - nbody_system%E_untracked) / abs(nbody_system%E_orbit_orig) + + nbody_system%L_orbit_error = norm2(L_orbit_now(:) - nbody_system%L_orbit_orig(:)) / norm2(nbody_system%L_total_orig(:)) + nbody_system%L_spin_error = norm2(L_spin_now(:) - nbody_system%L_spin_orig(:)) / norm2(nbody_system%L_total_orig(:)) + nbody_system%L_escape_error = norm2(nbody_system%L_escape(:)) / norm2(nbody_system%L_total_orig(:)) + nbody_system%L_total_error = norm2(L_total_now(:) - nbody_system%L_total_orig(:)) / norm2(nbody_system%L_total_orig(:)) nbody_system%Mescape_error = nbody_system%GMescape / nbody_system%GMtot_orig nbody_system%Mtot_error = (GMtot_now - nbody_system%GMtot_orig) / nbody_system%GMtot_orig - if (lterminal) write(display_unit, EGYTERMFMT) nbody_system%Ltot_error, nbody_system%Ecoll_error, nbody_system%Etot_error,nbody_system%Mtot_error + if (lterminal) write(display_unit, EGYTERMFMT) nbody_system%L_total_error, nbody_system%Ecoll_error, nbody_system%te_error,nbody_system%Mtot_error if (abs(nbody_system%Mtot_error) > 100 * epsilon(nbody_system%Mtot_error)) then write(*,*) "Severe error! Mass not conserved! Halting!" ! Save the frame of data to the bin file in the slot just after the present one for diagnostics - param%ioutput = param%ioutput + 1 call self%write_frame(nc, param) call nc%close() call util_exit(FAILURE) @@ -275,18 +274,27 @@ module subroutine swiftest_io_dump_storage(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i - integer(I8B) :: iloop_start + real(DP), dimension(:), allocatable :: tvals if (self%iframe == 0) return - iloop_start = max(param%iloop / param%istep_out - int(param%dump_cadence, kind=I8B) + 1_I8B,0_I8B) call self%make_index_map() associate(nc => param%system_history%nc) call nc%open(param) + ! Get the current time values in the file if this is an old file + if (nc%max_tslot > 0) then + allocate(tvals(nc%max_tslot)) + call netcdf_io_check( nf90_get_var(nc%id, nc%time_varid, tvals(:), start=[1]), "swiftest_io_dump_storage time_varid" ) + else + allocate(tvals(1)) + tvals(1) = -huge(1.0_DP) + end if do i = 1, self%iframe if (allocated(self%frame(i)%item)) then - param%ioutput = iloop_start + self%tmap(i) select type(nbody_system => self%frame(i)%item) class is (swiftest_nbody_system) + nc%tslot = findloc(tvals, nbody_system%t, dim=1) + if (nc%tslot == 0) nc%tslot = nc%max_tslot + 1 + nc%max_tslot = max(nc%max_tslot, nc%tslot) call nbody_system%write_frame(param) end select deallocate(self%frame(i)%item) @@ -523,16 +531,16 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, param) call netcdf_io_check( nf90_get_var(nc%id, nc%BE_varid, rtemp, start=[1], count=[1]), "netcdf_io_get_t0_values_system BE_varid" ) BE_orig = rtemp(1) - call netcdf_io_check( nf90_get_var(nc%id, nc%Ecollisions_varid, self%Ecollisions, start=[1]), "netcdf_io_get_t0_values_system Ecollisions_varid" ) - call netcdf_io_check( nf90_get_var(nc%id, nc%Euntracked_varid, self%Euntracked, start=[1]), "netcdf_io_get_t0_values_system Euntracked_varid" ) + call netcdf_io_check( nf90_get_var(nc%id, nc%E_collisions_varid, self%E_collisions, start=[1]), "netcdf_io_get_t0_values_system E_collisions_varid" ) + call netcdf_io_check( nf90_get_var(nc%id, nc%E_untracked_varid, self%E_untracked, start=[1]), "netcdf_io_get_t0_values_system E_untracked_varid" ) - self%Eorbit_orig = KE_orb_orig + KE_spin_orig + PE_orig + BE_orig + self%Ecollisions + self%Euntracked + self%E_orbit_orig = KE_orb_orig + KE_spin_orig + PE_orig + BE_orig + self%E_collisions + self%E_untracked - call netcdf_io_check( nf90_get_var(nc%id, nc%L_orb_varid, self%Lorbit_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_t0_values_system L_orb_varid" ) - call netcdf_io_check( nf90_get_var(nc%id, nc%Lspin_varid, self%Lspin_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_t0_values_system Lspin_varid" ) - call netcdf_io_check( nf90_get_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_t0_values_system L_escape_varid" ) + call netcdf_io_check( nf90_get_var(nc%id, nc%L_orbit_varid, self%L_orbit_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_t0_values_system L_orbit_varid" ) + call netcdf_io_check( nf90_get_var(nc%id, nc%L_spin_varid, self%L_spin_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_t0_values_system L_spin_varid" ) + call netcdf_io_check( nf90_get_var(nc%id, nc%L_escape_varid, self%L_escape(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_t0_values_system L_escape_varid" ) - self%Ltot_orig(:) = self%Lorbit_orig(:) + self%Lspin_orig(:) + self%Lescape(:) + self%L_total_orig(:) = self%L_orbit_orig(:) + self%L_spin_orig(:) + self%L_escape(:) call netcdf_io_check( nf90_get_var(nc%id, nc%Gmass_varid, vals, start=[1,1], count=[idmax,1]), "netcdf_io_get_t0_values_system Gmass_varid" ) call netcdf_io_check( nf90_get_var(nc%id, nc%GMescape_varid, self%GMescape, start=[1]), "netcdf_io_get_t0_values_system GMescape_varid" ) @@ -621,6 +629,7 @@ module subroutine swiftest_io_netcdf_initialize_output(self, param) ! Variables call netcdf_io_check( nf90_def_var(nc%id, nc%id_varname, NF90_INT, nc%name_dimid, nc%id_varid), "netcdf_io_initialize_output nf90_def_var id_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%status_varname, NF90_INT, [nc%name_dimid, nc%time_dimid], nc%status_varid), "netcdf_io_initialize_output nf90_def_var status_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%npl_varname, NF90_INT, nc%time_dimid, nc%npl_varid), "netcdf_io_initialize_output nf90_def_var npl_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%ntp_varname, NF90_INT, nc%time_dimid, nc%ntp_varid), "netcdf_io_initialize_output nf90_def_var ntp_varid" ) if (param%lmtiny_pl) call netcdf_io_check( nf90_def_var(nc%id, nc%nplm_varname, NF90_INT, nc%time_dimid, nc%nplm_varid), "netcdf_io_initialize_output nf90_def_var nplm_varid" ) @@ -652,7 +661,7 @@ module subroutine swiftest_io_netcdf_initialize_output(self, param) call netcdf_io_check( nf90_def_var(nc%id, nc%cape_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%cape_varid), "netcdf_io_initialize_output nf90_def_var cape_varid" ) end if - call netcdf_io_check( nf90_def_var(nc%id, nc%gmass_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%Gmass_varid), "netcdf_io_initialize_output nf90_def_var Gmass_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%Gmass_varid), "netcdf_io_initialize_output nf90_def_var Gmass_varid" ) if (param%lrhill_present) then call netcdf_io_check( nf90_def_var(nc%id, nc%rhill_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%rhill_varid), "netcdf_io_initialize_output nf90_def_var rhill_varid" ) @@ -689,11 +698,11 @@ module subroutine swiftest_io_netcdf_initialize_output(self, param) call netcdf_io_check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type, nc%time_dimid, nc%KE_spin_varid), "netcdf_io_initialize_output nf90_def_var KE_spin_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type, nc%time_dimid, nc%PE_varid), "netcdf_io_initialize_output nf90_def_var PE_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%be_varname, nc%out_type, nc%time_dimid, nc%BE_varid), "netcdf_io_initialize_output nf90_def_var PE_varid" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%L_orb_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%L_orb_varid), "netcdf_io_initialize_output nf90_def_var L_orb_varid" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%Lspin_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%Lspin_varid), "netcdf_io_initialize_output nf90_def_var Lspin_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%L_orbit_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%L_orbit_varid), "netcdf_io_initialize_output nf90_def_var L_orbit_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%L_spin_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%L_spin_varid), "netcdf_io_initialize_output nf90_def_var L_spin_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%L_escape_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%L_escape_varid), "netcdf_io_initialize_output nf90_def_var L_escape_varid" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%Ecollisions_varname, nc%out_type, nc%time_dimid, nc%Ecollisions_varid), "netcdf_io_initialize_output nf90_def_var Ecollisions_varid" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%Euntracked_varname, nc%out_type, nc%time_dimid, nc%Euntracked_varid), "netcdf_io_initialize_output nf90_def_var Euntracked_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%E_collisions_varname, nc%out_type, nc%time_dimid, nc%E_collisions_varid), "netcdf_io_initialize_output nf90_def_var E_collisions_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%E_untracked_varname, nc%out_type, nc%time_dimid, nc%E_untracked_varid), "netcdf_io_initialize_output nf90_def_var E_untracked_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%GMescape_varname, nc%out_type, nc%time_dimid, nc%GMescape_varid), "netcdf_io_initialize_output nf90_def_var GMescape_varid" ) end if @@ -707,7 +716,12 @@ module subroutine swiftest_io_netcdf_initialize_output(self, param) call netcdf_io_check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "netcdf_io_initialize_output nf90_inquire_variable" ) select case(vartype) case(NF90_INT) - call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "netcdf_io_initialize_output nf90_def_var_fill NF90_INT" ) + if (varid == nc%status_varid) then + ! Be sure the status variable fill value is the INACTIVE symbolic value + call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, INACTIVE), "netcdf_io_netcdf_initialize_output nf90_def_var_fill status variable" ) + else + call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "netcdf_io_initialize_output nf90_def_var_fill NF90_INT" ) + end if case(NF90_FLOAT) call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "netcdf_io_initialize_output nf90_def_var_fill NF90_FLOAT" ) case(NF90_DOUBLE) @@ -760,33 +774,34 @@ module subroutine swiftest_io_netcdf_open(self, param, readonly) associate(nc => self) - write(errmsg,*) "netcdf_io_open nf90_open ",trim(adjustl(nc%file_name)) + write(errmsg,*) "swiftest_io_netcdf_open nf90_open ",trim(adjustl(nc%file_name)) call netcdf_io_check( nf90_open(nc%file_name, mode, nc%id), errmsg) self%lfile_is_open = .true. ! Dimensions - call netcdf_io_check( nf90_inq_dimid(nc%id, nc%time_dimname, nc%time_dimid), "netcdf_io_open nf90_inq_dimid time_dimid" ) - call netcdf_io_check( nf90_inq_dimid(nc%id, nc%space_dimname, nc%space_dimid), "netcdf_io_open nf90_inq_dimid space_dimid" ) - call netcdf_io_check( nf90_inq_dimid(nc%id, nc%name_dimname, nc%name_dimid), "netcdf_io_open nf90_inq_dimid name_dimid" ) - call netcdf_io_check( nf90_inq_dimid(nc%id, nc%str_dimname, nc%str_dimid), "netcdf_io_open nf90_inq_dimid str_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%time_dimname, nc%time_dimid), "swiftest_io_netcdf_open nf90_inq_dimid time_dimid" ) + call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%time_dimid, nc%time_dimname, len=nc%max_tslot), "swiftest_io_netcdf_open nf90_inquire_dimension max_tslot" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%space_dimname, nc%space_dimid), "swiftest_io_netcdf_open nf90_inq_dimid space_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%name_dimname, nc%name_dimid), "swiftest_io_netcdf_open nf90_inq_dimid name_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%str_dimname, nc%str_dimid), "swiftest_io_netcdf_open nf90_inq_dimid str_dimid" ) ! Dimension coordinates - call netcdf_io_check( nf90_inq_varid(nc%id, nc%time_dimname, nc%time_varid), "netcdf_io_open nf90_inq_varid time_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%space_dimname, nc%space_varid), "netcdf_io_open nf90_inq_varid space_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%name_dimname, nc%name_varid), "netcdf_io_open nf90_inq_varid name_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%time_dimname, nc%time_varid), "swiftest_io_netcdf_open nf90_inq_varid time_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%space_dimname, nc%space_varid), "swiftest_io_netcdf_open nf90_inq_varid space_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%name_dimname, nc%name_varid), "swiftest_io_netcdf_open nf90_inq_varid name_varid" ) ! Required Variables - call netcdf_io_check( nf90_inq_varid(nc%id, nc%id_varname, nc%id_varid), "netcdf_io_open nf90_inq_varid name_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%gmass_varname, nc%Gmass_varid), "netcdf_io_open nf90_inq_varid Gmass_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%id_varname, nc%id_varid), "swiftest_io_netcdf_open nf90_inq_varid name_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%Gmass_varname, nc%Gmass_varid), "swiftest_io_netcdf_open nf90_inq_varid Gmass_varid" ) if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then - call netcdf_io_check( nf90_inq_varid(nc%id, nc%rh_varname, nc%rh_varid), "netcdf_io_open nf90_inq_varid rh_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%vh_varname, nc%vh_varid), "netcdf_io_open nf90_inq_varid vh_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%rh_varname, nc%rh_varid), "swiftest_io_netcdf_open nf90_inq_varid rh_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%vh_varname, nc%vh_varid), "swiftest_io_netcdf_open nf90_inq_varid vh_varid" ) if (param%lgr) then !! check if pseudovelocity vectors exist in this file. If they are, set the correct flag so we know whe should not do the conversion. status = nf90_inq_varid(nc%id, nc%gr_pseudo_vh_varname, nc%gr_pseudo_vh_varid) - nc%lpseudo_vel_exists = (status == nf90_noerr) + nc%lpseudo_vel_exists = (status == NF90_NOERR) if (param%lrestart .and. .not.nc%lpseudo_vel_exists) then write(*,*) "Warning! Pseudovelocity not found in input file for GR enabled run. If this is a restarted run, bit-identical trajectories are not guarunteed!" end if @@ -795,36 +810,37 @@ module subroutine swiftest_io_netcdf_open(self, param, readonly) end if if ((param%out_form == "EL") .or. (param%out_form == "XVEL")) then - call netcdf_io_check( nf90_inq_varid(nc%id, nc%a_varname, nc%a_varid), "netcdf_io_open nf90_inq_varid a_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%e_varname, nc%e_varid), "netcdf_io_open nf90_inq_varid e_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%inc_varname, nc%inc_varid), "netcdf_io_open nf90_inq_varid inc_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%capom_varname, nc%capom_varid), "netcdf_io_open nf90_inq_varid capom_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%omega_varname, nc%omega_varid), "netcdf_io_open nf90_inq_varid omega_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%capm_varname, nc%capm_varid), "netcdf_io_open nf90_inq_varid capm_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%a_varname, nc%a_varid), "swiftest_io_netcdf_open nf90_inq_varid a_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%e_varname, nc%e_varid), "swiftest_io_netcdf_open nf90_inq_varid e_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%inc_varname, nc%inc_varid), "swiftest_io_netcdf_open nf90_inq_varid inc_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%capom_varname, nc%capom_varid), "swiftest_io_netcdf_open nf90_inq_varid capom_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%omega_varname, nc%omega_varid), "swiftest_io_netcdf_open nf90_inq_varid omega_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%capm_varname, nc%capm_varid), "swiftest_io_netcdf_open nf90_inq_varid capm_varid" ) end if if (param%lclose) then - call netcdf_io_check( nf90_inq_varid(nc%id, nc%radius_varname, nc%radius_varid), "netcdf_io_open nf90_inq_varid radius_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%radius_varname, nc%radius_varid), "swiftest_io_netcdf_open nf90_inq_varid radius_varid" ) end if if (param%lrotation) then - call netcdf_io_check( nf90_inq_varid(nc%id, nc%Ip_varname, nc%Ip_varid), "netcdf_io_open nf90_inq_varid Ip_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%rot_varname, nc%rot_varid), "netcdf_io_open nf90_inq_varid rot_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%Ip_varname, nc%Ip_varid), "swiftest_io_netcdf_open nf90_inq_varid Ip_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%rot_varname, nc%rot_varid), "swiftest_io_netcdf_open nf90_inq_varid rot_varid" ) end if ! if (param%ltides) then - ! call netcdf_io_check( nf90_inq_varid(nc%id, nc%k2_varname, nc%k2_varid), "netcdf_io_open nf90_inq_varid k2_varid" ) - ! call netcdf_io_check( nf90_inq_varid(nc%id, nc%q_varname, nc%Q_varid), "netcdf_io_open nf90_inq_varid Q_varid" ) + ! call netcdf_io_check( nf90_inq_varid(nc%id, nc%k2_varname, nc%k2_varid), "swiftest_io_netcdf_open nf90_inq_varid k2_varid" ) + ! call netcdf_io_check( nf90_inq_varid(nc%id, nc%q_varname, nc%Q_varid), "swiftest_io_netcdf_open nf90_inq_varid Q_varid" ) ! end if ! Optional Variables if (param%lrhill_present) then status = nf90_inq_varid(nc%id, nc%rhill_varname, nc%rhill_varid) - if (status /= nf90_noerr) write(*,*) "Warning! RHILL variable not set in input file. Calculating." + if (status /= NF90_NOERR) write(*,*) "Warning! RHILL variable not set in input file. Calculating." end if ! Optional variables The User Doesn't Need to Know About status = nf90_inq_varid(nc%id, nc%npl_varname, nc%npl_varid) + status = nf90_inq_varid(nc%id, nc%status_varname, nc%status_varid) status = nf90_inq_varid(nc%id, nc%ntp_varname, nc%ntp_varid) status = nf90_inq_varid(nc%id, nc%j2rp2_varname, nc%j2rp2_varid) status = nf90_inq_varid(nc%id, nc%j4rp4_varname, nc%j4rp4_varid) @@ -853,11 +869,11 @@ module subroutine swiftest_io_netcdf_open(self, param, readonly) status = nf90_inq_varid(nc%id, nc%ke_spin_varname, nc%KE_spin_varid) status = nf90_inq_varid(nc%id, nc%pe_varname, nc%PE_varid) status = nf90_inq_varid(nc%id, nc%be_varname, nc%BE_varid) - status = nf90_inq_varid(nc%id, nc%L_orb_varname, nc%L_orb_varid) - status = nf90_inq_varid(nc%id, nc%Lspin_varname, nc%Lspin_varid) + status = nf90_inq_varid(nc%id, nc%L_orbit_varname, nc%L_orbit_varid) + status = nf90_inq_varid(nc%id, nc%L_spin_varname, nc%L_spin_varid) status = nf90_inq_varid(nc%id, nc%L_escape_varname, nc%L_escape_varid) - status = nf90_inq_varid(nc%id, nc%Ecollisions_varname, nc%Ecollisions_varid) - status = nf90_inq_varid(nc%id, nc%Euntracked_varname, nc%Euntracked_varid) + status = nf90_inq_varid(nc%id, nc%E_collisions_varname, nc%E_collisions_varid) + status = nf90_inq_varid(nc%id, nc%E_untracked_varname, nc%E_untracked_varid) status = nf90_inq_varid(nc%id, nc%GMescape_varname, nc%GMescape_varid) end if @@ -879,17 +895,16 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier ! Return integer(I4B) :: ierr !! Error code: returns 0 if the read is successful ! Internals - integer(I4B) :: i, tslot, idmax, npl_check, ntp_check, nplm_check, t_max, str_max, status + integer(I4B) :: i, idmax, npl_check, ntp_check, nplm_check, str_max, status real(DP), dimension(:), allocatable :: rtemp real(DP), dimension(:,:), allocatable :: vectemp 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) - associate(cb => self%cb, pl => self%pl, tp => self%tp, npl => self%pl%nbody, ntp => self%tp%nbody) + associate(cb => self%cb, pl => self%pl, tp => self%tp, npl => self%pl%nbody, ntp => self%tp%nbody, tslot => nc%tslot) call pl%setup(npl, param) call tp%setup(ntp, param) @@ -901,7 +916,13 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier allocate(validmask(idmax)) allocate(tpmask(idmax)) allocate(plmask(idmax)) - call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%time_dimid, len=t_max), "netcdf_io_read_frame_system nf90_inquire_dimension time_dimid" ) + call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%time_dimid, len=nc%max_tslot), "netcdf_io_read_frame_system nf90_inquire_dimension time_dimid" ) + if (tslot > nc%max_tslot) then + write(*,*) + write(*,*) "Error in reading frame from NetCDF file. Requested time index value (",tslot,") exceeds maximum time in file (",nc%max_tslot,"). " + call util_exit(FAILURE) + end if + call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%str_dimid, len=str_max), "netcdf_io_read_frame_system nf90_inquire_dimension str_dimid" ) ! First filter out only the id slots that contain valid bodies @@ -1065,14 +1086,14 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier ! end if status = nf90_inq_varid(nc%id, nc%j2rp2_varname, nc%j2rp2_varid) - if (status == nf90_noerr) then + if (status == NF90_NOERR) then call netcdf_io_check( nf90_get_var(nc%id, nc%j2rp2_varid, cb%j2rp2, start=[tslot]), "netcdf_io_read_frame_system nf90_getvar j2rp2_varid" ) else cb%j2rp2 = 0.0_DP end if status = nf90_inq_varid(nc%id, nc%j4rp4_varname, nc%j4rp4_varid) - if (status == nf90_noerr) then + if (status == NF90_NOERR) then call netcdf_io_check( nf90_get_var(nc%id, nc%j4rp4_varid, cb%j4rp4, start=[tslot]), "netcdf_io_read_frame_system nf90_getvar j4rp4_varid" ) else cb%j4rp4 = 0.0_DP @@ -1117,78 +1138,92 @@ module subroutine swiftest_io_netcdf_read_hdr_system(self, nc, param) class(swiftest_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 ! Internals - integer(I4B) :: tslot, status, idmax - real(DP), dimension(:), allocatable :: gmtemp + integer(I4B) :: status, idmax + real(DP), dimension(:), allocatable :: Gmtemp logical, dimension(:), allocatable :: plmask, tpmask, plmmask + integer(I4B), dimension(:), allocatable :: body_status + associate(tslot => nc%tslot) + call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%name_dimid, len=idmax), "netcdf_io_read_hdr_system nf90_inquire_dimension name_dimid" ) + call netcdf_io_check( nf90_get_var(nc%id, nc%time_varid, self%t, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar time_varid" ) - tslot = param%ioutput - call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%name_dimid, len=idmax), "netcdf_io_read_hdr_system nf90_inquire_dimension name_dimid" ) - call netcdf_io_check( nf90_get_var(nc%id, nc%time_varid, self%t, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar time_varid" ) - - allocate(gmtemp(idmax)) - allocate(tpmask(idmax)) - allocate(plmask(idmax)) - allocate(plmmask(idmax)) + allocate(Gmtemp(idmax)) + allocate(tpmask(idmax)) + allocate(plmask(idmax)) + allocate(plmmask(idmax)) + allocate(body_status(idmax)) - call netcdf_io_check( nf90_get_var(nc%id, nc%Gmass_varid, gmtemp, start=[1,1], count=[idmax,1]), "netcdf_io_read_hdr_system nf90_getvar Gmass_varid" ) + call netcdf_io_check( nf90_get_var(nc%id, nc%Gmass_varid, Gmtemp, start=[1,1], count=[idmax,1]), "netcdf_io_read_hdr_system nf90_getvar Gmass_varid" ) + status = nf90_inq_varid(nc%id, nc%status_varname, nc%status_varid) + if (status == NF90_NOERR) then + call netcdf_io_check( nf90_get_var(nc%id, nc%status_varid, body_status, start=[1,tslot], count=[idmax,1]), "netcdf_io_read_hdr_system nf90_getvar status_varid" ) + else + body_status(:) = ACTIVE + end if - plmask(:) = gmtemp(:) == gmtemp(:) - tpmask(:) = .not. plmask(:) - plmask(1) = .false. ! This is the central body - plmmask(:) = plmask(:) + plmask(:) = Gmtemp(:) == Gmtemp(:) + tpmask(:) = .not. plmask(:) + plmask(1) = .false. ! This is the central body - if (param%lmtiny_pl) then - where(plmask(:)) - plmmask(:) = gmtemp(:) > param%GMTINY - endwhere - end if + ! Select only active bodies + plmask(:) = plmask(:) .and. (body_status(:) /= INACTIVE) + tpmask(:) = tpmask(:) .and. (body_status(:) /= INACTIVE) - status = nf90_inq_varid(nc%id, nc%npl_varname, nc%npl_varid) - if (status == nf90_noerr) then - call netcdf_io_check( nf90_get_var(nc%id, nc%npl_varid, self%pl%nbody, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar npl_varid" ) - else - self%pl%nbody = count(plmask(:)) - end if + if (param%lmtiny_pl) then + where(plmask(:)) + plmmask(:) = Gmtemp(:) > param%GMTINY + endwhere + else + plmmask(:) = plmask(:) + end if - status = nf90_inq_varid(nc%id, nc%ntp_varname, nc%ntp_varid) - if (status == nf90_noerr) then - call netcdf_io_check( nf90_get_var(nc%id, nc%ntp_varid, self%tp%nbody, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar ntp_varid" ) - else - self%tp%nbody = count(tpmask(:)) - end if + status = nf90_inq_varid(nc%id, nc%npl_varname, nc%npl_varid) + if (status == NF90_NOERR) then + call netcdf_io_check( nf90_get_var(nc%id, nc%npl_varid, self%pl%nbody, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar npl_varid" ) + else + self%pl%nbody = count(plmask(:)) + end if - if (param%lmtiny_pl) then - status = nf90_inq_varid(nc%id, nc%nplm_varname, nc%nplm_varid) - if (status == nf90_noerr) then - call netcdf_io_check( nf90_get_var(nc%id, nc%nplm_varid, self%pl%nplm, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar nplm_varid" ) + status = nf90_inq_varid(nc%id, nc%ntp_varname, nc%ntp_varid) + if (status == NF90_NOERR) then + call netcdf_io_check( nf90_get_var(nc%id, nc%ntp_varid, self%tp%nbody, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar ntp_varid" ) else - self%pl%nplm = count(plmmask(:)) + self%tp%nbody = count(tpmask(:)) end if - end if - if (param%lenergy) then - status = nf90_inq_varid(nc%id, nc%ke_orb_varname, nc%KE_orb_varid) - if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%KE_orb_varid, self%ke_orbit, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar KE_orb_varid" ) - status = nf90_inq_varid(nc%id, nc%ke_spin_varname, nc%KE_spin_varid) - if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%KE_spin_varid, self%ke_spin, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar KE_spin_varid" ) - status = nf90_inq_varid(nc%id, nc%pe_varname, nc%PE_varid) - if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%PE_varid, self%pe, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar PE_varid" ) - status = nf90_inq_varid(nc%id, nc%be_varname, nc%BE_varid) - if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%BE_varid, self%be, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar BE_varid" ) - status = nf90_inq_varid(nc%id, nc%L_orb_varname, nc%L_orb_varid) - if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%L_orb_varid, self%Lorbit(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_read_hdr_system nf90_getvar L_orb_varid" ) - status = nf90_inq_varid(nc%id, nc%Lspin_varname, nc%Lspin_varid) - if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%Lspin_varid, self%Lspin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_read_hdr_system nf90_getvar Lspin_varid" ) - status = nf90_inq_varid(nc%id, nc%L_escape_varname, nc%L_escape_varid) - if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1, tslot], count=[NDIM,1]), "netcdf_io_read_hdr_system nf90_getvar L_escape_varid" ) - status = nf90_inq_varid(nc%id, nc%Ecollisions_varname, nc%Ecollisions_varid) - if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%Ecollisions_varid, self%Ecollisions, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar Ecollisions_varid" ) - status = nf90_inq_varid(nc%id, nc%Euntracked_varname, nc%Euntracked_varid) - if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%Euntracked_varid, self%Euntracked, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar Euntracked_varid" ) - status = nf90_inq_varid(nc%id, nc%GMescape_varname, nc%GMescape_varid) - if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%GMescape_varid, self%GMescape, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar GMescape_varid" ) - end if + if (param%lmtiny_pl) then + status = nf90_inq_varid(nc%id, nc%nplm_varname, nc%nplm_varid) + if (status == NF90_NOERR) then + call netcdf_io_check( nf90_get_var(nc%id, nc%nplm_varid, self%pl%nplm, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar nplm_varid" ) + else + self%pl%nplm = count(plmmask(:)) + end if + end if + + if (param%lenergy) then + status = nf90_inq_varid(nc%id, nc%ke_orb_varname, nc%KE_orb_varid) + if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%KE_orb_varid, self%ke_orbit, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar KE_orb_varid" ) + status = nf90_inq_varid(nc%id, nc%ke_spin_varname, nc%KE_spin_varid) + if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%KE_spin_varid, self%ke_spin, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar KE_spin_varid" ) + status = nf90_inq_varid(nc%id, nc%pe_varname, nc%PE_varid) + if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%PE_varid, self%pe, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar PE_varid" ) + status = nf90_inq_varid(nc%id, nc%be_varname, nc%BE_varid) + if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%BE_varid, self%be, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar BE_varid" ) + status = nf90_inq_varid(nc%id, nc%L_orbit_varname, nc%L_orbit_varid) + if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%L_orbit_varid, self%L_orbit(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_read_hdr_system nf90_getvar L_orbit_varid" ) + status = nf90_inq_varid(nc%id, nc%L_spin_varname, nc%L_spin_varid) + if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%L_spin_varid, self%L_spin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_read_hdr_system nf90_getvar L_spin_varid" ) + status = nf90_inq_varid(nc%id, nc%L_escape_varname, nc%L_escape_varid) + if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%L_escape_varid, self%L_escape(:), start=[1, tslot], count=[NDIM,1]), "netcdf_io_read_hdr_system nf90_getvar L_escape_varid" ) + status = nf90_inq_varid(nc%id, nc%E_collisions_varname, nc%E_collisions_varid) + if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%E_collisions_varid, self%E_collisions, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar E_collisions_varid" ) + status = nf90_inq_varid(nc%id, nc%E_untracked_varname, nc%E_untracked_varid) + if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%E_untracked_varid, self%E_untracked, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar E_untracked_varid" ) + status = nf90_inq_varid(nc%id, nc%GMescape_varname, nc%GMescape_varid) + if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%GMescape_varid, self%GMescape, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar GMescape_varid" ) + end if + + end associate return end subroutine swiftest_io_netcdf_read_hdr_system @@ -1260,7 +1295,7 @@ module subroutine swiftest_io_netcdf_read_particle_info_system(self, nc, param, end do status = nf90_get_var(nc%id, nc%ptype_varid, ctemp, count=[NAMELEN, idmax]) - if (status /= nf90_noerr) then ! Set default particle types + if (status /= NF90_NOERR) then ! Set default particle types call cb%info%set_value(particle_type=CB_TYPE_NAME) ! Handle semi-interacting bodies in SyMBA @@ -1289,7 +1324,7 @@ module subroutine swiftest_io_netcdf_read_particle_info_system(self, nc, param, if (param%lclose) then status = nf90_inq_varid(nc%id, nc%origin_type_varname, nc%origin_type_varid) - if (status == nf90_noerr) then + if (status == NF90_NOERR) then call netcdf_io_check( nf90_get_var(nc%id, nc%origin_type_varid, ctemp, count=[NAMELEN, idmax]), "netcdf_io_read_particle_info_system nf90_getvar origin_type_varid" ) else ctemp = "Initial Conditions" @@ -1304,7 +1339,7 @@ module subroutine swiftest_io_netcdf_read_particle_info_system(self, nc, param, end do status = nf90_inq_varid(nc%id, nc%origin_time_varname, nc%origin_time_varid) - if (status == nf90_noerr) then + if (status == NF90_NOERR) then call netcdf_io_check( nf90_get_var(nc%id, nc%origin_time_varid, rtemp), "netcdf_io_read_particle_info_system nf90_getvar origin_time_varid" ) else rtemp = param%t0 @@ -1319,7 +1354,7 @@ module subroutine swiftest_io_netcdf_read_particle_info_system(self, nc, param, end do status = nf90_inq_varid(nc%id, nc%origin_rh_varname, nc%origin_rh_varid) - if (status == nf90_noerr) then + if (status == NF90_NOERR) then call netcdf_io_check( nf90_get_var(nc%id, nc%origin_rh_varid, vectemp(:,:)), "netcdf_io_read_particle_info_system nf90_getvar origin_rh_varid" ) else if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then call netcdf_io_check( nf90_get_var(nc%id, nc%rh_varid, vectemp(:,:)), "netcdf_io_read_particle_info_system nf90_getvar rh_varid" ) @@ -1335,7 +1370,7 @@ module subroutine swiftest_io_netcdf_read_particle_info_system(self, nc, param, end do status = nf90_inq_varid(nc%id, nc%origin_vh_varname, nc%origin_vh_varid) - if (status == nf90_noerr) then + if (status == NF90_NOERR) then call netcdf_io_check( nf90_get_var(nc%id, nc%origin_vh_varid, vectemp(:,:)), "netcdf_io_read_particle_info_system nf90_getvar origin_vh_varid" ) else if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then call netcdf_io_check( nf90_get_var(nc%id, nc%vh_varid, vectemp(:,:)), "netcdf_io_read_particle_info_system nf90_getvar vh_varid" ) @@ -1351,7 +1386,7 @@ module subroutine swiftest_io_netcdf_read_particle_info_system(self, nc, param, end do status = nf90_inq_varid(nc%id, nc%collision_id_varname, nc%collision_id_varid) - if (status == nf90_noerr) then + if (status == NF90_NOERR) then call netcdf_io_check( nf90_get_var(nc%id, nc%collision_id_varid, itemp), "netcdf_io_read_particle_info_system nf90_getvar collision_id_varid" ) do i = 1, npl call pl%info(i)%set_value(collision_id=itemp(plind(i))) @@ -1362,7 +1397,7 @@ module subroutine swiftest_io_netcdf_read_particle_info_system(self, nc, param, end if status = nf90_inq_varid(nc%id, nc%discard_time_varname, nc%discard_time_varid) - if (status == nf90_noerr) then + if (status == NF90_NOERR) then call netcdf_io_check( nf90_get_var(nc%id, nc%discard_time_varid, rtemp), "netcdf_io_read_particle_info_system nf90_getvar discard_time_varid" ) call cb%info%set_value(discard_time=rtemp(1)) do i = 1, npl @@ -1374,7 +1409,7 @@ module subroutine swiftest_io_netcdf_read_particle_info_system(self, nc, param, end if status = nf90_inq_varid(nc%id, nc%discard_rh_varname, nc%discard_rh_varid) - if (status == nf90_noerr) then + if (status == NF90_NOERR) then call netcdf_io_check( nf90_get_var(nc%id, nc%discard_rh_varid, vectemp(:,:)), "netcdf_io_read_particle_info_system nf90_getvar discard_rh_varid" ) do i = 1, npl call pl%info(i)%set_value(discard_rh=vectemp(:,plind(i))) @@ -1385,7 +1420,7 @@ module subroutine swiftest_io_netcdf_read_particle_info_system(self, nc, param, end if status = nf90_inq_varid(nc%id, nc%discard_vh_varname, nc%discard_vh_varid) - if (status == nf90_noerr) then + if (status == NF90_NOERR) then call netcdf_io_check( nf90_get_var(nc%id, nc%discard_vh_varid, vectemp(:,:)), "netcdf_io_read_particle_info_system nf90_getvar discard_vh_varid" ) do i = 1, npl call pl%info(i)%set_value(discard_vh=vectemp(:,plind(i))) @@ -1409,25 +1444,24 @@ module subroutine swiftest_io_netcdf_write_frame_body(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_body), intent(in) :: self !! Swiftest base object + class(swiftest_body), intent(in) :: self !! Swiftest base object class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - integer(I4B) :: i, j, tslot, idslot, old_mode + integer(I4B) :: i, j,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 - tslot = param%ioutput call self%write_info(nc, param) - call netcdf_io_check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "netcdf_io_write_frame_body nf90_set_fill" ) + call netcdf_io_check( nf90_set_fill(nc%id, NF90_NOFILL, old_mode), "netcdf_io_write_frame_body nf90_set_fill" ) select type(self) class is (swiftest_body) select type (param) class is (swiftest_parameters) - associate(n => self%nbody) + associate(n => self%nbody, tslot => nc%tslot) if (n == 0) return call swiftest_util_sort(self%id(1:n), ind) @@ -1509,31 +1543,33 @@ module subroutine swiftest_io_netcdf_write_frame_cb(self, nc, param) !! Write a frame of output of the central body implicit none ! Arguments - class(swiftest_cb), intent(in) :: self !! Swiftest base object + class(swiftest_cb), intent(in) :: self !! Swiftest base object class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - integer(I4B) :: i, j, tslot, idslot, old_mode + integer(I4B) :: idslot, old_mode - tslot = param%ioutput - call self%write_info(nc, param) + associate(tslot => nc%tslot) + call self%write_info(nc, param) - call netcdf_io_check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "netcdf_io_write_frame_cb nf90_set_fill" ) + call netcdf_io_check( nf90_set_fill(nc%id, NF90_NOFILL, old_mode), "netcdf_io_write_frame_cb nf90_set_fill" ) - idslot = self%id + 1 - call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, self%id, start=[idslot]), "netcdf_io_write_frame_cb nf90_put_var cb id_varid" ) + idslot = self%id + 1 + call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, self%id, start=[idslot]), "netcdf_io_write_frame_cb nf90_put_var cb id_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%status_varid, ACTIVE, start=[idslot, tslot]), "netcdf_io_write_frame_cb nf90_put_var cb id_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%Gmass_varid, self%Gmass, start=[idslot, tslot]), "netcdf_io_write_frame_cb nf90_put_var cb Gmass_varid" ) - if (param%lclose) call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, self%radius, start=[idslot, tslot]), "netcdf_io_write_frame_cb nf90_put_var cb radius_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%j2rp2_varid, self%j2rp2, start=[tslot]), "netcdf_io_write_frame_cb nf90_put_var cb j2rp2_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%j4rp4_varid, self%j4rp4, start=[tslot]), "netcdf_io_write_frame_cb nf90_put_var cb j4rp4_varid" ) - if (param%lrotation) then - call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, self%Ip(:), start=[1, idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_cb nf90_put_var cb Ip_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, self%rot(:), start=[1, idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_cby nf90_put_var cb rot_varid" ) - end if + call netcdf_io_check( nf90_put_var(nc%id, nc%Gmass_varid, self%Gmass, start=[idslot, tslot]), "netcdf_io_write_frame_cb nf90_put_var cb Gmass_varid" ) + if (param%lclose) call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, self%radius, start=[idslot, tslot]), "netcdf_io_write_frame_cb nf90_put_var cb radius_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%j2rp2_varid, self%j2rp2, start=[tslot]), "netcdf_io_write_frame_cb nf90_put_var cb j2rp2_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%j4rp4_varid, self%j4rp4, start=[tslot]), "netcdf_io_write_frame_cb nf90_put_var cb j4rp4_varid" ) + if (param%lrotation) then + call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, self%Ip(:), start=[1, idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_cb nf90_put_var cb Ip_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, self%rot(:), start=[1, idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_cby nf90_put_var cb rot_varid" ) + end if - call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode), "netcdf_io_write_frame_cb nf90_set_fill old_mode" ) + call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode), "netcdf_io_write_frame_cb nf90_set_fill old_mode" ) + end associate return end subroutine swiftest_io_netcdf_write_frame_cb @@ -1569,28 +1605,28 @@ module subroutine swiftest_io_netcdf_write_hdr_system(self, nc, param) class(swiftest_nbody_system), intent(in) :: self !! Swiftest nbody system object class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: tslot - - tslot = param%ioutput - - call netcdf_io_check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var time_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%npl_varid, self%pl%nbody, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var npl_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%ntp_varid, self%tp%nbody, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var ntp_varid" ) - if (param%lmtiny_pl) call netcdf_io_check( nf90_put_var(nc%id, nc%nplm_varid, self%pl%nplm, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var nplm_varid" ) - - if (param%lenergy) then - call netcdf_io_check( nf90_put_var(nc%id, nc%KE_orb_varid, self%ke_orbit, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var KE_orb_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%KE_spin_varid, self%ke_spin, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var KE_spin_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%PE_varid, self%pe, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var PE_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%BE_varid, self%be, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var BE_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%L_orb_varid, self%Lorbit(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_write_hdr_system nf90_put_var L_orb_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%Lspin_varid, self%Lspin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_write_hdr_system nf90_put_var Lspin_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_write_hdr_system nf90_put_var L_escape_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%Ecollisions_varid, self%Ecollisions, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var Ecollisions_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%Euntracked_varid, self%Euntracked, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var Euntracked_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%GMescape_varid, self%GMescape, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var GMescape_varid" ) - end if + + associate(tslot => nc%tslot) + + call netcdf_io_check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var time_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%npl_varid, self%pl%nbody, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var npl_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%ntp_varid, self%tp%nbody, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var ntp_varid" ) + if (param%lmtiny_pl) call netcdf_io_check( nf90_put_var(nc%id, nc%nplm_varid, self%pl%nplm, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var nplm_varid" ) + + if (param%lenergy) then + call netcdf_io_check( nf90_put_var(nc%id, nc%KE_orb_varid, self%ke_orbit, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var KE_orb_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%KE_spin_varid, self%ke_spin, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var KE_spin_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%PE_varid, self%pe, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var PE_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%BE_varid, self%be, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var BE_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%L_orbit_varid, self%L_orbit(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_write_hdr_system nf90_put_var L_orbit_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%L_spin_varid, self%L_spin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_write_hdr_system nf90_put_var L_spin_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%L_escape_varid, self%L_escape(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_write_hdr_system nf90_put_var L_escape_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%E_collisions_varid, self%E_collisions, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var E_collisions_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%E_untracked_varid, self%E_untracked, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var E_untracked_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%GMescape_varid, self%GMescape, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var GMescape_varid" ) + end if + + end associate return end subroutine swiftest_io_netcdf_write_hdr_system @@ -1599,23 +1635,23 @@ end subroutine swiftest_io_netcdf_write_hdr_system module subroutine swiftest_io_netcdf_write_info_body(self, nc, param) !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton !! - !! Write all current particle to file + !! Write all current particle information metadata to file implicit none ! Arguments - class(swiftest_body), intent(in) :: self !! Swiftest particle object + class(swiftest_body), intent(in) :: self !! Swiftest particle object class(swiftest_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_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i, j, idslot, old_mode integer(I4B), dimension(:), allocatable :: ind character(len=:), allocatable :: charstring ! This string of spaces of length NAMELEN is used to clear out any old data left behind inside the string variables - call netcdf_io_check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "netcdf_io_write_info_body nf90_set_fill nf90_nofill" ) + call netcdf_io_check( nf90_set_fill(nc%id, NF90_NOFILL, old_mode), "netcdf_io_write_info_body nf90_set_fill NF90_NOFILL" ) select type(self) class is (swiftest_body) - associate(n => self%nbody) + associate(n => self%nbody, tslot => nc%tslot) if (n == 0) return call swiftest_util_sort(self%id(1:n), ind) @@ -1623,6 +1659,7 @@ module subroutine swiftest_io_netcdf_write_info_body(self, nc, param) j = ind(i) idslot = self%id(j) + 1 call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, self%id(j), start=[idslot]), "netcdf_io_write_info_body nf90_put_var id_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%status_varid, self%status(j), start=[idslot,tslot]), "netcdf_io_write_info_body nf90_put_var status_varid" ) charstring = trim(adjustl(self%info(j)%name)) call netcdf_io_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "netcdf_io_write_info_body nf90_put_var name_varid" ) @@ -1665,7 +1702,7 @@ module subroutine swiftest_io_netcdf_write_info_cb(self, nc, param) character(len=:), allocatable :: charstring ! This string of spaces of length NAMELEN is used to clear out any old data left behind inside the string variables - call netcdf_io_check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "netcdf_io_write_info_cb nf90_set_fill nf90_nofill" ) + call netcdf_io_check( nf90_set_fill(nc%id, NF90_NOFILL, old_mode), "netcdf_io_write_info_cb nf90_set_fill NF90_NOFILL" ) idslot = self%id + 1 call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, self%id, start=[idslot]), "netcdf_io_write_info_cb nf90_put_var id_varid" ) @@ -1852,43 +1889,43 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i call swiftest_io_toupper(param_value) if (param_value == "NO" .or. param_value == 'F') param%lfirstenergy = .false. case("EORBIT_ORIG") - read(param_value, *, err = 667, iomsg = iomsg) param%Eorbit_orig + read(param_value, *, err = 667, iomsg = iomsg) param%E_orbit_orig case("GMTOT_ORIG") read(param_value, *, err = 667, iomsg = iomsg) param%GMtot_orig case("LTOT_ORIG") - read(param_value, *, err = 667, iomsg = iomsg) param%Ltot_orig(1) + read(param_value, *, err = 667, iomsg = iomsg) param%L_total_orig(1) do i = 2, NDIM ifirst = ilast + 2 param_value = swiftest_io_get_token(line, ifirst, ilast, iostat) - read(param_value, *, err = 667, iomsg = iomsg) param%Ltot_orig(i) + read(param_value, *, err = 667, iomsg = iomsg) param%L_total_orig(i) end do case("LORBIT_ORIG") - read(param_value, *, err = 667, iomsg = iomsg) param%Lorbit_orig(1) + read(param_value, *, err = 667, iomsg = iomsg) param%L_orbit_orig(1) do i = 2, NDIM ifirst = ilast + 2 param_value = swiftest_io_get_token(line, ifirst, ilast, iostat) - read(param_value, *, err = 667, iomsg = iomsg) param%Lorbit_orig(i) + read(param_value, *, err = 667, iomsg = iomsg) param%L_orbit_orig(i) end do case("LSPIN_ORIG") - read(param_value, *, err = 667, iomsg = iomsg) param%Lspin_orig(1) + read(param_value, *, err = 667, iomsg = iomsg) param%L_spin_orig(1) do i = 2, NDIM ifirst = ilast + 2 param_value = swiftest_io_get_token(line, ifirst, ilast, iostat) - read(param_value, *, err = 667, iomsg = iomsg) param%Lspin_orig(i) + read(param_value, *, err = 667, iomsg = iomsg) param%L_spin_orig(i) end do case("LESCAPE") - read(param_value, *, err = 667, iomsg = iomsg) param%Lescape(1) + read(param_value, *, err = 667, iomsg = iomsg) param%L_escape(1) do i = 2, NDIM ifirst = ilast + 2 param_value = swiftest_io_get_token(line, ifirst, ilast, iostat) - read(param_value, *, err = 667, iomsg = iomsg) param%Lescape(i) + read(param_value, *, err = 667, iomsg = iomsg) param%L_escape(i) end do case("GMESCAPE") read(param_value, *, err = 667, iomsg = iomsg) param%GMescape case("ECOLLISIONS") - read(param_value, *, err = 667, iomsg = iomsg) param%Ecollisions + read(param_value, *, err = 667, iomsg = iomsg) param%E_collisions case("EUNTRACKED") - read(param_value, *, err = 667, iomsg = iomsg) param%Euntracked + read(param_value, *, err = 667, iomsg = iomsg) param%E_untracked case ("MAXID") read(param_value, *, err = 667, iomsg = iomsg) param%maxid case ("MAXID_COLLISION") @@ -2026,7 +2063,6 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i ! Calculate the G for the nbody_system units param%GU = GC / (param%DU2M**3 / (param%MU2KG * param%TU2S**2)) - ! A minimal log of collision outcomes is stored in the following log file ! More complete data on collisions is stored in the NetCDF output files call swiftest_io_log_start(param, COLLISION_LOG_OUT, "Collision logfile") @@ -2163,6 +2199,8 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i end select iostat = 0 + + call param%set_display(param%display_style) ! Print the contents of the parameter file to standard output call param%writer(unit = param%display_unit, iotype = "none", v_list = [0], iostat = iostat, iomsg = iomsg) @@ -2584,14 +2622,14 @@ module subroutine swiftest_io_read_in_system(self, param) call self%pl%read_in(param) call self%tp%read_in(param) ! Copy over param file variable inputs - self%Eorbit_orig = param%Eorbit_orig + self%E_orbit_orig = param%E_orbit_orig self%GMtot_orig = param%GMtot_orig - self%Ltot_orig(:) = param%Ltot_orig(:) - self%Lorbit_orig(:) = param%Lorbit_orig(:) - self%Lspin_orig(:) = param%Lspin_orig(:) - self%Lescape(:) = param%Lescape(:) - self%Ecollisions = param%Ecollisions - self%Euntracked = param%Euntracked + self%L_total_orig(:) = param%L_total_orig(:) + self%L_orbit_orig(:) = param%L_orbit_orig(:) + self%L_spin_orig(:) = param%L_spin_orig(:) + self%L_escape(:) = param%L_escape(:) + self%E_collisions = param%E_collisions + self%E_untracked = param%E_untracked else allocate(tmp_param, source=param) tmp_param%system_history%nc%file_name = param%nc_in @@ -2737,7 +2775,6 @@ module subroutine swiftest_io_read_in_param(self, param_file_name) ! Read in name of parameter file self%param_file_name = trim(adjustl(param_file_name)) - write(self%display_unit, *) 'Parameter input file is ' // self%param_file_name !! todo: Currently this procedure does not work in user-defined derived-type input mode !! as the newline characters are ignored in the input file when compiled in ifort. @@ -2762,14 +2799,16 @@ module subroutine swiftest_io_set_display_param(self, display_style) class(swiftest_parameters), intent(inout) :: self !! Current run configuration parameters character(*), intent(in) :: display_style !! Style of the output display ! Internals - character(STRMAX) :: errmsg + character(STRMAX) :: errmsg + logical :: fileExists select case(display_style) case ('STANDARD') self%display_unit = OUTPUT_UNIT !! stdout from iso_fortran_env self%log_output = .false. case ('COMPACT', 'PROGRESS') - if (self%lrestart) then + inquire(SWIFTEST_LOG_FILE, exist=fileExists) + if (self%lrestart.and.fileExists) then open(unit=SWIFTEST_LOG_OUT, file=SWIFTEST_LOG_FILE, status="OLD", position="APPEND", err = 667, iomsg = errmsg) else open(unit=SWIFTEST_LOG_OUT, file=SWIFTEST_LOG_FILE, status="REPLACE", err = 667, iomsg = errmsg) @@ -2858,8 +2897,6 @@ module subroutine swiftest_io_write_frame_system(self, param) logical :: fileExists associate (nc => param%system_history%nc, pl => self%pl, tp => self%tp, npl => self%pl%nbody, ntp => self%tp%nbody) - nc%name_chunk = npl + ntp - nc%time_chunk = max(param%dump_cadence / param%istep_out, 1) nc%file_name = param%outfile if (lfirst) then inquire(file=param%outfile, exist=fileExists) diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index 019768166..f41dcbef7 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -334,36 +334,36 @@ module swiftest real(DP) :: be = 0.0_DP !! nbody_system binding energy of all bodies real(DP) :: te = 0.0_DP !! nbody_system total energy real(DP) :: oblpot = 0.0_DP !! nbody_system potential energy due to oblateness of the central body - real(DP), dimension(NDIM) :: Lorbit = 0.0_DP !! nbody_system orbital angular momentum vector - real(DP), dimension(NDIM) :: Lspin = 0.0_DP !! nbody_system spin angular momentum vector - real(DP), dimension(NDIM) :: Ltot = 0.0_DP !! nbody_system angular momentum vector + real(DP), dimension(NDIM) :: L_orbit = 0.0_DP !! nbody_system orbital angular momentum vector + real(DP), dimension(NDIM) :: L_spin = 0.0_DP !! nbody_system spin angular momentum vector + real(DP), dimension(NDIM) :: L_total = 0.0_DP !! nbody_system angular momentum vector real(DP) :: ke_orbit_orig = 0.0_DP !! Initial orbital kinetic energy real(DP) :: ke_spin_orig = 0.0_DP !! Initial spin kinetic energy real(DP) :: pe_orig = 0.0_DP !! Initial potential energy real(DP) :: be_orig = 0.0_DP !! Initial binding energy - real(DP) :: Eorbit_orig = 0.0_DP !! Initial orbital energy + real(DP) :: E_orbit_orig = 0.0_DP !! Initial orbital energy real(DP) :: GMtot_orig = 0.0_DP !! Initial nbody_system mass - real(DP), dimension(NDIM) :: Ltot_orig = 0.0_DP !! Initial total angular momentum vector - real(DP), dimension(NDIM) :: Lorbit_orig = 0.0_DP !! Initial orbital angular momentum - real(DP), dimension(NDIM) :: Lspin_orig = 0.0_DP !! Initial spin angular momentum vector - real(DP), dimension(NDIM) :: Lescape = 0.0_DP !! Angular momentum of bodies that escaped the nbody_system (used for bookeeping) + real(DP), dimension(NDIM) :: L_total_orig = 0.0_DP !! Initial total angular momentum vector + real(DP), dimension(NDIM) :: L_orbit_orig = 0.0_DP !! Initial orbital angular momentum + real(DP), dimension(NDIM) :: L_spin_orig = 0.0_DP !! Initial spin angular momentum vector + real(DP), dimension(NDIM) :: L_escape = 0.0_DP !! Angular momentum of bodies that escaped the nbody_system (used for bookeeping) real(DP) :: GMescape = 0.0_DP !! Mass of bodies that escaped the nbody_system (used for bookeeping) - real(DP) :: Ecollisions = 0.0_DP !! Energy lost from nbody_system due to collisions - real(DP) :: Euntracked = 0.0_DP !! Energy gained from nbody_system due to escaped bodies + real(DP) :: E_collisions = 0.0_DP !! Energy lost from nbody_system due to collisions + real(DP) :: E_untracked = 0.0_DP !! Energy gained from nbody_system due to escaped bodies ! Energy, momentum, and mass errors (used in error reporting) real(DP) :: ke_orbit_error = 0.0_DP real(DP) :: ke_spin_error = 0.0_DP real(DP) :: pe_error = 0.0_DP real(DP) :: be_error = 0.0_DP - real(DP) :: Eorbit_error = 0.0_DP + real(DP) :: E_orbit_error = 0.0_DP real(DP) :: Ecoll_error = 0.0_DP - real(DP) :: Euntracked_error = 0.0_DP - real(DP) :: Etot_error = 0.0_DP - real(DP) :: Lorbit_error = 0.0_DP - real(DP) :: Lspin_error = 0.0_DP - real(DP) :: Lescape_error = 0.0_DP - real(DP) :: Ltot_error = 0.0_DP + real(DP) :: E_untracked_error = 0.0_DP + real(DP) :: te_error = 0.0_DP + real(DP) :: L_orbit_error = 0.0_DP + real(DP) :: L_spin_error = 0.0_DP + real(DP) :: L_escape_error = 0.0_DP + real(DP) :: L_total_error = 0.0_DP real(DP) :: Mtot_error = 0.0_DP real(DP) :: Mescape_error = 0.0_DP diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 9d34b1af5..da750e956 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -1242,9 +1242,9 @@ module subroutine swiftest_util_get_energy_momentum_system(self, param) real(DP) :: hx, hy, hz associate(nbody_system => self, pl => self%pl, npl => self%pl%nbody, cb => self%cb) - nbody_system%Lorbit(:) = 0.0_DP - nbody_system%Lspin(:) = 0.0_DP - nbody_system%Ltot(:) = 0.0_DP + nbody_system%L_orbit(:) = 0.0_DP + nbody_system%L_spin(:) = 0.0_DP + nbody_system%L_total(:) = 0.0_DP nbody_system%ke_orbit = 0.0_DP nbody_system%ke_spin = 0.0_DP @@ -1312,20 +1312,20 @@ module subroutine swiftest_util_get_energy_momentum_system(self, param) nbody_system%ke_orbit = 0.5_DP * (kecb + sum(kepl(1:npl), pl%lmask(1:npl))) if (param%lrotation) nbody_system%ke_spin = 0.5_DP * (kespincb + sum(kespinpl(1:npl), pl%lmask(1:npl))) - nbody_system%Lorbit(1) = Lcborbit(1) + sum(Lplorbitx(1:npl), pl%lmask(1:npl)) - nbody_system%Lorbit(2) = Lcborbit(2) + sum(Lplorbity(1:npl), pl%lmask(1:npl)) - nbody_system%Lorbit(3) = Lcborbit(3) + sum(Lplorbitz(1:npl), pl%lmask(1:npl)) + nbody_system%L_orbit(1) = Lcborbit(1) + sum(Lplorbitx(1:npl), pl%lmask(1:npl)) + nbody_system%L_orbit(2) = Lcborbit(2) + sum(Lplorbity(1:npl), pl%lmask(1:npl)) + nbody_system%L_orbit(3) = Lcborbit(3) + sum(Lplorbitz(1:npl), pl%lmask(1:npl)) if (param%lrotation) then - nbody_system%Lspin(1) = Lcbspin(1) + sum(Lplspinx(1:npl), pl%lmask(1:npl)) - nbody_system%Lspin(2) = Lcbspin(2) + sum(Lplspiny(1:npl), pl%lmask(1:npl)) - nbody_system%Lspin(3) = Lcbspin(3) + sum(Lplspinz(1:npl), pl%lmask(1:npl)) + nbody_system%L_spin(1) = Lcbspin(1) + sum(Lplspinx(1:npl), pl%lmask(1:npl)) + nbody_system%L_spin(2) = Lcbspin(2) + sum(Lplspiny(1:npl), pl%lmask(1:npl)) + nbody_system%L_spin(3) = Lcbspin(3) + sum(Lplspinz(1:npl), pl%lmask(1:npl)) end if nbody_system%be = sum(-3*pl%Gmass(1:npl)*pl%mass(1:npl)/(5*pl%radius(1:npl)), pl%lmask(1:npl)) nbody_system%te = nbody_system%ke_orbit + nbody_system%ke_spin + nbody_system%pe + nbody_system%be - nbody_system%Ltot(:) = nbody_system%Lorbit(:) + nbody_system%Lspin(:) + nbody_system%L_total(:) = nbody_system%L_orbit(:) + nbody_system%L_spin(:) end associate return @@ -1717,10 +1717,6 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param) end where end if - call nc%open(param) - call pl%write_info(nc, param) - call nc%close() - ! Reindex the new list of bodies call pl%sort("mass", ascending=.false.) call pl%flatten(param) @@ -2608,7 +2604,11 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param) call encounter_history%reset() select type(nc => encounter_history%nc) class is (encounter_netcdf_parameters) - nc%file_number = param%iloop / param%dump_cadence + nc%file_name = ENCOUNTER_OUTFILE + if (.not.param%lrestart) then + call nc%initialize(param) + call nc%close() + end if end select allocate(nbody_system%encounter_history, source=encounter_history) end if @@ -2617,7 +2617,11 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param) call collision_history%reset() select type(nc => collision_history%nc) class is (collision_netcdf_parameters) - nc%file_number = param%iloop / param%dump_cadence + nc%file_name = COLLISION_OUTFILE + if (.not.param%lrestart) then + call nc%initialize(param) + call nc%close() + end if end select allocate(nbody_system%collision_history, source=collision_history) diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 216ab5f28..88d674fb5 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -116,7 +116,7 @@ subroutine symba_discard_conserve_mtm(pl, nbody_system, param, ipl, lescape_body integer(I4B), intent(in) :: ipl logical, intent(in) :: lescape_body ! Internals - real(DP), dimension(NDIM) :: Lpl, Ltot, Lcb, xcom, vcom + real(DP), dimension(NDIM) :: Lpl, L_total, Lcb, xcom, vcom real(DP) :: pe, be, ke_orbit, ke_spin integer(I4B) :: i, oldstat @@ -142,12 +142,12 @@ subroutine symba_discard_conserve_mtm(pl, nbody_system, param, ipl, lescape_body pe = pe - pl%Gmass(i) * pl%mass(ipl) / norm2(pl%rb(:, ipl) - pl%rb(:, i)) end do - Ltot(:) = 0.0_DP + L_total(:) = 0.0_DP do i = 1, pl%nbody Lpl(:) = pL%mass(i) * (pl%rb(:,i) .cross. pl%vb(:, i)) - Ltot(:) = Ltot(:) + Lpl(:) + L_total(:) = L_total(:) + Lpl(:) end do - Ltot(:) = Ltot(:) + cb%mass * (cb%rb(:) .cross. cb%vb(:)) + L_total(:) = L_total(:) + cb%mass * (cb%rb(:) .cross. cb%vb(:)) call pl%b2h(cb) oldstat = pl%status(ipl) pl%status(ipl) = INACTIVE @@ -156,11 +156,11 @@ subroutine symba_discard_conserve_mtm(pl, nbody_system, param, ipl, lescape_body do i = 1, pl%nbody if (i == ipl) cycle Lpl(:) = pl%mass(i) * (pl%rb(:,i) .cross. pl%vb(:, i)) - Ltot(:) = Ltot(:) - Lpl(:) + L_total(:) = L_total(:) - Lpl(:) end do - Ltot(:) = Ltot(:) - cb%mass * (cb%rb(:) .cross. cb%vb(:)) - nbody_system%Lescape(:) = nbody_system%Lescape(:) + Ltot(:) - if (param%lrotation) nbody_system%Lescape(:) = nbody_system%Lescape + pl%mass(ipl) * pl%radius(ipl)**2 & + L_total(:) = L_total(:) - cb%mass * (cb%rb(:) .cross. cb%vb(:)) + nbody_system%L_escape(:) = nbody_system%L_escape(:) + L_total(:) + if (param%lrotation) nbody_system%L_escape(:) = nbody_system%L_escape + pl%mass(ipl) * pl%radius(ipl)**2 & * pl%Ip(3, ipl) * pl%rot(:, ipl) else @@ -195,8 +195,8 @@ subroutine symba_discard_conserve_mtm(pl, nbody_system, param, ipl, lescape_body ! We must do this for proper book-keeping, since we can no longer track this body's contribution to energy directly if (lescape_body) then - nbody_system%Ecollisions = nbody_system%Ecollisions + ke_orbit + ke_spin + pe + be - nbody_system%Euntracked = nbody_system%Euntracked - (ke_orbit + ke_spin + pe + be) + nbody_system%E_collisions = nbody_system%E_collisions + ke_orbit + ke_spin + pe + be + nbody_system%E_untracked = nbody_system%E_untracked - (ke_orbit + ke_spin + pe + be) end if end select @@ -301,7 +301,7 @@ subroutine symba_discard_peri_pl(pl, nbody_system, param) logical, save :: lfirst = .true. logical :: lfirst_orig integer(I4B) :: i - character(len=STRMAX) :: timestr, idstr + character(len=STRMAX) :: timestr, idstr, message lfirst_orig = pl%lfirst @@ -320,8 +320,9 @@ subroutine symba_discard_peri_pl(pl, nbody_system, param) pl%status(i) = DISCARDED_PERI write(timestr, *) nbody_system%t write(idstr, *) pl%id(i) - write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // & + write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // & ") perihelion distance too small at t = " // trim(adjustl(timestr)) + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) call pl%info(i)%set_value(status="DISCARDED_PERI", discard_time=nbody_system%t, & discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), discard_body_id=nbody_system%cb%id) end if @@ -345,7 +346,7 @@ module subroutine symba_discard_pl(self, nbody_system, param) class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - real(DP) :: Eorbit_before, Eorbit_after + real(DP) :: E_orbit_before, E_orbit_after select type(nbody_system) class is (symba_nbody_system) @@ -362,7 +363,7 @@ module subroutine symba_discard_pl(self, nbody_system, param) if (param%lenergy) then call nbody_system%get_energy_and_momentum(param) - Eorbit_before = nbody_system%te + E_orbit_before = nbody_system%te end if call symba_discard_nonplpl_conservation(self, nbody_system, param) @@ -374,8 +375,8 @@ module subroutine symba_discard_pl(self, nbody_system, param) if (param%lenergy) then call nbody_system%get_energy_and_momentum(param) - Eorbit_after = nbody_system%te - nbody_system%Ecollisions = nbody_system%Ecollisions + (Eorbit_after - Eorbit_before) + E_orbit_after = nbody_system%te + nbody_system%E_collisions = nbody_system%E_collisions + (E_orbit_after - E_orbit_before) end if end associate diff --git a/src/whm/whm_drift.f90 b/src/whm/whm_drift.f90 index 21a092793..5d507c42b 100644 --- a/src/whm/whm_drift.f90 +++ b/src/whm/whm_drift.f90 @@ -27,6 +27,7 @@ module subroutine whm_drift_pl(self, nbody_system, param, dt) ! Internals integer(I4B) :: i integer(I4B), dimension(:), allocatable :: iflag + character(len=STRMAX) :: message if (self%nbody == 0) return @@ -41,11 +42,12 @@ module subroutine whm_drift_pl(self, nbody_system, param, dt) end where do i = 1, npl if (iflag(i) /= 0) then - write(*, *) " Planet ", pl%id(i), " is lost!!!!!!!!!!!!" - WRITE(*, *) pl%muj(i), dt - WRITE(*, *) pl%xj(:,i) - WRITE(*, *) pl%vj(:,i) - WRITE(*, *) " STOPPING " + write(message, *) " Planet ", pl%id(i), " is lost!!!!!!!!!!!!",new_line('A'), & + pl%muj(i), dt, new_line('A'), & + pl%xj(:,i),new_line('A'), & + pl%vj(:,i),new_line('A'), & + " STOPPING " + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) end if end do call util_exit(FAILURE)