diff --git a/.gitignore b/.gitignore index da1f14c2a..2d50953b7 100644 --- a/.gitignore +++ b/.gitignore @@ -2,7 +2,6 @@ * # Now, whitelist anything that's a directory !*/ -!*/python # whitelist only the files that ever need to be tracked !*.f90 !*.sh @@ -13,46 +12,13 @@ !paper/paper.md !paper/paper.bib !README.swifter -!*.in dump* -!example/cleanup -!example/swifter_symba_omp -!Makefile -src/*/Makefile -!Makefile.Defines -src/*/Makefile.Defines -!.gitignore +!**/.gitignore !*.py !*.ipynb +!examples/** *ipynb_checkpoints - -#fxdr library -!fxdr*/*.c -!fxdr*/*.F -!fxdr*/*.3f -!fxdr*/*.h -!fxdr*/*.inc -!fxdr*/cxdrreal64 -!fxdr*/test.orig.xdr -!fxdr*/test_read_only.xdr -!fxdr*/configure/ -!fxdr*/README* -!fxdr*/Makefile.bak -!fxdr*/Makefile.fxdr -!fxdr*/Makefile.tmpl -!fxdr*/Defines.* - -#collresolve -!collresolve/*.c -!collresolve/*.h -!collresolve/*.py -!collresolve/Makefile.am -!collresolve/configure.ac -!collresolve/cambioni2019/*.h -!collresolve/cambioni2019/*.c -collresolve/config.h - - +**/.DS_Store #Documentation !docs/* @@ -60,25 +26,5 @@ collresolve/config.h !docs/*/*/* !README_figs/* -python/swiftest/tests/convert_code_type/swift2swifter/pl.swift2swifter.in - -python/swiftest/tests/convert_code_type/swift2swifter/tp.swift2swifter.in - -python/swiftest/tests/convert_code_type/swift2swiftest/cb.swift2swiftest.in - -python/swiftest/tests/convert_code_type/swift2swiftest/pl.swift2swiftest.in - -python/swiftest/tests/convert_code_type/swift2swiftest/tp.swift2swiftest.in - -python/swiftest/tests/convert_code_type/swifter2swiftest/cb.swifter2swiftest.in - -python/swiftest/tests/convert_code_type/cb.swifter2swiftest.in - -python/swiftest/tests/convert_code_type/swifter2swiftest/pl.swifter2swiftest.in - -python/swiftest/tests/convert_code_type/swifter2swiftest/tp.swifter2swiftest.in - -!python/swiftest/requirements.txt - bin/ build/* diff --git a/examples/.gitignore b/examples/.gitignore index 4aceee79f..7c5c72692 100644 --- a/examples/.gitignore +++ b/examples/.gitignore @@ -1,4 +1,6 @@ -*/param.* -*/simdata/* -*/*/simdata/* - +* +!.gitignore +!Basic_Simulation/* +!Fragmentation/* +!helio_gr_test/* +!whm_gr_test/* \ No newline at end of file diff --git a/examples/Basic_Simulation/.gitignore b/examples/Basic_Simulation/.gitignore new file mode 100644 index 000000000..0a4af3872 --- /dev/null +++ b/examples/Basic_Simulation/.gitignore @@ -0,0 +1,6 @@ +* +!.gitignore +!initial_conditions.py +!output_reader.py +!run_from_file.py +!read_old_run.py diff --git a/examples/Fragmentation/.gitignore b/examples/Fragmentation/.gitignore new file mode 100644 index 000000000..ff09bd225 --- /dev/null +++ b/examples/Fragmentation/.gitignore @@ -0,0 +1,3 @@ +* +!.gitignore +!Fragmentation_Movie.py diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 5d660bcec..80f75a86a 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -40,8 +40,8 @@ movie_titles = dict(zip(available_movie_styles, movie_title_list)) # These initial conditions were generated by trial and error -pos_vectors = {"disruption_headon" : [np.array([1.0, -1.807993e-05, 0.0]), - np.array([1.0, 1.807993e-05 ,0.0])], +pos_vectors = {"disruption_headon" : [np.array([1.0, -2.807993e-05, 0.0]), + np.array([1.0, 2.807993e-05 ,0.0])], "supercatastrophic_off_axis": [np.array([1.0, -4.2e-05, 0.0]), np.array([1.0, 4.2e-05, 0.0])], "hitandrun" : [np.array([1.0, -2.0e-05, 0.0]), @@ -84,7 +84,9 @@ class AnimatedScatter(object): """An animated scatter plot using matplotlib.animations.FuncAnimation.""" def __init__(self, sim, animfile, title, nskip=1): - nframes = int(sim.data['time'].size) + tgood = sim.enc.where(sim.enc['Gmass'] > 9e-8).time + self.ds = sim.enc.sel(time=tgood) + nframes = int(self.ds['time'].size) self.sim = sim self.title = title self.body_color_list = {'Initial conditions': 'xkcd:windows blue', @@ -106,13 +108,14 @@ def setup_plot(self): fig = plt.figure(figsize=figsize, dpi=300) plt.tight_layout(pad=0) - # Calculate the distance along the y-axis between the colliding bodies at the start of the simulation. # This will be used to scale the axis limits on the movie. - scale_frame = abs(sim.data['xhy'].isel(time=0, name=1).values) \ - + abs( sim.data['xhy'].isel(time=0, name=2).values) + rhy1 = self.ds['rh'].sel(name="Body1",space='y').isel(time=0).values[()] + rhy2 = self.ds['rh'].sel(name="Body2",space='y').isel(time=0).values[()] + + scale_frame = abs(rhy1) + abs(rhy2) ax = plt.Axes(fig, [0.1, 0.1, 0.8, 0.8]) - self.ax_pt_size = figsize[0] * 0.8 * 72 / scale_frame + self.ax_pt_size = figsize[0] * 0.8 * 72 / (2 * scale_frame) ax.set_xlim(-scale_frame, scale_frame) ax.set_ylim(-scale_frame, scale_frame) ax.set_xticks([]) @@ -131,19 +134,20 @@ def center(Gmass, x, y): x = x[~np.isnan(x)] y = y[~np.isnan(y)] Gmass = Gmass[~np.isnan(Gmass)] + x = x[Gmass] x_com = np.sum(Gmass * x) / np.sum(Gmass) - y_com = np.sum(Gmass * y) / np.sum(Gmass) + y_com = #np.sum(Gmass * y) / np.sum(Gmass) return x_com, y_com Gmass, rh, point_rad = next(self.data_stream(frame)) - x_com, y_com = center(Gmass, rh[:,1], rh[:,2]) - self.scatter_artist.set_offsets(np.c_[x - x_com, y - y_com]) - self.scatter_artist.set_sizes(point_rad) + #x_com, y_com = center(Gmass, rh[:,0], rh[:,1]) + self.scatter_artist.set_offsets(np.c_[rh[:,0] - x_com, rh[:,1] - y_com]) + self.scatter_artist.set_sizes(point_rad**2) return self.scatter_artist, def data_stream(self, frame=0): while True: - ds = self.sim.data.isel(time=frame) + ds = self.ds.isel(time=frame) ds = ds.where(ds['name'] != "Sun", drop=True) radius = ds['radius'].values Gmass = ds['Gmass'].values @@ -153,47 +157,30 @@ def data_stream(self, frame=0): if __name__ == "__main__": - print("Select a fragmentation movie to generate.") - print("1. Head-on disruption") - print("2. Off-axis supercatastrophic") - print("3. Hit and run") - print("4. All of the above") - user_selection = int(input("? ")) - - if user_selection > 0 and user_selection < 4: - movie_styles = [available_movie_styles[user_selection-1]] - else: - print("Generating all movie styles") - movie_styles = available_movie_styles.copy() - - for style in movie_styles: - bin_file = Path(style) / "bin.nc" - if bin_file.exists(): - user_selection = input(f"An older simulation of {movie_titles[style]} exists. Do you want to re-run it? [y/N] ") - if user_selection == "": - run_new = False - else: - try: - run_new = swiftest.io.str2bool(user_selection) - except: - run_new = False - else: - run_new = True - - movie_filename = f"{style}.mp4" - - # Pull in the Swiftest output data from the parameter file and store it as a Xarray dataset. - if run_new: - sim = swiftest.Simulation(simdir=style, rotation=True, init_cond_format = "XV", compute_conservation_values=True) - sim.add_solar_system_body("Sun") - sim.add_body(Gmass=body_Gmass[style], radius=body_radius[style], rh=pos_vectors[style], vh=vel_vectors[style]) #, rot=rot_vectors[style]) - - # Set fragmentation parameters - minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body - gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades - sim.set_parameter(fragmentation = True, gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) - sim.run(dt=2e-5, tstop=2.e-5) - else: - sim = swiftest.Simulation(param_file=param_file, read_old_output_file=True) - - anim = AnimatedScatter(sim,movie_filename,movie_titles[style],nskip=10) + print("Select a fragmentation movie to generate.") + print("1. Head-on disruption") + print("2. Off-axis supercatastrophic") + print("3. Hit and run") + print("4. All of the above") + user_selection = int(input("? ")) + + if user_selection > 0 and user_selection < 4: + movie_styles = [available_movie_styles[user_selection-1]] + else: + print("Generating all movie styles") + movie_styles = available_movie_styles.copy() + + for style in movie_styles: + movie_filename = f"{style}.mp4" + # Pull in the Swiftest output data from the parameter file and store it as a Xarray dataset. + sim = swiftest.Simulation(simdir=style, rotation=True, init_cond_format = "XV", compute_conservation_values=True) + sim.add_solar_system_body("Sun") + sim.add_body(Gmass=body_Gmass[style], radius=body_radius[style], rh=pos_vectors[style], vh=vel_vectors[style]) #, rot=rot_vectors[style]) + # + # Set fragmentation parameters + minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body + gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades + sim.set_parameter(fragmentation=True, fragmentation_save="TRAJECTORY", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) + sim.run(dt=1e-4, tstop=2.0e-3, istep_out=1, dump_cadence=0) + + anim = AnimatedScatter(sim,movie_filename,movie_titles[style],nskip=1) \ No newline at end of file diff --git a/examples/helio_gr_test/.gitignore b/examples/helio_gr_test/.gitignore new file mode 100644 index 000000000..8968b5dd7 --- /dev/null +++ b/examples/helio_gr_test/.gitignore @@ -0,0 +1,3 @@ +* +!.gitignore +!helio_gr_test.py \ No newline at end of file diff --git a/examples/whm_gr_test/.gitignore b/examples/whm_gr_test/.gitignore new file mode 100644 index 000000000..1463c046c --- /dev/null +++ b/examples/whm_gr_test/.gitignore @@ -0,0 +1,3 @@ +* +!.gitignore +!whm_gr_test.py \ No newline at end of file diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 40ce7ae33..f03fb622f 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -31,7 +31,11 @@ "INTERACTION_LOOPS", "ENCOUNTER_CHECK", "TSTART", - "DUMP_CADENCE") + "DUMP_CADENCE", + "ENCOUNTER_SAVE", + "FRAGMENTATION_SAVE") + + # This list defines features that are booleans, so must be converted to/from string when writing/reading from file bool_param = ["RESTART", @@ -51,7 +55,10 @@ float_param = ["T0", "TSTART", "TSTOP", "DT", "CHK_RMIN", "CHK_RMAX", "CHK_EJECT", "CHK_QMIN", "DU2M", "MU2KG", "TU2S", "MIN_GMFRAG", "GMTINY"] -upper_str_param = ["OUT_TYPE","OUT_FORM","OUT_STAT","IN_TYPE","IN_FORM"] +upper_str_param = ["OUT_TYPE","OUT_FORM","OUT_STAT","IN_TYPE","IN_FORM","ENCOUNTER_SAVE","FRAGMENTATION_SAVE", "CHK_QMIN_COORD"] +lower_str_param = ["NC_IN", "PL_IN", "TP_IN", "CB_IN", "CHK_QMIN_RANGE"] + +param_keys = ['! VERSION'] + int_param + float_param + upper_str_param + lower_str_param+ bool_param # This defines Xarray Dataset variables that are strings, which must be processed due to quirks in how NetCDF-Fortran # handles strings differently than Python's Xarray. @@ -417,52 +424,22 @@ def write_labeled_param(param, param_file_name): Prints a text file containing the parameter information. """ outfile = open(param_file_name, 'w') - keylist = ['! VERSION', - 'T0', - 'TSTART', - 'TSTOP', - 'DT', - 'ISTEP_OUT', - 'DUMP_CADENCE', - 'NC_IN', - 'PL_IN', - 'TP_IN', - 'CB_IN', - 'IN_TYPE', - 'IN_FORM', - 'BIN_OUT', - 'OUT_FORM', - 'OUT_TYPE', - 'OUT_STAT', - 'CHK_QMIN', - 'CHK_RMIN', - 'CHK_RMAX', - 'CHK_EJECT', - 'CHK_QMIN_COORD', - 'CHK_QMIN_RANGE', - 'MU2KG', - 'TU2S', - 'DU2M', - 'GMTINY', - 'FRAGMENTATION', - 'MIN_GMFRAG', - 'RESTART'] ptmp = param.copy() # Print the list of key/value pairs in the preferred order - for key in keylist: + for key in param_keys: val = ptmp.pop(key, None) if val is not None: if type(val) is bool: - print(f"{key:<16} {bool2yesno(val)}", file=outfile) + print(f"{key:<32} {bool2yesno(val)}", file=outfile) else: - print(f"{key:<16} {val}", file=outfile) + print(f"{key:<32} {val}", file=outfile) # Print the remaining key/value pairs in whatever order for key, val in ptmp.items(): if val != "": if type(val) is bool: - print(f"{key:<16} {bool2yesno(val)}", file=outfile) + print(f"{key:<32} {bool2yesno(val)}", file=outfile) else: - print(f"{key:<16} {val}", file=outfile) + print(f"{key:<32} {val}", file=outfile) outfile.close() return @@ -823,6 +800,36 @@ def swifter2xr(param, verbose=True): if verbose: print(f"Successfully converted {ds.sizes['time']} output frames.") return ds +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 + ---------- + ds : Xarray dataset + + Returns + ------- + ds : xarray dataset + """ + + if param['OUT_TYPE'] == "NETCDF_DOUBLE": + ds = fix_types(ds,ftype=np.float64) + elif param['OUT_TYPE'] == "NETCDF_FLOAT": + ds = fix_types(ds,ftype=np.float32) + ds = ds.where(ds.id >=0 ,drop=True) + # Check if the name variable contains unique values. If so, make name the dimension instead of id + if "id" in ds.dims: + if len(np.unique(ds['name'])) == len(ds['name']): + ds = ds.swap_dims({"id" : "name"}) + ds = ds.reset_coords("id") + + return ds + def swiftest2xr(param, verbose=True): """ Converts a Swiftest binary data file into an xarray DataSet. @@ -837,17 +844,11 @@ def swiftest2xr(param, verbose=True): xarray dataset """ + if ((param['OUT_TYPE'] == 'NETCDF_DOUBLE') or (param['OUT_TYPE'] == 'NETCDF_FLOAT')): if verbose: print('\nCreating Dataset from NetCDF file') ds = xr.open_dataset(param['BIN_OUT'], mask_and_scale=False) - if param['OUT_TYPE'] == "NETCDF_DOUBLE": - ds = fix_types(ds,ftype=np.float64) - elif param['OUT_TYPE'] == "NETCDF_FLOAT": - ds = fix_types(ds,ftype=np.float32) - # Check if the name variable contains unique values. If so, make name the dimension instead of id - if len(np.unique(ds['name'])) == len(ds['name']): - ds = ds.swap_dims({"id" : "name"}) - ds = ds.reset_coords("id") + ds = process_netcdf_input(ds, param) else: print(f"Error encountered. OUT_TYPE {param['OUT_TYPE']} not recognized.") return None @@ -855,14 +856,14 @@ def swiftest2xr(param, verbose=True): return ds -def xstrip(a): +def xstrip_nonstr(a): """ Cleans up the string values in the DataSet to remove extra white space Parameters ---------- a : xarray dataset - + Returns ------- da : xarray dataset with the strings cleaned up @@ -870,6 +871,22 @@ def xstrip(a): func = lambda x: np.char.strip(x) return xr.apply_ufunc(func, a.str.decode(encoding='utf-8'),dask='parallelized') +def xstrip_str(a): + """ + Cleans up the string values in the DataSet to remove extra white space + + Parameters + ---------- + a : xarray dataset + + Returns + ------- + da : xarray dataset with the strings cleaned up + """ + func = lambda x: np.char.strip(x) + return xr.apply_ufunc(func, a,dask='parallelized') + + def string_converter(da): """ Converts a string to a unicode string @@ -883,9 +900,10 @@ def string_converter(da): da : xarray dataset with the strings cleaned up """ if da.dtype == np.dtype(object): - da = da.astype(' 0 and "param_file" not in kwargs or len(kwargs) > 1 and "param_file" in kwargs: + self.set_parameter(verbose=False, **kwargs) # Let the user know that there was a problem reading an old parameter file and we're going to create a new one if read_param and not param_file_found: @@ -480,6 +492,9 @@ def run(self,**kwargs): warnings.warn(msg,stacklevel=2) return + if not self.restart: + self.clean() + print(f"Running a {self.codename} {self.integrator} run from tstart={self.param['TSTART']} {self.TU_name} to tstop={self.param['TSTOP']} {self.TU_name}") self._run_swiftest_driver() @@ -775,6 +790,8 @@ def set_parameter(self, verbose: bool = True, **kwargs): "encounter_check_loops": "TRIANGULAR", "ephemeris_date": "MBCL", "restart": False, + "encounter_save" : "NONE", + "fragmentation_save" : "NONE" } param_file = kwargs.pop("param_file",None) @@ -1010,6 +1027,8 @@ def set_feature(self, tides: bool | None = None, interaction_loops: Literal["TRIANGULAR", "FLAT", "ADAPTIVE"] | None = None, encounter_check_loops: Literal["TRIANGULAR", "SORTSWEEP", "ADAPTIVE"] | None = None, + encounter_save: Literal["NONE", "TRAJECTORY", "CLOSEST"] | None = None, + fragmentation_save: Literal["NONE", "TRAJECTORY", "CLOSEST"] | None = None, verbose: bool | None = None, **kwargs: Any ): @@ -1021,11 +1040,22 @@ def set_feature(self, close_encounter_check : bool, optional Check for close encounters between bodies. If set to True, then the radii of massive bodies must be included in initial conditions. + encounter_save : {"NONE","TRAJECTORY","CLOSEST"}, default "NONE" + Indicate if and how encounter data should be saved. If set to "TRAJECTORY" the full close encounter + trajectories are saved to file. If set to "CLOSEST" only the trajectories at the time of closest approach + are saved. If set to "NONE" no trajectory information is saved. + *WARNING*: Enabling this feature could lead to very large files. general_relativity : bool, optional Include the post-Newtonian correction in acceleration calculations. fragmentation : bool, optional If set to True, this turns on the Fraggle fragment generation code and `rotation` must also be True. This argument only applies to Swiftest-SyMBA simulations. It will be ignored otherwise. + fragmentation_save : {"NONE","TRAJECTORY","CLOSEST"}, default "NONE" + Indicate if and how fragmentation data should be saved. If set to "TRAJECTORY" the full close encounter + trajectories associated with each collision are saved to file. If set to "CLOSEST" only the trajectories + at a the time the collision occurs are saved. If set to "NONE" no trajectory information is saved (collision + details are still logged fraggle.log). + *WARNING*: Enabling this feature could lead to very large files. minimum_fragment_gmass : float, optional If fragmentation is turned on, this sets the mimimum G*mass of a collisional fragment that can be generated. *Note.* Only set one of minimum_fragment_gmass or minimum_fragment_mass @@ -1179,6 +1209,33 @@ def set_feature(self, self.param["ENCOUNTER_CHECK"] = encounter_check_loops update_list.append("encounter_check_loops") + if encounter_save is not None: + valid_vals = ["NONE", "TRAJECTORY", "CLOSEST"] + encounter_save = encounter_save.upper() + if encounter_save not in valid_vals: + msg = f"{encounter_save} is not a valid option for encounter_save." + msg += f"\nMust be one of {valid_vals}" + warnings.warn(msg,stacklevel=2) + if "ENCOUNTER_SAVE" not in self.param: + self.param["ENCOUNTER_SAVE"] = valid_vals[0] + else: + self.param["ENCOUNTER_SAVE"] = encounter_save + update_list.append("encounter_save") + + + if fragmentation_save is not None: + fragmentation_save = fragmentation_save.upper() + valid_vals = ["NONE", "TRAJECTORY", "CLOSEST"] + if fragmentation_save not in valid_vals: + msg = f"{fragmentation_save} is not a valid option for fragmentation_save." + msg += f"\nMust be one of {valid_vals}" + warnings.warn(msg,stacklevel=2) + if "FRAGMENTATION_SAVE" not in self.param: + self.param["FRAGMENTATION_SAVE"] = valid_vals[0] + else: + self.param["FRAGMENTATION_SAVE"] = fragmentation_save + update_list.append("fragmentation_save") + self.param["TIDES"] = False feature_dict = self.get_feature(update_list, verbose) @@ -1210,6 +1267,8 @@ def get_feature(self, arg_list: str | List[str] | None = None, verbose: bool | N valid_var = {"close_encounter_check": "CHK_CLOSE", "fragmentation": "FRAGMENTATION", + "encounter_save": "ENCOUNTER_SAVE", + "fragmentation_save": "FRAGMENTATION_SAVE", "minimum_fragment_gmass": "MIN_GMFRAG", "rotation": "ROTATION", "general_relativity": "GR", @@ -2670,17 +2729,27 @@ def read_output_file(self,read_init_cond : bool = True): # Make a temporary copy of the parameter dictionary so we can supply the absolute path of the binary file # This is done to handle cases where the method is called from a different working directory than the simulation # results + + if "ENCOUNTER_SAVE" in self.param or "FRAGMENTATION_SAVE" in self.param: + read_encounter = self.param["ENCOUNTER_SAVE"] != "NONE" or self.param["FRAGMENTATION_SAVE"] != "NONE" + else: + read_encounter = False param_tmp = self.param.copy() param_tmp['BIN_OUT'] = os.path.join(self.simdir, self.param['BIN_OUT']) if self.codename == "Swiftest": self.data = io.swiftest2xr(param_tmp, verbose=self.verbose) if self.verbose: print('Swiftest simulation data stored as xarray DataSet .data') if read_init_cond: + if self.verbose: + print("Reading initial conditions file as .ic") if "NETCDF" in self.param['IN_TYPE']: - param_tmp['BIN_OUT'] = os.path.join(self.simdir, self.param['NC_IN']) + param_tmp['BIN_OUT'] = self.simdir / self.param['NC_IN'] + self.ic = io.swiftest2xr(param_tmp, verbose=self.verbose) else: self.ic = self.data.isel(time=0) + if read_encounter: + self.read_encounter() elif self.codename == "Swifter": self.data = io.swifter2xr(param_tmp, verbose=self.verbose) @@ -2691,6 +2760,22 @@ 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_encounter(self): + if self.verbose: + print("Reading encounter history file as .enc") + enc_files = self.simdir.glob("**/encounter_*.nc") + + # 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.enc = xr.open_mfdataset(enc_files,parallel=True,preprocess=partial_func,mask_and_scale=True) + self.enc = io.process_netcdf_input(self.enc, self.param) + + return + + def follow(self, codestyle="Swifter"): """ An implementation of the Swift tool_follow algorithm. Under development. Currently only for Swift simulations. @@ -2849,3 +2934,25 @@ def initial_conditions_from_bin(self, framenum=-1, new_param=None, new_param_fil self.write_param(new_param_file, param=new_param) return frame + + def clean(self): + """ + Cleans up simulation directory by deleting old files (dump, logs, and output files). + """ + old_files = [self.simdir / self.param['BIN_OUT'], + self.simdir / "fraggle.log", + self.simdir / "swiftest.log", + ] + glob_files = [self.simdir.glob("**/dump_param?.in")] \ + + [self.simdir.glob("**/dump_bin?.nc")] \ + + [self.simdir.glob("**/enc*.nc")] \ + + [self.simdir.glob("**/frag*.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/CMakeLists.txt b/src/CMakeLists.txt index 668b07c47..7f45ddd5a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -28,7 +28,6 @@ SET(FAST_MATH_FILES ${SRC}/discard/discard.f90 ${SRC}/drift/drift.f90 ${SRC}/encounter/encounter_check.f90 - ${SRC}/encounter/encounter_io.f90 ${SRC}/encounter/encounter_setup.f90 ${SRC}/encounter/encounter_util.f90 ${SRC}/fraggle/fraggle_generate.f90 @@ -81,6 +80,7 @@ SET(FAST_MATH_FILES ${SRC}/util/util_minimize_bfgs.f90 ${SRC}/util/util_peri.f90 ${SRC}/util/util_rescale.f90 + ${SRC}/util/util_reset.f90 ${SRC}/util/util_resize.f90 ${SRC}/util/util_set.f90 ${SRC}/util/util_solve.f90 diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index 41ece554b..1e0be68bd 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -38,8 +38,8 @@ module subroutine discard_system(self, param) call tp%discard(system, param) ltp_discards = (tp_discards%nbody > 0) end if - if (ltp_discards) call tp_discards%write_info(param%nciu, param) - if (lpl_discards) call pl_discards%write_info(param%nciu, param) + if (ltp_discards) call tp_discards%write_info(param%nc, param) + if (lpl_discards) call pl_discards%write_info(param%nc, param) if (lpl_discards .and. param%lenergy) call self%conservation_report(param, lterminal=.false.) if (lpl_check) call pl_discards%setup(0,param) if (ltp_check) call tp_discards%setup(0,param) diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 deleted file mode 100644 index f41c25f85..000000000 --- a/src/encounter/encounter_io.f90 +++ /dev/null @@ -1,172 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -submodule (encounter_classes) s_encounter_io - use swiftest - use netcdf -contains - - module subroutine encounter_io_dump_storage_list(self, param) - !! author: David A. Minton - !! - !! Dumps the time history of an encounter to file. - implicit none - ! Arguments - class(encounter_storage(*)), intent(inout) :: self !! Encounter storage object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: i - - ! Most of this is just temporary test code just to get something working. Eventually this should get cleaned up. - call self%nciu%initialize(param) - do i = 1, self%nframes - if (allocated(self%frame(i)%item)) then - select type(plplenc_list => self%frame(i)%item) - class is (symba_plplenc) - self%nciu%ienc_frame = i - call plplenc_list%write_frame(self%nciu,param) - end select - end if - end do - call self%nciu%close() - - - return - end subroutine encounter_io_dump_storage_list - - - module subroutine encounter_io_initialize_output(self, param) - !! author: David A. Minton - !! - !! Initialize a NetCDF encounter file system. This is a simplified version of the main simulation output NetCDF file, but with fewer variables. - use, intrinsic :: ieee_arithmetic - implicit none - ! Arguments - class(encounter_io_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: nvar, varid, vartype - real(DP) :: dfill - real(SP) :: sfill - logical :: fileExists - character(len=STRMAX) :: errmsg - integer(I4B) :: ndims - - - associate(nciu => self) - dfill = ieee_value(dfill, IEEE_QUIET_NAN) - sfill = ieee_value(sfill, IEEE_QUIET_NAN) - - select case (param%out_type) - case("NETCDF_FLOAT") - self%out_type = NF90_FLOAT - case("NETCDF_DOUBLE") - self%out_type = NF90_DOUBLE - end select - - - ! Check if the file exists, and if it does, delete it - inquire(file=nciu%enc_file, exist=fileExists) - if (fileExists) then - open(unit=LUN, file=nciu%enc_file, status="old", err=667, iomsg=errmsg) - close(unit=LUN, status="delete") - end if - - call check( nf90_create(nciu%enc_file, NF90_NETCDF4, nciu%id), "encounter_io_initialize_output nf90_create" ) - - ! Dimensions - call check( nf90_def_dim(nciu%id, nciu%time_dimname, NF90_UNLIMITED, nciu%time_dimid), "encounter_io_initialize_output nf90_def_dim time_dimid" ) ! Simulation time dimension - call check( nf90_def_dim(nciu%id, nciu%space_dimname, NDIM, nciu%space_dimid), "encounter_io_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension - call check( nf90_def_dim(nciu%id, nciu%id_dimname, NF90_UNLIMITED, nciu%id_dimid), "encounter_io_initialize_output nf90_def_dim id_dimid" ) ! dimension to store particle id numbers - call check( nf90_def_dim(nciu%id, nciu%str_dimname, NAMELEN, nciu%str_dimid), "encounter_io_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) - - ! Dimension coordinates - call check( nf90_def_var(nciu%id, nciu%time_dimname, nciu%out_type, nciu%time_dimid, nciu%time_varid), "encounter_io_initialize_output nf90_def_var time_varid" ) - call check( nf90_def_var(nciu%id, nciu%space_dimname, NF90_CHAR, nciu%space_dimid, nciu%space_varid), "encounter_io_initialize_output nf90_def_var space_varid" ) - call check( nf90_def_var(nciu%id, nciu%id_dimname, NF90_INT, nciu%id_dimid, nciu%id_varid), "encounter_io_initialize_output nf90_def_var id_varid" ) - call check( nf90_def_var(nciu%id, nciu%name_varname, NF90_CHAR, [nciu%str_dimid, nciu%id_dimid], nciu%name_varid), "encounter_io_initialize_output nf90_def_var name_varid" ) - - ! Variables - call check( nf90_def_var(nciu%id, nciu%name_varname, NF90_CHAR, [nciu%str_dimid, nciu%id_dimid], nciu%name_varid), "encounter_io_initialize_output nf90_def_var name_varid" ) - call check( nf90_def_var(nciu%id, nciu%ptype_varname, NF90_CHAR, [nciu%str_dimid, nciu%id_dimid], nciu%ptype_varid), "encounter_io_initialize_output nf90_def_var ptype_varid" ) - call check( nf90_def_var(nciu%id, nciu%rh_varname, nciu%out_type, [nciu%space_dimid, nciu%id_dimid, nciu%time_dimid], nciu%rh_varid), "encounter_io_initialize_output nf90_def_var rh_varid" ) - call check( nf90_def_var(nciu%id, nciu%vh_varname, nciu%out_type, [nciu%space_dimid, nciu%id_dimid, nciu%time_dimid], nciu%vh_varid), "encounter_io_initialize_output nf90_def_var vh_varid" ) - call check( nf90_def_var(nciu%id, nciu%gmass_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%Gmass_varid), "encounter_io_initialize_output nf90_def_var Gmass_varid" ) - if (param%lclose) then - call check( nf90_def_var(nciu%id, nciu%radius_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%radius_varid), "encounter_io_initialize_output nf90_def_var radius_varid" ) - end if - if (param%lrotation) then - call check( nf90_def_var(nciu%id, nciu%Ip_varname, nciu%out_type, [nciu%space_dimid, nciu%id_dimid, nciu%time_dimid], nciu%Ip_varid), "encounter_io_initialize_output nf90_def_var Ip_varid" ) - call check( nf90_def_var(nciu%id, nciu%rot_varname, nciu%out_type, [nciu%space_dimid, nciu%id_dimid, nciu%time_dimid], nciu%rot_varid), "encounter_io_initialize_output nf90_def_var rot_varid" ) - end if - - call check( nf90_inquire(nciu%id, nVariables=nvar), "encounter_io_initialize_output nf90_inquire nVariables" ) - do varid = 1, nvar - call check( nf90_inquire_variable(nciu%id, varid, xtype=vartype, ndims=ndims), "encounter_io_initialize_output nf90_inquire_variable" ) - select case(vartype) - case(NF90_INT) - call check( nf90_def_var_fill(nciu%id, varid, 0, NF90_FILL_INT), "encounter_io_initialize_output nf90_def_var_fill NF90_INT" ) - case(NF90_FLOAT) - call check( nf90_def_var_fill(nciu%id, varid, 0, sfill), "encounter_io_initialize_output nf90_def_var_fill NF90_FLOAT" ) - case(NF90_DOUBLE) - call check( nf90_def_var_fill(nciu%id, varid, 0, dfill), "encounter_io_initialize_output nf90_def_var_fill NF90_DOUBLE" ) - case(NF90_CHAR) - call check( nf90_def_var_fill(nciu%id, varid, 0, 0), "encounter_io_initialize_output nf90_def_var_fill NF90_CHAR" ) - end select - end do - - ! Take the file out of define mode - call check( nf90_enddef(nciu%id), "encounter_io_initialize_output nf90_enddef" ) - end associate - - return - - 667 continue - write(*,*) "Error creating encounter output file. " // trim(adjustl(errmsg)) - call util_exit(FAILURE) - end subroutine encounter_io_initialize_output - - - module subroutine encounter_io_write_frame(self, nciu, param) - !! author: David A. Minton - !! - !! Write a frame of output of an encounter list structure. - implicit none - ! Arguments - class(encounter_list), intent(in) :: self !! Swiftest encounter structure - class(encounter_io_parameters), intent(inout) :: nciu !! Parameters used to identify a particular encounter io NetCDF dataset - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: tslot,i,old_mode, n - character(len=NAMELEN) :: charstring - - tslot = nciu%ienc_frame - call check( nf90_set_fill(nciu%id, nf90_nofill, old_mode), "encounter_io_write_frame_base nf90_set_fill" ) - call check( nf90_put_var(nciu%id, nciu%time_varid, self%t, start=[tslot]), "encounter_io_write_frame nf90_put_var time_varid" ) - - ! charstring = trim(adjustl(self%info(j)%name)) - ! call check( nf90_put_var(nciu%id, nciu%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "encounter_io_write_info_base nf90_put_var name_varid" ) - - ! charstring = trim(adjustl(self%info(j)%particle_type)) - ! call check( nf90_put_var(nciu%id, nciu%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "encounter_io_write_info_base nf90_put_var particle_type_varid" ) - - - ! call check( nf90_put_var(nciu%id, nciu%rh_varid, self%rh(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame_base nf90_put_var rh_varid" ) - ! call check( nf90_put_var(nciu%id, nciu%vh_varid, self%vh(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame_base nf90_put_var vh_varid" ) - ! call check( nf90_put_var(nciu%id, nciu%Gmass_varid, self%Gmass(j), start=[idslot, tslot]), "encounter_io_write_frame_base nf90_put_var body Gmass_varid" ) - ! if (param%lclose) call check( nf90_put_var(nciu%id, nciu%radius_varid, self%radius(j), start=[idslot, tslot]), "encounter_io_write_frame_base nf90_put_var body radius_varid" ) - ! if (param%lrotation) then - ! call check( nf90_put_var(nciu%id, nciu%Ip_varid, self%Ip(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame_base nf90_put_var body Ip_varid" ) - ! call check( nf90_put_var(nciu%id, nciu%rot_varid, self%rot(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame_base nf90_put_var body rotx_varid" ) - ! end if - - return - end subroutine encounter_io_write_frame - -end submodule s_encounter_io \ No newline at end of file diff --git a/src/io/io.f90 b/src/io/io.f90 index 0f8c6761c..92924d458 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -179,8 +179,8 @@ module subroutine io_conservation_report(self, param, lterminal) 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(param%nciu, param) - call param%nciu%close() + call self%write_frame(param%nc, param) + call param%nc%close() call util_exit(FAILURE) end if end if @@ -248,8 +248,8 @@ module subroutine io_dump_system(self, param) dump_param%out_stat = 'APPEND' dump_param%in_type = "NETCDF_DOUBLE" dump_param%in_netcdf = trim(adjustl(DUMP_NC_FILE(idx))) - dump_param%nciu%id_chunk = self%pl%nbody + self%tp%nbody - dump_param%nciu%time_chunk = 1 + dump_param%nc%id_chunk = self%pl%nbody + self%tp%nbody + dump_param%nc%time_chunk = 1 dump_param%tstart = self%t call dump_param%dump(param_file_name) @@ -257,11 +257,11 @@ module subroutine io_dump_system(self, param) dump_param%out_form = "XV" dump_param%outfile = trim(adjustl(DUMP_NC_FILE(idx))) dump_param%ioutput = 1 - call dump_param%nciu%initialize(dump_param) - call self%write_frame(dump_param%nciu, dump_param) - call dump_param%nciu%close() + call dump_param%nc%initialize(dump_param) + call self%write_frame(dump_param%nc, dump_param) + call dump_param%nc%close() ! Syncrhonize the disk and memory buffer of the NetCDF file (e.g. commit the frame files stored in memory to disk) - call param%nciu%flush(param) + call param%nc%flush(param) idx = idx + 1 if (idx > NDUMPFILES) idx = 1 @@ -671,7 +671,8 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) else if (param_value == "YES" .or. param_value == 'T') then param%lrestart = .true. end if - case ("NPLMAX", "NTPMAX", "GMTINY", "MIN_GMFRAG", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters + ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters + case ("NPLMAX", "NTPMAX", "GMTINY", "MIN_GMFRAG", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP", "ENCOUNTER_SAVE", "FRAGMENTATION_SAVE") case default write(*,*) "Ignoring unknown parameter -> ",param_name end select @@ -1308,7 +1309,7 @@ module subroutine io_read_in_system(self, param) ! Turn off energy computation so we don't have to feed it into the initial conditions tmp_param%lenergy = .false. end if - ierr = self%read_frame(tmp_param%nciu, tmp_param) + ierr = self%read_frame(tmp_param%nc, tmp_param) deallocate(tmp_param) if (ierr /=0) call util_exit(FAILURE) end if @@ -1537,8 +1538,8 @@ module subroutine io_write_frame_system(self, param) character(len=STRMAX) :: errmsg logical :: fileExists - param%nciu%id_chunk = self%pl%nbody + self%tp%nbody - param%nciu%time_chunk = max(param%dump_cadence / param%istep_out, 1) + param%nc%id_chunk = self%pl%nbody + self%tp%nbody + param%nc%time_chunk = max(param%dump_cadence / param%istep_out, 1) if (lfirst) then inquire(file=param%outfile, exist=fileExists) @@ -1553,15 +1554,15 @@ module subroutine io_write_frame_system(self, param) errmsg = param%outfile // " Alread Exists! You must specify OUT_STAT = APPEND, REPLACE, or UNKNOWN" goto 667 end if - call param%nciu%initialize(param) + call param%nc%initialize(param) case('REPLACE', 'UNKNOWN') - call param%nciu%initialize(param) + call param%nc%initialize(param) end select lfirst = .false. end if - call self%write_frame(param%nciu, param) + call self%write_frame(param%nc, param) return diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index bc08a1a1b..4ff06d4a8 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -18,26 +18,26 @@ program swiftest_driver use swiftest implicit none - class(swiftest_nbody_system), allocatable :: nbody_system !! Polymorphic object containing the nbody system to be integrated - class(swiftest_parameters), allocatable :: param !! Run configuration parameters - character(len=:), allocatable :: integrator !! Integrator type code (see swiftest_globals for symbolic names) - character(len=:),allocatable :: param_file_name !! Name of the file containing user-defined parameters - character(len=:), allocatable :: display_style !! Style of the output display {"STANDARD", "COMPACT", "PROGRESS"}). Default is "STANDARD" - integer(I8B) :: istart !! Starting index for loop counter - integer(I8B) :: nloops !! Number of steps to take in the simulation - integer(I4B) :: iout !! Output cadence counter - integer(I4B) :: idump !! Dump cadence counter - type(walltimer) :: integration_timer !! Object used for computing elapsed wall time - real(DP) :: tfrac - type(progress_bar) :: pbar !! Object used to print out a progress bar - character(*), parameter :: statusfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // & - '"; Number of active pl, tp = ", I6, ", ", I6)' - character(*), parameter :: symbastatfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // & - '"; Number of active plm, pl, tp = ", I6, ", ", I6, ", ", I6)' - character(*), parameter :: pbarfmt = '("Time = ", ES12.5," of ",ES12.5)' - character(len=64) :: pbarmessage - - character(*), parameter :: symbacompactfmt = '(";NPLM",ES22.15,$)' + class(swiftest_nbody_system), allocatable :: nbody_system !! Polymorphic object containing the nbody system to be integrated + class(swiftest_parameters), allocatable :: param !! Run configuration parameters + character(len=:), allocatable :: integrator !! Integrator type code (see swiftest_globals for symbolic names) + character(len=:),allocatable :: param_file_name !! Name of the file containing user-defined parameters + character(len=:), allocatable :: display_style !! Style of the output display {"STANDARD", "COMPACT", "PROGRESS"}). Default is "STANDARD" + integer(I8B) :: istart !! Starting index for loop counter + integer(I8B) :: nloops !! Number of steps to take in the simulation + integer(I4B) :: iout !! Output cadence counter + integer(I4B) :: idump !! Dump cadence counter + type(walltimer) :: integration_timer !! Object used for computing elapsed wall time + real(DP) :: tfrac + type(progress_bar) :: pbar !! Object used to print out a progress bar + character(*), parameter :: statusfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // & + '"; Number of active pl, tp = ", I6, ", ", I6)' + character(*), parameter :: symbastatfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // & + '"; Number of active plm, pl, tp = ", I6, ", ", I6, ", ", I6)' + character(*), parameter :: pbarfmt = '("Time = ", ES12.5," of ",ES12.5)' + character(len=64) :: pbarmessage + + character(*), parameter :: symbacompactfmt = '(";NPLM",ES22.15,$)' type(swiftest_storage(nframes=:)), allocatable :: system_history call io_get_args(integrator, param_file_name, display_style) diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index 3b775370e..b04fca1a1 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -38,35 +38,9 @@ module encounter_classes procedure :: dealloc => encounter_util_dealloc_list !! Deallocates all allocatables procedure :: spill => encounter_util_spill_list !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) procedure :: resize => encounter_util_resize_list !! Checks the current size of the encounter list against the required size and extends it by a factor of 2 more than requested if it is too small. - procedure :: write_frame => encounter_io_write_frame !! Writes a frame of encounter data to file final :: encounter_util_final_list !! Finalize the encounter list - deallocates all allocatables end type encounter_list - !! NetCDF dimension and variable names for the enounter save object - type, extends(netcdf_parameters) :: encounter_io_parameters - integer(I4B) :: COLLIDER_DIM_SIZE = 2 !! Size of collider dimension - integer(I4B) :: ienc_frame !! Current frame number for the encounter history - character(STRMAX) :: enc_file = "encounter.nc" !! Encounter output file name - - character(NAMELEN) :: eid_dimname = "encounter" !! The index of the encountering pair in the encounter list - integer(I4B) :: eid_dimid !! ID for the encounter pair index dimension - character(NAMELEN) :: collider_dimname = "collider" !! Dimension that defines the colliding bodies (bodies 1 and 2 are at dimension coordinates 1 and 2, respectively) - integer(I4B) :: collider_dimid !! ID for the collider dimension - character(NAMELEN) :: nenc_varname = "nenc" !! Total number of encounters - integer(I4B) :: nenc_varid !! ID for the number of encounters variable - character(NAMELEN) :: level_varname = "level" !! Recursion depth - integer(I4B) :: level_varid !! ID for the recursion level variable - contains - procedure :: initialize => encounter_io_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object - end type encounter_io_parameters - - type, extends(swiftest_storage) :: encounter_storage - !! A class that that is used to store simulation history data between file output - type(encounter_io_parameters) :: nciu - contains - procedure :: dump => encounter_io_dump_storage_list !! Dumps contents of encounter history to file - end type encounter_storage - type encounter_bounding_box_1D integer(I4B) :: n !! Number of bodies with extents integer(I4B), dimension(:), allocatable :: ind !! Sorted minimum/maximum extent indices (value > n indicates an ending index) @@ -199,25 +173,6 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical array indicating which pairs are approaching end subroutine encounter_check_sweep_aabb_single_list - module subroutine encounter_io_dump_storage_list(self, param) - implicit none - class(encounter_storage(*)), intent(inout) :: self !! Encounter storage object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine encounter_io_dump_storage_list - - module subroutine encounter_io_initialize_output(self, param) - implicit none - class(encounter_io_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset - class(swiftest_parameters), intent(in) :: param - end subroutine encounter_io_initialize_output - - module subroutine encounter_io_write_frame(self, nciu, param) - implicit none - class(encounter_list), intent(in) :: self !! Swiftest encounter structure - class(encounter_io_parameters), intent(inout) :: nciu !! Parameters used to identify a particular encounter io NetCDF dataset - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine encounter_io_write_frame - module subroutine encounter_setup_aabb(self, n, n_last) implicit none class(encounter_bounding_box), intent(inout) :: self !! Swiftest encounter structure diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 4f87b7707..a3e109bb5 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -150,8 +150,8 @@ module swiftest_classes !> User defined parameters that are read in from the parameters input file. !> Each paramter is initialized to a default values. type :: swiftest_parameters - character(len=:), allocatable :: integrator !! Symbolic name of the nbody integrator used - character(len=:), allocatable :: param_file_name !! The name of the parameter file + character(len=:), allocatable :: integrator !! Symbolic 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 real(DP) :: t0 = 0.0_DP !! Integration reference time @@ -229,7 +229,7 @@ module swiftest_classes logical :: lyarkovsky = .false. !! Turn on Yarkovsky effect logical :: lyorp = .false. !! Turn on YORP effect - type(netcdf_parameters) :: nciu !! Object containing NetCDF parameters + type(netcdf_parameters) :: nc !! Object containing NetCDF parameters contains procedure :: reader => io_param_reader procedure :: writer => io_param_writer @@ -543,11 +543,13 @@ module swiftest_classes end type type :: swiftest_storage(nframes) - integer(I4B), len :: nframes - !! An abstract class that establishes the pattern for various storage objects - type(swiftest_storage_frame), dimension(nframes) :: frame + !! An class that establishes the pattern for various storage objects + integer(I4B), len :: nframes = 2048 !! Total number of frames that can be stored + type(swiftest_storage_frame), dimension(nframes) :: frame !! Array of stored frames + integer(I4B) :: iframe = 0 !! The current frame number contains - procedure :: dump => io_dump_storage + procedure :: dump => io_dump_storage !! Dumps storage object contents to file + procedure :: reset => util_reset_storage !! Resets a storage object by deallocating all items and resetting the frame counter to 0 end type swiftest_storage abstract interface @@ -1022,55 +1024,55 @@ module subroutine netcdf_sync(self) class(netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset end subroutine netcdf_sync - module function netcdf_read_frame_system(self, nciu, param) result(ierr) + module function netcdf_read_frame_system(self, nc, param) result(ierr) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(netcdf_parameters), intent(inout) :: nciu !! Parameters used to for reading a NetCDF dataset to file + class(netcdf_parameters), intent(inout) :: nc !! Parameters used to for reading a NetCDF dataset to file class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters integer(I4B) :: ierr !! Error code: returns 0 if the read is successful end function netcdf_read_frame_system - module subroutine netcdf_read_hdr_system(self, nciu, param) + module subroutine netcdf_read_hdr_system(self, nc, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object - class(netcdf_parameters), intent(inout) :: nciu !! Parameters used to for reading a NetCDF dataset to file + class(netcdf_parameters), intent(inout) :: nc !! Parameters used to for reading a NetCDF dataset to file class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine netcdf_read_hdr_system - module subroutine netcdf_read_particle_info_system(self, nciu, param, plmask, tpmask) + module subroutine netcdf_read_particle_info_system(self, nc, param, plmask, tpmask) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object - class(netcdf_parameters), intent(inout) :: nciu !! Parameters used to identify a particular NetCDF dataset + class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters logical, dimension(:), intent(in) :: plmask !! Logical array indicating which index values belong to massive bodies logical, dimension(:), intent(in) :: tpmask !! Logical array indicating which index values belong to test particles end subroutine netcdf_read_particle_info_system - module subroutine netcdf_write_frame_base(self, nciu, param) + module subroutine netcdf_write_frame_base(self, nc, param) implicit none class(swiftest_base), intent(in) :: self !! Swiftest base object - class(netcdf_parameters), intent(inout) :: nciu !! Parameters used to for writing a NetCDF dataset to file + class(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 end subroutine netcdf_write_frame_base - module subroutine netcdf_write_frame_system(self, nciu, param) + module subroutine netcdf_write_frame_system(self, nc, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(netcdf_parameters), intent(inout) :: nciu !! Parameters used to for writing a NetCDF dataset to file + class(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 end subroutine netcdf_write_frame_system - module subroutine netcdf_write_hdr_system(self, nciu, param) + module subroutine netcdf_write_hdr_system(self, nc, param) implicit none class(swiftest_nbody_system), intent(in) :: self !! Swiftest nbody system object - class(netcdf_parameters), intent(inout) :: nciu !! Parameters used to for writing a NetCDF dataset to file + class(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 end subroutine netcdf_write_hdr_system - module subroutine netcdf_write_info_base(self, nciu, param) + module subroutine netcdf_write_info_base(self, nc, param) implicit none class(swiftest_base), intent(in) :: self !! Swiftest particle object - class(netcdf_parameters), intent(inout) :: nciu !! Parameters used to identify a particular NetCDF dataset + class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine netcdf_write_info_base @@ -1525,13 +1527,17 @@ module subroutine util_peri_tp(self, system, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine util_peri_tp - module subroutine util_rescale_system(self, param, mscale, dscale, tscale) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters. Returns with new values of the scale vactors and GU real(DP), intent(in) :: mscale, dscale, tscale !! Scale factors for mass, distance, and time units, respectively. end subroutine util_rescale_system + + module subroutine util_reset_storage(self) + implicit none + class(swiftest_storage(*)), intent(inout) :: self !! Swiftest storage object + end subroutine util_reset_storage end interface diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 432126e4e..68b94c373 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -13,10 +13,10 @@ module symba_classes !! Definition of classes and methods specific to the SyMBA integrator !! Adapted from David E. Kaufmann's Swifter routine: module_symba.f90 use swiftest_globals - use swiftest_classes, only : swiftest_parameters, swiftest_base, swiftest_particle_info, netcdf_parameters + use swiftest_classes, only : swiftest_parameters, swiftest_base, swiftest_particle_info, swiftest_storage, netcdf_parameters use helio_classes, only : helio_cb, helio_pl, helio_tp, helio_nbody_system use fraggle_classes, only : fraggle_colliders, fraggle_fragments - use encounter_classes, only : encounter_list, encounter_storage + use encounter_classes, only : encounter_list implicit none public @@ -26,11 +26,13 @@ module symba_classes real(DP), private, parameter :: RSHELL = 0.48075_DP type, extends(swiftest_parameters) :: symba_parameters - real(DP) :: GMTINY = -1.0_DP !! Smallest G*mass that is fully gravitating - real(DP) :: min_GMfrag = -1.0_DP !! Smallest G*mass that can be produced in a fragmentation event - integer(I4B), dimension(:), allocatable :: seed !! Random seeds - logical :: lfragmentation = .false. !! Do fragmentation modeling instead of simple merger. - character(STRMAX) :: encounter_save = "NONE" !! Indicate how encounter and/or fragmentation data should be saved + real(DP) :: GMTINY = -1.0_DP !! Smallest G*mass that is fully gravitating + real(DP) :: min_GMfrag = -1.0_DP !! Smallest G*mass that can be produced in a fragmentation event + integer(I4B), dimension(:), allocatable :: seed !! Random seeds + logical :: lfragmentation = .false. !! Do fragmentation modeling instead of simple merger. + character(STRMAX) :: encounter_save = "NONE" !! Indicate if and how encounter data should be saved + character(STRMAX) :: fragmentation_save = "NONE" !! Indicate if and how fragmentation data should be saved + logical :: lencounter_save = .false. !! Turns on encounter saving contains procedure :: reader => symba_io_param_reader procedure :: writer => symba_io_param_writer @@ -177,24 +179,39 @@ module symba_classes procedure :: resolve_collision => symba_collision_resolve_plplenc !! Process the pl-pl collision list, then modifiy the massive bodies based on the outcome of the c end type symba_plplenc - type, extends(helio_nbody_system) :: symba_system_snapshot + + !! NetCDF dimension and variable names for the enounter save object + type, extends(netcdf_parameters) :: symba_io_encounter_parameters + integer(I4B) :: COLLIDER_DIM_SIZE = 2 !! Size of collider dimension + integer(I4B) :: ienc_frame = 1 !! Current frame number for the encounter history + character(STRMAX) :: enc_file = "encounter.nc" !! Encounter output file name + + character(NAMELEN) :: level_varname = "level" !! Recursion depth + integer(I4B) :: level_varid !! ID for the recursion level variable + contains + procedure :: initialize => symba_io_encounter_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object + end type symba_io_encounter_parameters + + type, extends(swiftest_storage) :: symba_encounter_storage + !! A class that that is used to store simulation history data between file output + type(symba_io_encounter_parameters) :: nc !! NetCDF parameter object + real(DP), dimension(nframes) :: tvals !! Stored time values for snapshots contains - procedure :: snapshot => symba_util_take_system_snapshot - final :: symba_util_final_snapshot - end type + procedure :: dump => symba_io_encounter_dump !! Dumps contents of encounter history to file + final :: symba_util_final_encounter_storage + end type symba_encounter_storage + !******************************************************************************************************************************** ! symba_nbody_system class definitions and method interfaces !******************************************************************************************************************************** type, extends(helio_nbody_system) :: symba_nbody_system - class(symba_merger), allocatable :: pl_adds !! List of added bodies in mergers or collisions - class(symba_pltpenc), allocatable :: pltpenc_list !! List of massive body-test particle encounters in a single step - class(symba_plplenc), allocatable :: plplenc_list !! List of massive body-massive body encounters in a single step - class(symba_plplenc), allocatable :: plplcollision_list !! List of massive body-massive body collisions in a single step - integer(I4B) :: irec !! System recursion level - type(encounter_storage(nframes=:)), allocatable :: encounter_history !! Stores encounter history for later retrieval and saving to file - integer(I4B) :: ienc_frame = 0 !! Encounter history frame number - type(symba_system_snapshot) :: snapshot + class(symba_merger), allocatable :: pl_adds !! List of added bodies in mergers or collisions + class(symba_pltpenc), allocatable :: pltpenc_list !! List of massive body-test particle encounters in a single step + class(symba_plplenc), allocatable :: plplenc_list !! List of massive body-massive body encounters in a single step + class(symba_plplenc), allocatable :: plplcollision_list !! List of massive body-massive body collisions in a single step + integer(I4B) :: irec !! System recursion level + type(symba_encounter_storage(nframes=:)), allocatable :: encounter_history !! Stores encounter history for later retrieval and saving to file contains procedure :: write_discard => symba_io_write_discard !! Write out information about discarded and merged planets and test particles in SyMBA procedure :: initialize => symba_setup_initialize_system !! Performs SyMBA-specific initilization steps @@ -204,11 +221,22 @@ module symba_classes procedure :: recursive_step => symba_step_recur_system !! Step interacting planets and active test particles ahead in democratic heliocentric coordinates at the current recursion level, if applicable, and descend to the next deeper level if necessary procedure :: reset => symba_step_reset_system !! Resets pl, tp,and encounter structures at the start of a new step procedure :: dealloc => symba_util_dealloc_system !! Deallocates all allocatable arrays - procedure :: resize_storage => symba_util_resize_storage + procedure :: resize_storage => symba_util_resize_storage !! Resizes the encounter history storage object so that it contains enough spaces for the number of snapshots needed + procedure :: snapshot => symba_util_take_encounter_snapshot !! Take a minimal snapshot of the system through an encounter + procedure :: start_encounter => symba_io_start_encounter !! Initializes the new encounter history + procedure :: stop_encounter => symba_io_stop_encounter !! Saves the encounter and/or fragmentation data to file(s) final :: symba_util_final_system !! Finalizes the SyMBA nbody system object - deallocates all allocatables end type symba_nbody_system + type, extends(symba_nbody_system) :: symba_encounter_snapshot + integer(I4B) :: tslot !! The index for the time array in the final NetCDF file + contains + procedure :: write_encounter_frame => symba_io_encounter_write_frame !! Writes a frame of encounter data to file + generic :: write_frame => write_encounter_frame + final :: symba_util_final_encounter_snapshot + end type symba_encounter_snapshot + interface module function symba_collision_check_encounter(self, system, param, t, dt, irec) result(lany_collision) @@ -383,14 +411,32 @@ module subroutine symba_util_set_renc(self, scale) integer(I4B), intent(in) :: scale !! Current recursion depth end subroutine symba_util_set_renc - module subroutine symba_util_take_system_snapshot(self, system, param, t) + module subroutine symba_util_take_encounter_snapshot(self, param, t) use swiftest_classes, only : swiftest_parameters implicit none - class(symba_system_snapshot), intent(inout) :: self !! SyMBA nbody system snapshot object - class(symba_nbody_system), intent(in) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! current time - end subroutine symba_util_take_system_snapshot + class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! current time + end subroutine symba_util_take_encounter_snapshot + + module subroutine symba_io_encounter_dump(self, param) + implicit none + class(symba_encounter_storage(*)), intent(inout) :: self !! Encounter storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine symba_io_encounter_dump + + module subroutine symba_io_encounter_initialize_output(self, param) + implicit none + class(symba_io_encounter_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(in) :: param + end subroutine symba_io_encounter_initialize_output + + module subroutine symba_io_encounter_write_frame(self, nc, param) + implicit none + class(symba_encounter_snapshot), intent(in) :: self !! Swiftest encounter structure + class(symba_io_encounter_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine symba_io_encounter_write_frame module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) implicit none @@ -414,6 +460,20 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 end subroutine symba_io_param_writer + module subroutine symba_io_start_encounter(self, param, t) + implicit none + class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current simulation time + end subroutine symba_io_start_encounter + + module subroutine symba_io_stop_encounter(self, param, t) + implicit none + class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current simulation time + end subroutine symba_io_stop_encounter + module subroutine symba_io_write_discard(self, param) use swiftest_classes, only : swiftest_parameters implicit none @@ -522,8 +582,8 @@ recursive module subroutine symba_step_recur_system(self, param, t, ireci) implicit none class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - real(DP), intent(in) :: t - integer(I4B), intent(in) :: ireci !! input recursion level + real(DP), intent(in) :: t !! Current simulation time + integer(I4B), intent(in) :: ireci !! input recursion level end subroutine symba_step_recur_system module subroutine symba_step_reset_system(self, param) @@ -651,6 +711,16 @@ module subroutine symba_util_final_encounter_list(self) type(symba_encounter), intent(inout) :: self !! SyMBA encounter list object end subroutine symba_util_final_encounter_list + module subroutine symba_util_final_encounter_snapshot(self) + implicit none + type(symba_encounter_snapshot), intent(inout) :: self !! SyMBA nbody system object + end subroutine symba_util_final_encounter_snapshot + + module subroutine symba_util_final_encounter_storage(self) + implicit none + type(symba_encounter_storage(*)), intent(inout) :: self !! SyMBA nbody system object + end subroutine symba_util_final_encounter_storage + module subroutine symba_util_final_kin(self) implicit none type(symba_kinship), intent(inout) :: self !! SyMBA kinship object @@ -666,11 +736,6 @@ module subroutine symba_util_final_pl(self) type(symba_pl), intent(inout) :: self !! SyMBA massive body object end subroutine symba_util_final_pl - module subroutine symba_util_final_snapshot(self) - implicit none - type(symba_system_snapshot), intent(inout) :: self !! SyMBA nbody system object - end subroutine symba_util_final_snapshot - module subroutine symba_util_final_system(self) implicit none type(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index c7bd9a256..31de061de 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -80,38 +80,38 @@ module function netcdf_get_old_t_final_system(self, param) result(old_t_final) real(DP), dimension(NDIM) :: vectemp, rot0, Ip0, Lnow real(DP) :: KE_orb_orig, KE_spin_orig, PE_orig - call param%nciu%open(param) - call check( nf90_inquire_dimension(param%nciu%id, param%nciu%time_dimid, len=itmax), "netcdf_get_old_t_final_system time_dimid" ) - call check( nf90_inquire_dimension(param%nciu%id, param%nciu%id_dimid, len=idmax), "netcdf_get_old_t_final_system id_dimid" ) + call param%nc%open(param) + call check( nf90_inquire_dimension(param%nc%id, param%nc%time_dimid, len=itmax), "netcdf_get_old_t_final_system time_dimid" ) + call check( nf90_inquire_dimension(param%nc%id, param%nc%id_dimid, len=idmax), "netcdf_get_old_t_final_system id_dimid" ) allocate(vals(idmax)) - call check( nf90_get_var(param%nciu%id, param%nciu%time_varid, rtemp, start=[1], count=[1]), "netcdf_get_old_t_final_system time_varid" ) + call check( nf90_get_var(param%nc%id, param%nc%time_varid, rtemp, start=[1], count=[1]), "netcdf_get_old_t_final_system time_varid" ) !old_t_final = rtemp(1) old_t_final = param%t0 ! For NetCDF it is safe to overwrite the final t value on a restart if (param%lenergy) then - call check( nf90_get_var(param%nciu%id, param%nciu%KE_orb_varid, rtemp, start=[1], count=[1]), "netcdf_get_old_t_final_system KE_orb_varid" ) + call check( nf90_get_var(param%nc%id, param%nc%KE_orb_varid, rtemp, start=[1], count=[1]), "netcdf_get_old_t_final_system KE_orb_varid" ) KE_orb_orig = rtemp(1) - call check( nf90_get_var(param%nciu%id, param%nciu%KE_spin_varid, rtemp, start=[1], count=[1]), "netcdf_get_old_t_final_system KE_spin_varid" ) + call check( nf90_get_var(param%nc%id, param%nc%KE_spin_varid, rtemp, start=[1], count=[1]), "netcdf_get_old_t_final_system KE_spin_varid" ) KE_spin_orig = rtemp(1) - call check( nf90_get_var(param%nciu%id, param%nciu%PE_varid, rtemp, start=[1], count=[1]), "netcdf_get_old_t_final_system PE_varid" ) + call check( nf90_get_var(param%nc%id, param%nc%PE_varid, rtemp, start=[1], count=[1]), "netcdf_get_old_t_final_system PE_varid" ) PE_orig = rtemp(1) - call check( nf90_get_var(param%nciu%id, param%nciu%Ecollisions_varid, self%Ecollisions, start=[1]), "netcdf_get_old_t_final_system Ecollisions_varid" ) - call check( nf90_get_var(param%nciu%id, param%nciu%Euntracked_varid, self%Euntracked, start=[1]), "netcdf_get_old_t_final_system Euntracked_varid" ) + call check( nf90_get_var(param%nc%id, param%nc%Ecollisions_varid, self%Ecollisions, start=[1]), "netcdf_get_old_t_final_system Ecollisions_varid" ) + call check( nf90_get_var(param%nc%id, param%nc%Euntracked_varid, self%Euntracked, start=[1]), "netcdf_get_old_t_final_system Euntracked_varid" ) self%Eorbit_orig = KE_orb_orig + KE_spin_orig + PE_orig + self%Ecollisions + self%Euntracked - call check( nf90_get_var(param%nciu%id, param%nciu%L_orb_varid, self%Lorbit_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_get_old_t_final_system L_orb_varid" ) - call check( nf90_get_var(param%nciu%id, param%nciu%L_spin_varid, self%Lspin_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_get_old_t_final_system L_spin_varid" ) - call check( nf90_get_var(param%nciu%id, param%nciu%L_escape_varid, self%Lescape(:), start=[1,1], count=[NDIM,1]), "netcdf_get_old_t_final_system L_escape_varid" ) + call check( nf90_get_var(param%nc%id, param%nc%L_orb_varid, self%Lorbit_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_get_old_t_final_system L_orb_varid" ) + call check( nf90_get_var(param%nc%id, param%nc%L_spin_varid, self%Lspin_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_get_old_t_final_system L_spin_varid" ) + call check( nf90_get_var(param%nc%id, param%nc%L_escape_varid, self%Lescape(:), start=[1,1], count=[NDIM,1]), "netcdf_get_old_t_final_system L_escape_varid" ) self%Ltot_orig(:) = self%Lorbit_orig(:) + self%Lspin_orig(:) + self%Lescape(:) - call check( nf90_get_var(param%nciu%id, param%nciu%Gmass_varid, vals, start=[1,1], count=[idmax,1]), "netcdf_get_old_t_final_system Gmass_varid" ) - call check( nf90_get_var(param%nciu%id, param%nciu%GMescape_varid, self%GMescape, start=[1]), "netcdf_get_old_t_final_system GMescape_varid" ) + call check( nf90_get_var(param%nc%id, param%nc%Gmass_varid, vals, start=[1,1], count=[idmax,1]), "netcdf_get_old_t_final_system Gmass_varid" ) + call check( nf90_get_var(param%nc%id, param%nc%GMescape_varid, self%GMescape, start=[1]), "netcdf_get_old_t_final_system GMescape_varid" ) self%GMtot_orig = vals(1) + sum(vals(2:idmax), vals(2:idmax) == vals(2:idmax)) + self%GMescape select type(cb => self%cb) @@ -119,13 +119,13 @@ module function netcdf_get_old_t_final_system(self, param) result(old_t_final) cb%GM0 = vals(1) cb%dGM = cb%Gmass - cb%GM0 - call check( nf90_get_var(param%nciu%id, param%nciu%radius_varid, rtemp, start=[1,1], count=[1,1]), "netcdf_get_old_t_final_system radius_varid" ) + call check( nf90_get_var(param%nc%id, param%nc%radius_varid, rtemp, start=[1,1], count=[1,1]), "netcdf_get_old_t_final_system radius_varid" ) cb%R0 = rtemp(1) if (param%lrotation) then - call check( nf90_get_var(param%nciu%id, param%nciu%rot_varid, rot0, start=[1,1,1], count=[NDIM,1,1]), "netcdf_get_old_t_final_system rot_varid" ) - call check( nf90_get_var(param%nciu%id, param%nciu%Ip_varid, Ip0, start=[1,1,1], count=[NDIM,1,1]), "netcdf_get_old_t_final_system Ip_varid" ) + call check( nf90_get_var(param%nc%id, param%nc%rot_varid, rot0, start=[1,1,1], count=[NDIM,1,1]), "netcdf_get_old_t_final_system rot_varid" ) + call check( nf90_get_var(param%nc%id, param%nc%Ip_varid, Ip0, start=[1,1,1], count=[NDIM,1,1]), "netcdf_get_old_t_final_system Ip_varid" ) cb%L0(:) = Ip0(3) * cb%GM0 * cb%R0**2 * rot0(:) @@ -159,16 +159,16 @@ module subroutine netcdf_initialize_output(self, param) character(len=STRMAX) :: errmsg integer(I4B) :: ndims - associate(nciu => self) + associate(nc => self) dfill = ieee_value(dfill, IEEE_QUIET_NAN) sfill = ieee_value(sfill, IEEE_QUIET_NAN) select case (param%out_type) case("NETCDF_FLOAT") - nciu%out_type = NF90_FLOAT + nc%out_type = NF90_FLOAT case("NETCDF_DOUBLE") - nciu%out_type = NF90_DOUBLE + nc%out_type = NF90_DOUBLE end select ! Check if the file exists, and if it does, delete it @@ -179,120 +179,121 @@ module subroutine netcdf_initialize_output(self, param) end if ! Create the file - call check( nf90_create(param%outfile, NF90_NETCDF4, nciu%id), "netcdf_initialize_output nf90_create" ) + call check( nf90_create(param%outfile, NF90_NETCDF4, nc%id), "netcdf_initialize_output nf90_create" ) ! Dimensions - call check( nf90_def_dim(nciu%id, nciu%time_dimname, NF90_UNLIMITED, nciu%time_dimid), "netcdf_initialize_output nf90_def_dim time_dimid" ) ! Simulation time dimension - call check( nf90_def_dim(nciu%id, nciu%space_dimname, NDIM, nciu%space_dimid), "netcdf_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension - call check( nf90_def_dim(nciu%id, nciu%id_dimname, NF90_UNLIMITED, nciu%id_dimid), "netcdf_initialize_output nf90_def_dim id_dimid" ) ! dimension to store particle id numbers - call check( nf90_def_dim(nciu%id, nciu%str_dimname, NAMELEN, nciu%str_dimid), "netcdf_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) + call check( nf90_def_dim(nc%id, nc%time_dimname, NF90_UNLIMITED, nc%time_dimid), "netcdf_initialize_output nf90_def_dim time_dimid" ) ! Simulation time dimension + call check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "netcdf_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension + call check( nf90_def_dim(nc%id, nc%id_dimname, NF90_UNLIMITED, nc%id_dimid), "netcdf_initialize_output nf90_def_dim id_dimid" ) ! dimension to store particle id numbers + call check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "netcdf_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) ! Dimension coordinates - call check( nf90_def_var(nciu%id, nciu%time_dimname, nciu%out_type, nciu%time_dimid, nciu%time_varid), "netcdf_initialize_output nf90_def_var time_varid" ) - call check( nf90_def_var(nciu%id, nciu%space_dimname, NF90_CHAR, nciu%space_dimid, nciu%space_varid), "netcdf_initialize_output nf90_def_var space_varid" ) - call check( nf90_def_var(nciu%id, nciu%id_dimname, NF90_INT, nciu%id_dimid, nciu%id_varid), "netcdf_initialize_output nf90_def_var id_varid" ) + call check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, nc%time_dimid, nc%time_varid), "netcdf_initialize_output nf90_def_var time_varid" ) + call check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "netcdf_initialize_output nf90_def_var space_varid" ) + call check( nf90_def_var(nc%id, nc%id_dimname, NF90_INT, nc%id_dimid, nc%id_varid), "netcdf_initialize_output nf90_def_var id_varid" ) ! Variables - call check( nf90_def_var(nciu%id, nciu%npl_varname, NF90_INT, nciu%time_dimid, nciu%npl_varid), "netcdf_initialize_output nf90_def_var npl_varid" ) - call check( nf90_def_var(nciu%id, nciu%ntp_varname, NF90_INT, nciu%time_dimid, nciu%ntp_varid), "netcdf_initialize_output nf90_def_var ntp_varid" ) - if (param%integrator == SYMBA) call check( nf90_def_var(nciu%id, nciu%nplm_varname, NF90_INT, nciu%time_dimid, nciu%nplm_varid), "netcdf_initialize_output nf90_def_var nplm_varid" ) - call check( nf90_def_var(nciu%id, nciu%name_varname, NF90_CHAR, [nciu%str_dimid, nciu%id_dimid], nciu%name_varid), "netcdf_initialize_output nf90_def_var name_varid" ) - call check( nf90_def_var(nciu%id, nciu%ptype_varname, NF90_CHAR, [nciu%str_dimid, nciu%id_dimid], nciu%ptype_varid), "netcdf_initialize_output nf90_def_var ptype_varid" ) + call check( nf90_def_var(nc%id, nc%npl_varname, NF90_INT, nc%time_dimid, nc%npl_varid), "netcdf_initialize_output nf90_def_var npl_varid" ) + call check( nf90_def_var(nc%id, nc%ntp_varname, NF90_INT, nc%time_dimid, nc%ntp_varid), "netcdf_initialize_output nf90_def_var ntp_varid" ) + if (param%integrator == SYMBA) call check( nf90_def_var(nc%id, nc%nplm_varname, NF90_INT, nc%time_dimid, nc%nplm_varid), "netcdf_initialize_output nf90_def_var nplm_varid" ) + call check( nf90_def_var(nc%id, nc%name_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], nc%name_varid), "netcdf_initialize_output nf90_def_var name_varid" ) + call check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], nc%ptype_varid), "netcdf_initialize_output nf90_def_var ptype_varid" ) if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then - call check( nf90_def_var(nciu%id, nciu%rh_varname, nciu%out_type, [nciu%space_dimid, nciu%id_dimid, nciu%time_dimid], nciu%rh_varid), "netcdf_initialize_output nf90_def_var rh_varid" ) - call check( nf90_def_var(nciu%id, nciu%vh_varname, nciu%out_type, [nciu%space_dimid, nciu%id_dimid, nciu%time_dimid], nciu%vh_varid), "netcdf_initialize_output nf90_def_var vh_varid" ) + call check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%rh_varid), "netcdf_initialize_output nf90_def_var rh_varid" ) + call check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%vh_varid), "netcdf_initialize_output nf90_def_var vh_varid" ) !! When GR is enabled, we need to save the pseudovelocity vectors in addition to the true heliocentric velocity vectors, otherwise !! we cannnot expect bit-identical runs from restarted runs with GR enabled due to floating point errors during the conversion. if (param%lgr) then - call check( nf90_def_var(nciu%id, nciu%gr_pseudo_vh_varname, nciu%out_type, [nciu%space_dimid, nciu%id_dimid, nciu%time_dimid], nciu%gr_pseudo_vh_varid), "netcdf_initialize_output nf90_def_var gr_psuedo_vh_varid" ) - nciu%lpseudo_vel_exists = .true. + call check( nf90_def_var(nc%id, nc%gr_pseudo_vh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%gr_pseudo_vh_varid), "netcdf_initialize_output nf90_def_var gr_psuedo_vh_varid" ) + nc%lpseudo_vel_exists = .true. end if end if if ((param%out_form == "EL") .or. (param%out_form == "XVEL")) then - call check( nf90_def_var(nciu%id, nciu%a_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%a_varid), "netcdf_initialize_output nf90_def_var a_varid" ) - call check( nf90_def_var(nciu%id, nciu%e_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%e_varid), "netcdf_initialize_output nf90_def_var e_varid" ) - call check( nf90_def_var(nciu%id, nciu%inc_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%inc_varid), "netcdf_initialize_output nf90_def_var inc_varid" ) - call check( nf90_def_var(nciu%id, nciu%capom_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%capom_varid), "netcdf_initialize_output nf90_def_var capom_varid" ) - call check( nf90_def_var(nciu%id, nciu%omega_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%omega_varid), "netcdf_initialize_output nf90_def_var omega_varid" ) - call check( nf90_def_var(nciu%id, nciu%capm_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%capm_varid), "netcdf_initialize_output nf90_def_var capm_varid" ) - call check( nf90_def_var(nciu%id, nciu%varpi_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%varpi_varid), "netcdf_initialize_output nf90_def_var varpi_varid" ) - call check( nf90_def_var(nciu%id, nciu%lam_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%lam_varid), "netcdf_initialize_output nf90_def_var lam_varid" ) - call check( nf90_def_var(nciu%id, nciu%f_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%f_varid), "netcdf_initialize_output nf90_def_var f_varid" ) - call check( nf90_def_var(nciu%id, nciu%cape_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%cape_varid), "netcdf_initialize_output nf90_def_var cape_varid" ) + call check( nf90_def_var(nc%id, nc%a_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%a_varid), "netcdf_initialize_output nf90_def_var a_varid" ) + call check( nf90_def_var(nc%id, nc%e_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%e_varid), "netcdf_initialize_output nf90_def_var e_varid" ) + call check( nf90_def_var(nc%id, nc%inc_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%inc_varid), "netcdf_initialize_output nf90_def_var inc_varid" ) + call check( nf90_def_var(nc%id, nc%capom_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%capom_varid), "netcdf_initialize_output nf90_def_var capom_varid" ) + call check( nf90_def_var(nc%id, nc%omega_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%omega_varid), "netcdf_initialize_output nf90_def_var omega_varid" ) + call check( nf90_def_var(nc%id, nc%capm_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%capm_varid), "netcdf_initialize_output nf90_def_var capm_varid" ) + call check( nf90_def_var(nc%id, nc%varpi_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%varpi_varid), "netcdf_initialize_output nf90_def_var varpi_varid" ) + call check( nf90_def_var(nc%id, nc%lam_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%lam_varid), "netcdf_initialize_output nf90_def_var lam_varid" ) + call check( nf90_def_var(nc%id, nc%f_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%f_varid), "netcdf_initialize_output nf90_def_var f_varid" ) + call check( nf90_def_var(nc%id, nc%cape_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%cape_varid), "netcdf_initialize_output nf90_def_var cape_varid" ) end if - call check( nf90_def_var(nciu%id, nciu%gmass_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%Gmass_varid), "netcdf_initialize_output nf90_def_var Gmass_varid" ) + call check( nf90_def_var(nc%id, nc%gmass_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%Gmass_varid), "netcdf_initialize_output nf90_def_var Gmass_varid" ) if (param%lrhill_present) then - call check( nf90_def_var(nciu%id, nciu%rhill_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%rhill_varid), "netcdf_initialize_output nf90_def_var rhill_varid" ) + call check( nf90_def_var(nc%id, nc%rhill_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%rhill_varid), "netcdf_initialize_output nf90_def_var rhill_varid" ) end if if (param%lclose) then - call check( nf90_def_var(nciu%id, nciu%radius_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%radius_varid), "netcdf_initialize_output nf90_def_var radius_varid" ) - - call check( nf90_def_var(nciu%id, nciu%origin_time_varname, nciu%out_type, nciu%id_dimid, nciu%origin_time_varid), "netcdf_initialize_output nf90_def_var origin_time_varid" ) - call check( nf90_def_var(nciu%id, nciu%origin_type_varname, NF90_CHAR, [nciu%str_dimid, nciu%id_dimid], & - nciu%origin_type_varid), "netcdf_initialize_output nf90_create" ) - call check( nf90_def_var(nciu%id, nciu%origin_rh_varname, nciu%out_type, [nciu%space_dimid, nciu%id_dimid], nciu%origin_rh_varid), "netcdf_initialize_output nf90_def_var origin_rh_varid" ) - call check( nf90_def_var(nciu%id, nciu%origin_vh_varname, nciu%out_type, [nciu%space_dimid, nciu%id_dimid], nciu%origin_vh_varid), "netcdf_initialize_output nf90_def_var origin_vh_varid" ) - - call check( nf90_def_var(nciu%id, nciu%collision_id_varname, NF90_INT, nciu%id_dimid, nciu%collision_id_varid), "netcdf_initialize_output nf90_def_var collision_id_varid" ) - call check( nf90_def_var(nciu%id, nciu%discard_time_varname, nciu%out_type, nciu%id_dimid, nciu%discard_time_varid), "netcdf_initialize_output nf90_def_var discard_time_varid" ) - call check( nf90_def_var(nciu%id, nciu%discard_rh_varname, nciu%out_type, [nciu%space_dimid, nciu%id_dimid], nciu%discard_rh_varid), "netcdf_initialize_output nf90_def_var discard_rh_varid" ) - call check( nf90_def_var(nciu%id, nciu%discard_vh_varname, nciu%out_type, [nciu%space_dimid, nciu%id_dimid], nciu%discard_vh_varid), "netcdf_initialize_output nf90_def_var discard_vh_varid" ) - call check( nf90_def_var(nciu%id, nciu%discard_body_id_varname, NF90_INT, nciu%id_dimid, nciu%discard_body_id_varid), "netcdf_initialize_output nf90_def_var discard_body_id_varid" ) + call check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%radius_varid), "netcdf_initialize_output nf90_def_var radius_varid" ) + + call check( nf90_def_var(nc%id, nc%origin_time_varname, nc%out_type, nc%id_dimid, nc%origin_time_varid), "netcdf_initialize_output nf90_def_var origin_time_varid" ) + call check( nf90_def_var(nc%id, nc%origin_type_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], & + nc%origin_type_varid), "netcdf_initialize_output nf90_create" ) + call check( nf90_def_var(nc%id, nc%origin_rh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid], nc%origin_rh_varid), "netcdf_initialize_output nf90_def_var origin_rh_varid" ) + call check( nf90_def_var(nc%id, nc%origin_vh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid], nc%origin_vh_varid), "netcdf_initialize_output nf90_def_var origin_vh_varid" ) + + call check( nf90_def_var(nc%id, nc%collision_id_varname, NF90_INT, nc%id_dimid, nc%collision_id_varid), "netcdf_initialize_output nf90_def_var collision_id_varid" ) + call check( nf90_def_var(nc%id, nc%discard_time_varname, nc%out_type, nc%id_dimid, nc%discard_time_varid), "netcdf_initialize_output nf90_def_var discard_time_varid" ) + call check( nf90_def_var(nc%id, nc%discard_rh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid], nc%discard_rh_varid), "netcdf_initialize_output nf90_def_var discard_rh_varid" ) + call check( nf90_def_var(nc%id, nc%discard_vh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid], nc%discard_vh_varid), "netcdf_initialize_output nf90_def_var discard_vh_varid" ) + call check( nf90_def_var(nc%id, nc%discard_body_id_varname, NF90_INT, nc%id_dimid, nc%discard_body_id_varid), "netcdf_initialize_output nf90_def_var discard_body_id_varid" ) end if if (param%lrotation) then - call check( nf90_def_var(nciu%id, nciu%Ip_varname, nciu%out_type, [nciu%space_dimid, nciu%id_dimid, nciu%time_dimid], nciu%Ip_varid), "netcdf_initialize_output nf90_def_var Ip_varid" ) - call check( nf90_def_var(nciu%id, nciu%rot_varname, nciu%out_type, [nciu%space_dimid, nciu%id_dimid, nciu%time_dimid], nciu%rot_varid), "netcdf_initialize_output nf90_def_var rot_varid" ) + call check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%Ip_varid), "netcdf_initialize_output nf90_def_var Ip_varid" ) + call check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%rot_varid), "netcdf_initialize_output nf90_def_var rot_varid" ) end if ! if (param%ltides) then - ! call check( nf90_def_var(nciu%id, nciu%k2_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%k2_varid), "netcdf_initialize_output nf90_def_var k2_varid" ) - ! call check( nf90_def_var(nciu%id, nciu%q_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%Q_varid), "netcdf_initialize_output nf90_def_var Q_varid" ) + ! call check( nf90_def_var(nc%id, nc%k2_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%k2_varid), "netcdf_initialize_output nf90_def_var k2_varid" ) + ! call check( nf90_def_var(nc%id, nc%q_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%Q_varid), "netcdf_initialize_output nf90_def_var Q_varid" ) ! end if if (param%lenergy) then - call check( nf90_def_var(nciu%id, nciu%ke_orb_varname, nciu%out_type, nciu%time_dimid, nciu%KE_orb_varid), "netcdf_initialize_output nf90_def_var KE_orb_varid" ) - call check( nf90_def_var(nciu%id, nciu%ke_spin_varname, nciu%out_type, nciu%time_dimid, nciu%KE_spin_varid), "netcdf_initialize_output nf90_def_var KE_spin_varid" ) - call check( nf90_def_var(nciu%id, nciu%pe_varname, nciu%out_type, nciu%time_dimid, nciu%PE_varid), "netcdf_initialize_output nf90_def_var PE_varid" ) - call check( nf90_def_var(nciu%id, nciu%L_orb_varname, nciu%out_type, [nciu%space_dimid, nciu%time_dimid], nciu%L_orb_varid), "netcdf_initialize_output nf90_def_var L_orb_varid" ) - call check( nf90_def_var(nciu%id, nciu%L_spin_varname, nciu%out_type, [nciu%space_dimid, nciu%time_dimid], nciu%L_spin_varid), "netcdf_initialize_output nf90_def_var L_spin_varid" ) - call check( nf90_def_var(nciu%id, nciu%L_escape_varname, nciu%out_type, [nciu%space_dimid, nciu%time_dimid], nciu%L_escape_varid), "netcdf_initialize_output nf90_def_var L_escape_varid" ) - call check( nf90_def_var(nciu%id, nciu%Ecollisions_varname, nciu%out_type, nciu%time_dimid, nciu%Ecollisions_varid), "netcdf_initialize_output nf90_def_var Ecollisions_varid" ) - call check( nf90_def_var(nciu%id, nciu%Euntracked_varname, nciu%out_type, nciu%time_dimid, nciu%Euntracked_varid), "netcdf_initialize_output nf90_def_var Euntracked_varid" ) - call check( nf90_def_var(nciu%id, nciu%GMescape_varname, nciu%out_type, nciu%time_dimid, nciu%GMescape_varid), "netcdf_initialize_output nf90_def_var GMescape_varid" ) + call check( nf90_def_var(nc%id, nc%ke_orb_varname, nc%out_type, nc%time_dimid, nc%KE_orb_varid), "netcdf_initialize_output nf90_def_var KE_orb_varid" ) + call check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type, nc%time_dimid, nc%KE_spin_varid), "netcdf_initialize_output nf90_def_var KE_spin_varid" ) + call check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type, nc%time_dimid, nc%PE_varid), "netcdf_initialize_output nf90_def_var PE_varid" ) + call check( nf90_def_var(nc%id, nc%L_orb_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%L_orb_varid), "netcdf_initialize_output nf90_def_var L_orb_varid" ) + call check( nf90_def_var(nc%id, nc%L_spin_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%L_spin_varid), "netcdf_initialize_output nf90_def_var L_spin_varid" ) + call check( nf90_def_var(nc%id, nc%L_escape_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%L_escape_varid), "netcdf_initialize_output nf90_def_var L_escape_varid" ) + call check( nf90_def_var(nc%id, nc%Ecollisions_varname, nc%out_type, nc%time_dimid, nc%Ecollisions_varid), "netcdf_initialize_output nf90_def_var Ecollisions_varid" ) + call check( nf90_def_var(nc%id, nc%Euntracked_varname, nc%out_type, nc%time_dimid, nc%Euntracked_varid), "netcdf_initialize_output nf90_def_var Euntracked_varid" ) + call check( nf90_def_var(nc%id, nc%GMescape_varname, nc%out_type, nc%time_dimid, nc%GMescape_varid), "netcdf_initialize_output nf90_def_var GMescape_varid" ) end if - call check( nf90_def_var(nciu%id, nciu%j2rp2_varname, nciu%out_type, nciu%time_dimid, nciu%j2rp2_varid), "netcdf_initialize_output nf90_def_var j2rp2_varid" ) - call check( nf90_def_var(nciu%id, nciu%j4rp4_varname, nciu%out_type, nciu%time_dimid, nciu%j4rp4_varid), "netcdf_initialize_output nf90_def_var j4rp4_varid" ) + call check( nf90_def_var(nc%id, nc%j2rp2_varname, nc%out_type, nc%time_dimid, nc%j2rp2_varid), "netcdf_initialize_output nf90_def_var j2rp2_varid" ) + call check( nf90_def_var(nc%id, nc%j4rp4_varname, nc%out_type, nc%time_dimid, nc%j4rp4_varid), "netcdf_initialize_output nf90_def_var j4rp4_varid" ) ! Set fill mode to NaN for all variables - call check( nf90_inquire(nciu%id, nVariables=nvar), "netcdf_initialize_output nf90_inquire nVariables" ) + call check( nf90_inquire(nc%id, nVariables=nvar), "netcdf_initialize_output nf90_inquire nVariables" ) do varid = 1, nvar - call check( nf90_inquire_variable(nciu%id, varid, xtype=vartype, ndims=ndims), "netcdf_initialize_output nf90_inquire_variable" ) + call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "netcdf_initialize_output nf90_inquire_variable" ) select case(vartype) case(NF90_INT) - call check( nf90_def_var_fill(nciu%id, varid, 0, NF90_FILL_INT), "netcdf_initialize_output nf90_def_var_fill NF90_INT" ) + call check( nf90_def_var_fill(nc%id, varid, 0, NF90_FILL_INT), "netcdf_initialize_output nf90_def_var_fill NF90_INT" ) case(NF90_FLOAT) - call check( nf90_def_var_fill(nciu%id, varid, 0, sfill), "netcdf_initialize_output nf90_def_var_fill NF90_FLOAT" ) + call check( nf90_def_var_fill(nc%id, varid, 0, sfill), "netcdf_initialize_output nf90_def_var_fill NF90_FLOAT" ) case(NF90_DOUBLE) - call check( nf90_def_var_fill(nciu%id, varid, 0, dfill), "netcdf_initialize_output nf90_def_var_fill NF90_DOUBLE" ) + call check( nf90_def_var_fill(nc%id, varid, 0, dfill), "netcdf_initialize_output nf90_def_var_fill NF90_DOUBLE" ) case(NF90_CHAR) - call check( nf90_def_var_fill(nciu%id, varid, 0, 0), "netcdf_initialize_output nf90_def_var_fill NF90_CHAR" ) + call check( nf90_def_var_fill(nc%id, varid, 0, 0), "netcdf_initialize_output nf90_def_var_fill NF90_CHAR" ) end select end do ! Take the file out of define mode - call check( nf90_enddef(nciu%id), "netcdf_initialize_output nf90_enddef" ) + call check( nf90_enddef(nc%id), "netcdf_initialize_output nf90_enddef" ) - call check( nf90_put_var(nciu%id, nciu%space_varid, nciu%space_coords, start=[1], count=[NDIM]), "netcdf_initialize_output nf90_put_var space" ) + ! Add in the space dimension coordinates + call check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "netcdf_initialize_output nf90_put_var space" ) end associate return @@ -321,35 +322,35 @@ module subroutine netcdf_open(self, param, readonly) if (readonly) mode = NF90_NOWRITE end if - associate(nciu => self) + associate(nc => self) write(errmsg,*) "netcdf_open nf90_open ",trim(adjustl(param%outfile)) - call check( nf90_open(param%outfile, mode, nciu%id), errmsg) + call check( nf90_open(param%outfile, mode, nc%id), errmsg) ! Dimensions - call check( nf90_inq_dimid(nciu%id, nciu%time_dimname, nciu%time_dimid), "netcdf_open nf90_inq_dimid time_dimid" ) - call check( nf90_inq_dimid(nciu%id, nciu%space_dimname, nciu%space_dimid), "netcdf_open nf90_inq_dimid space_dimid" ) - call check( nf90_inq_dimid(nciu%id, nciu%id_dimname, nciu%id_dimid), "netcdf_open nf90_inq_dimid id_dimid" ) - call check( nf90_inq_dimid(nciu%id, nciu%str_dimname, nciu%str_dimid), "netcdf_open nf90_inq_dimid str_dimid" ) + call check( nf90_inq_dimid(nc%id, nc%time_dimname, nc%time_dimid), "netcdf_open nf90_inq_dimid time_dimid" ) + call check( nf90_inq_dimid(nc%id, nc%space_dimname, nc%space_dimid), "netcdf_open nf90_inq_dimid space_dimid" ) + call check( nf90_inq_dimid(nc%id, nc%id_dimname, nc%id_dimid), "netcdf_open nf90_inq_dimid id_dimid" ) + call check( nf90_inq_dimid(nc%id, nc%str_dimname, nc%str_dimid), "netcdf_open nf90_inq_dimid str_dimid" ) ! Dimension coordinates - call check( nf90_inq_varid(nciu%id, nciu%time_dimname, nciu%time_varid), "netcdf_open nf90_inq_varid time_varid" ) - call check( nf90_inq_varid(nciu%id, nciu%space_dimname, nciu%space_varid), "netcdf_open nf90_inq_varid space_varid" ) - call check( nf90_inq_varid(nciu%id, nciu%id_dimname, nciu%id_varid), "netcdf_open nf90_inq_varid id_varid" ) + call check( nf90_inq_varid(nc%id, nc%time_dimname, nc%time_varid), "netcdf_open nf90_inq_varid time_varid" ) + call check( nf90_inq_varid(nc%id, nc%space_dimname, nc%space_varid), "netcdf_open nf90_inq_varid space_varid" ) + call check( nf90_inq_varid(nc%id, nc%id_dimname, nc%id_varid), "netcdf_open nf90_inq_varid id_varid" ) ! Required Variables - call check( nf90_inq_varid(nciu%id, nciu%name_varname, nciu%name_varid), "netcdf_open nf90_inq_varid name_varid" ) - call check( nf90_inq_varid(nciu%id, nciu%gmass_varname, nciu%Gmass_varid), "netcdf_open nf90_inq_varid Gmass_varid" ) + call check( nf90_inq_varid(nc%id, nc%name_varname, nc%name_varid), "netcdf_open nf90_inq_varid name_varid" ) + call check( nf90_inq_varid(nc%id, nc%gmass_varname, nc%Gmass_varid), "netcdf_open nf90_inq_varid Gmass_varid" ) if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then - call check( nf90_inq_varid(nciu%id, nciu%rh_varname, nciu%rh_varid), "netcdf_open nf90_inq_varid rh_varid" ) - call check( nf90_inq_varid(nciu%id, nciu%vh_varname, nciu%vh_varid), "netcdf_open nf90_inq_varid vh_varid" ) + call check( nf90_inq_varid(nc%id, nc%rh_varname, nc%rh_varid), "netcdf_open nf90_inq_varid rh_varid" ) + call check( nf90_inq_varid(nc%id, nc%vh_varname, nc%vh_varid), "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(nciu%id, nciu%gr_pseudo_vh_varname, nciu%gr_pseudo_vh_varid) - nciu%lpseudo_vel_exists = (status == nf90_noerr) - if (.not.nciu%lpseudo_vel_exists) then + status = nf90_inq_varid(nc%id, nc%gr_pseudo_vh_varname, nc%gr_pseudo_vh_varid) + 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 @@ -357,71 +358,71 @@ module subroutine netcdf_open(self, param, readonly) end if if ((param%out_form == "EL") .or. (param%out_form == "XVEL")) then - call check( nf90_inq_varid(nciu%id, nciu%a_varname, nciu%a_varid), "netcdf_open nf90_inq_varid a_varid" ) - call check( nf90_inq_varid(nciu%id, nciu%e_varname, nciu%e_varid), "netcdf_open nf90_inq_varid e_varid" ) - call check( nf90_inq_varid(nciu%id, nciu%inc_varname, nciu%inc_varid), "netcdf_open nf90_inq_varid inc_varid" ) - call check( nf90_inq_varid(nciu%id, nciu%capom_varname, nciu%capom_varid), "netcdf_open nf90_inq_varid capom_varid" ) - call check( nf90_inq_varid(nciu%id, nciu%omega_varname, nciu%omega_varid), "netcdf_open nf90_inq_varid omega_varid" ) - call check( nf90_inq_varid(nciu%id, nciu%capm_varname, nciu%capm_varid), "netcdf_open nf90_inq_varid capm_varid" ) + call check( nf90_inq_varid(nc%id, nc%a_varname, nc%a_varid), "netcdf_open nf90_inq_varid a_varid" ) + call check( nf90_inq_varid(nc%id, nc%e_varname, nc%e_varid), "netcdf_open nf90_inq_varid e_varid" ) + call check( nf90_inq_varid(nc%id, nc%inc_varname, nc%inc_varid), "netcdf_open nf90_inq_varid inc_varid" ) + call check( nf90_inq_varid(nc%id, nc%capom_varname, nc%capom_varid), "netcdf_open nf90_inq_varid capom_varid" ) + call check( nf90_inq_varid(nc%id, nc%omega_varname, nc%omega_varid), "netcdf_open nf90_inq_varid omega_varid" ) + call check( nf90_inq_varid(nc%id, nc%capm_varname, nc%capm_varid), "netcdf_open nf90_inq_varid capm_varid" ) end if if (param%lclose) then - call check( nf90_inq_varid(nciu%id, nciu%radius_varname, nciu%radius_varid), "netcdf_open nf90_inq_varid radius_varid" ) + call check( nf90_inq_varid(nc%id, nc%radius_varname, nc%radius_varid), "netcdf_open nf90_inq_varid radius_varid" ) end if if (param%lrotation) then - call check( nf90_inq_varid(nciu%id, nciu%Ip_varname, nciu%Ip_varid), "netcdf_open nf90_inq_varid Ip_varid" ) - call check( nf90_inq_varid(nciu%id, nciu%rot_varname, nciu%rot_varid), "netcdf_open nf90_inq_varid rot_varid" ) + call check( nf90_inq_varid(nc%id, nc%Ip_varname, nc%Ip_varid), "netcdf_open nf90_inq_varid Ip_varid" ) + call check( nf90_inq_varid(nc%id, nc%rot_varname, nc%rot_varid), "netcdf_open nf90_inq_varid rot_varid" ) end if ! if (param%ltides) then - ! call check( nf90_inq_varid(nciu%id, nciu%k2_varname, nciu%k2_varid), "netcdf_open nf90_inq_varid k2_varid" ) - ! call check( nf90_inq_varid(nciu%id, nciu%q_varname, nciu%Q_varid), "netcdf_open nf90_inq_varid Q_varid" ) + ! call check( nf90_inq_varid(nc%id, nc%k2_varname, nc%k2_varid), "netcdf_open nf90_inq_varid k2_varid" ) + ! call check( nf90_inq_varid(nc%id, nc%q_varname, nc%Q_varid), "netcdf_open nf90_inq_varid Q_varid" ) ! end if ! Optional Variables if (param%lrhill_present) then - status = nf90_inq_varid(nciu%id, nciu%rhill_varname, nciu%rhill_varid) + 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." end if ! Optional variables The User Doesn't Need to Know About - status = nf90_inq_varid(nciu%id, nciu%npl_varname, nciu%npl_varid) - status = nf90_inq_varid(nciu%id, nciu%ntp_varname, nciu%ntp_varid) - status = nf90_inq_varid(nciu%id, nciu%j2rp2_varname, nciu%j2rp2_varid) - status = nf90_inq_varid(nciu%id, nciu%j4rp4_varname, nciu%j4rp4_varid) - status = nf90_inq_varid(nciu%id, nciu%ptype_varname, nciu%ptype_varid) - status = nf90_inq_varid(nciu%id, nciu%varpi_varname, nciu%varpi_varid) - status = nf90_inq_varid(nciu%id, nciu%lam_varname, nciu%lam_varid) - status = nf90_inq_varid(nciu%id, nciu%f_varname, nciu%f_varid) - status = nf90_inq_varid(nciu%id, nciu%cape_varname, nciu%cape_varid) + status = nf90_inq_varid(nc%id, nc%npl_varname, nc%npl_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) + status = nf90_inq_varid(nc%id, nc%ptype_varname, nc%ptype_varid) + status = nf90_inq_varid(nc%id, nc%varpi_varname, nc%varpi_varid) + status = nf90_inq_varid(nc%id, nc%lam_varname, nc%lam_varid) + status = nf90_inq_varid(nc%id, nc%f_varname, nc%f_varid) + status = nf90_inq_varid(nc%id, nc%cape_varname, nc%cape_varid) if (param%integrator == SYMBA) then - status = nf90_inq_varid(nciu%id, nciu%nplm_varname, nciu%nplm_varid) + status = nf90_inq_varid(nc%id, nc%nplm_varname, nc%nplm_varid) end if if (param%lclose) then - status = nf90_inq_varid(nciu%id, nciu%origin_type_varname, nciu%origin_type_varid) - status = nf90_inq_varid(nciu%id, nciu%origin_time_varname, nciu%origin_time_varid) - status = nf90_inq_varid(nciu%id, nciu%origin_rh_varname, nciu%origin_rh_varid) - status = nf90_inq_varid(nciu%id, nciu%origin_vh_varname, nciu%origin_vh_varid) - status = nf90_inq_varid(nciu%id, nciu%collision_id_varname, nciu%collision_id_varid) - status = nf90_inq_varid(nciu%id, nciu%discard_time_varname, nciu%discard_time_varid) - status = nf90_inq_varid(nciu%id, nciu%discard_rh_varname, nciu%discard_rh_varid) - status = nf90_inq_varid(nciu%id, nciu%discard_vh_varname, nciu%discard_vh_varid) - status = nf90_inq_varid(nciu%id, nciu%discard_body_id_varname, nciu%discard_body_id_varid) + status = nf90_inq_varid(nc%id, nc%origin_type_varname, nc%origin_type_varid) + status = nf90_inq_varid(nc%id, nc%origin_time_varname, nc%origin_time_varid) + status = nf90_inq_varid(nc%id, nc%origin_rh_varname, nc%origin_rh_varid) + status = nf90_inq_varid(nc%id, nc%origin_vh_varname, nc%origin_vh_varid) + status = nf90_inq_varid(nc%id, nc%collision_id_varname, nc%collision_id_varid) + status = nf90_inq_varid(nc%id, nc%discard_time_varname, nc%discard_time_varid) + status = nf90_inq_varid(nc%id, nc%discard_rh_varname, nc%discard_rh_varid) + status = nf90_inq_varid(nc%id, nc%discard_vh_varname, nc%discard_vh_varid) + status = nf90_inq_varid(nc%id, nc%discard_body_id_varname, nc%discard_body_id_varid) end if if (param%lenergy) then - status = nf90_inq_varid(nciu%id, nciu%ke_orb_varname, nciu%KE_orb_varid) - status = nf90_inq_varid(nciu%id, nciu%ke_spin_varname, nciu%KE_spin_varid) - status = nf90_inq_varid(nciu%id, nciu%pe_varname, nciu%PE_varid) - status = nf90_inq_varid(nciu%id, nciu%L_orb_varname, nciu%L_orb_varid) - status = nf90_inq_varid(nciu%id, nciu%L_spin_varname, nciu%L_spin_varid) - status = nf90_inq_varid(nciu%id, nciu%L_escape_varname, nciu%L_escape_varid) - status = nf90_inq_varid(nciu%id, nciu%Ecollisions_varname, nciu%Ecollisions_varid) - status = nf90_inq_varid(nciu%id, nciu%Euntracked_varname, nciu%Euntracked_varid) - status = nf90_inq_varid(nciu%id, nciu%GMescape_varname, nciu%GMescape_varid) + status = nf90_inq_varid(nc%id, nc%ke_orb_varname, nc%KE_orb_varid) + 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%L_orb_varname, nc%L_orb_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%GMescape_varname, nc%GMescape_varid) end if end associate @@ -430,14 +431,14 @@ module subroutine netcdf_open(self, param, readonly) end subroutine netcdf_open - module function netcdf_read_frame_system(self, nciu, param) result(ierr) + module function netcdf_read_frame_system(self, nc, param) result(ierr) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! !! Read a frame (header plus records for each massive body and active test particle) from an output binary file implicit none ! Arguments class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(netcdf_parameters), intent(inout) :: nciu !! Parameters used to identify a particular NetCDF dataset + class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Return integer(I4B) :: ierr !! Error code: returns 0 if the read is successful @@ -448,8 +449,8 @@ module function netcdf_read_frame_system(self, nciu, param) result(ierr) integer(I4B), dimension(:), allocatable :: itemp logical, dimension(:), allocatable :: validmask, tpmask, plmask - call nciu%open(param, readonly=.true.) - call self%read_hdr(nciu, param) + 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) @@ -458,27 +459,27 @@ module function netcdf_read_frame_system(self, nciu, param) result(ierr) tslot = param%ioutput - call check( nf90_inquire_dimension(nciu%id, nciu%id_dimid, len=idmax), "netcdf_read_frame_system nf90_inquire_dimension id_dimid" ) + call check( nf90_inquire_dimension(nc%id, nc%id_dimid, len=idmax), "netcdf_read_frame_system nf90_inquire_dimension id_dimid" ) allocate(rtemp(idmax)) allocate(vectemp(NDIM,idmax)) allocate(itemp(idmax)) allocate(validmask(idmax)) allocate(tpmask(idmax)) allocate(plmask(idmax)) - call check( nf90_inquire_dimension(nciu%id, nciu%time_dimid, len=t_max), "netcdf_read_frame_system nf90_inquire_dimension time_dimid" ) - call check( nf90_inquire_dimension(nciu%id, nciu%str_dimid, len=str_max), "netcdf_read_frame_system nf90_inquire_dimension str_dimid" ) + call check( nf90_inquire_dimension(nc%id, nc%time_dimid, len=t_max), "netcdf_read_frame_system nf90_inquire_dimension time_dimid" ) + call check( nf90_inquire_dimension(nc%id, nc%str_dimid, len=str_max), "netcdf_read_frame_system nf90_inquire_dimension str_dimid" ) ! First filter out only the id slots that contain valid bodies if (param%in_form == "XV") then - call check( nf90_get_var(nciu%id, nciu%rh_varid, vectemp(:,:), start=[1, 1, tslot]), "netcdf_read_frame_system filter pass nf90_getvar rh_varid" ) + call check( nf90_get_var(nc%id, nc%rh_varid, vectemp(:,:), start=[1, 1, tslot]), "netcdf_read_frame_system filter pass nf90_getvar rh_varid" ) validmask(:) = vectemp(1,:) == vectemp(1,:) else - call check( nf90_get_var(nciu%id, nciu%a_varid, rtemp(:), start=[1, tslot]), "netcdf_read_frame_system filter pass nf90_getvar a_varid" ) + call check( nf90_get_var(nc%id, nc%a_varid, rtemp(:), start=[1, tslot]), "netcdf_read_frame_system filter pass nf90_getvar a_varid" ) validmask(:) = rtemp(:) == rtemp(:) end if ! Next, filter only bodies that don't have mass (test particles) - call check( nf90_get_var(nciu%id, nciu%Gmass_varid, rtemp(:), start=[1, tslot]), "netcdf_read_frame_system nf90_getvar tp finder Gmass_varid" ) + call check( nf90_get_var(nc%id, nc%Gmass_varid, rtemp(:), start=[1, tslot]), "netcdf_read_frame_system nf90_getvar tp finder Gmass_varid" ) plmask(:) = rtemp(:) == rtemp(:) .and. validmask(:) tpmask(:) = .not. plmask(:) .and. validmask(:) plmask(1) = .false. ! This is the central body @@ -511,20 +512,20 @@ module function netcdf_read_frame_system(self, nciu, param) result(ierr) ! Now read in each variable and split the outputs by body type if ((param%in_form == "XV") .or. (param%in_form == "XVEL")) then - call check( nf90_get_var(nciu%id, nciu%rh_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "netcdf_read_frame_system nf90_getvar rh_varid" ) + call check( nf90_get_var(nc%id, nc%rh_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "netcdf_read_frame_system nf90_getvar rh_varid" ) do i = 1, NDIM if (npl > 0) pl%rh(i,:) = pack(vectemp(i,:), plmask(:)) if (ntp > 0) tp%rh(i,:) = pack(vectemp(i,:), tpmask(:)) end do - if (param%lgr .and. nciu%lpseudo_vel_exists) then - call check( nf90_get_var(nciu%id, nciu%gr_pseudo_vh_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "netcdf_read_frame_system nf90_getvar gr_pseudo_vh_varid" ) + if (param%lgr .and. nc%lpseudo_vel_exists) then + call check( nf90_get_var(nc%id, nc%gr_pseudo_vh_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "netcdf_read_frame_system nf90_getvar gr_pseudo_vh_varid" ) do i = 1, NDIM if (npl > 0) pl%vh(i,:) = pack(vectemp(i,:), plmask(:)) if (ntp > 0) tp%vh(i,:) = pack(vectemp(i,:), tpmask(:)) end do else - call check( nf90_get_var(nciu%id, nciu%vh_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "netcdf_read_frame_system nf90_getvar vh_varid" ) + call check( nf90_get_var(nc%id, nc%vh_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "netcdf_read_frame_system nf90_getvar vh_varid" ) do i = 1, NDIM if (npl > 0) pl%vh(i,:) = pack(vectemp(i,:), plmask(:)) if (ntp > 0) tp%vh(i,:) = pack(vectemp(i,:), tpmask(:)) @@ -533,40 +534,40 @@ module function netcdf_read_frame_system(self, nciu, param) result(ierr) end if if ((param%in_form == "EL") .or. (param%in_form == "XVEL")) then - call check( nf90_get_var(nciu%id, nciu%a_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar a_varid" ) + call check( nf90_get_var(nc%id, nc%a_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar a_varid" ) if (.not.allocated(pl%a)) allocate(pl%a(npl)) if (.not.allocated(tp%a)) allocate(tp%a(ntp)) if (npl > 0) pl%a(:) = pack(rtemp, plmask) if (ntp > 0) tp%a(:) = pack(rtemp, tpmask) - call check( nf90_get_var(nciu%id, nciu%e_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar e_varid" ) + call check( nf90_get_var(nc%id, nc%e_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar e_varid" ) if (.not.allocated(pl%e)) allocate(pl%e(npl)) if (.not.allocated(tp%e)) allocate(tp%e(ntp)) if (npl > 0) pl%e(:) = pack(rtemp, plmask) if (ntp > 0) tp%e(:) = pack(rtemp, tpmask) - call check( nf90_get_var(nciu%id, nciu%inc_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar inc_varid" ) + call check( nf90_get_var(nc%id, nc%inc_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar inc_varid" ) rtemp = rtemp * DEG2RAD if (.not.allocated(pl%inc)) allocate(pl%inc(npl)) if (.not.allocated(tp%inc)) allocate(tp%inc(ntp)) if (npl > 0) pl%inc(:) = pack(rtemp, plmask) if (ntp > 0) tp%inc(:) = pack(rtemp, tpmask) - call check( nf90_get_var(nciu%id, nciu%capom_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar capom_varid" ) + call check( nf90_get_var(nc%id, nc%capom_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar capom_varid" ) rtemp = rtemp * DEG2RAD if (.not.allocated(pl%capom)) allocate(pl%capom(npl)) if (.not.allocated(tp%capom)) allocate(tp%capom(ntp)) if (npl > 0) pl%capom(:) = pack(rtemp, plmask) if (ntp > 0) tp%capom(:) = pack(rtemp, tpmask) - call check( nf90_get_var(nciu%id, nciu%omega_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar omega_varid" ) + call check( nf90_get_var(nc%id, nc%omega_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar omega_varid" ) rtemp = rtemp * DEG2RAD if (.not.allocated(pl%omega)) allocate(pl%omega(npl)) if (.not.allocated(tp%omega)) allocate(tp%omega(ntp)) if (npl > 0) pl%omega(:) = pack(rtemp, plmask) if (ntp > 0) tp%omega(:) = pack(rtemp, tpmask) - call check( nf90_get_var(nciu%id, nciu%capm_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar capm_varid" ) + call check( nf90_get_var(nc%id, nc%capm_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar capm_varid" ) rtemp = rtemp * DEG2RAD if (.not.allocated(pl%capm)) allocate(pl%capm(npl)) if (.not.allocated(tp%capm)) allocate(tp%capm(ntp)) @@ -575,7 +576,7 @@ module function netcdf_read_frame_system(self, nciu, param) result(ierr) end if - call check( nf90_get_var(nciu%id, nciu%Gmass_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar Gmass_varid" ) + call check( nf90_get_var(nc%id, nc%Gmass_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar Gmass_varid" ) cb%Gmass = rtemp(1) cb%mass = cb%Gmass / param%GU @@ -591,13 +592,13 @@ module function netcdf_read_frame_system(self, nciu, param) result(ierr) pl%mass(:) = pl%Gmass(:) / param%GU if (param%lrhill_present) then - call check( nf90_get_var(nciu%id, nciu%rhill_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar rhill_varid" ) + call check( nf90_get_var(nc%id, nc%rhill_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar rhill_varid" ) pl%rhill(:) = pack(rtemp, plmask) end if end if if (param%lclose) then - call check( nf90_get_var(nciu%id, nciu%radius_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar radius_varid" ) + call check( nf90_get_var(nc%id, nc%radius_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar radius_varid" ) cb%radius = rtemp(1) ! Set initial central body radius for SyMBA bookkeeping @@ -612,13 +613,13 @@ module function netcdf_read_frame_system(self, nciu, param) result(ierr) end if if (param%lrotation) then - call check( nf90_get_var(nciu%id, nciu%Ip_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "netcdf_read_frame_system nf90_getvar Ip_varid" ) + call check( nf90_get_var(nc%id, nc%Ip_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "netcdf_read_frame_system nf90_getvar Ip_varid" ) cb%Ip(:) = vectemp(:,1) do i = 1, NDIM if (npl > 0) pl%Ip(i,:) = pack(vectemp(i,:), plmask(:)) end do - call check( nf90_get_var(nciu%id, nciu%rot_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "netcdf_read_frame_system nf90_getvar rot_varid" ) + call check( nf90_get_var(nc%id, nc%rot_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "netcdf_read_frame_system nf90_getvar rot_varid" ) cb%rot(:) = vectemp(:,1) do i = 1, NDIM if (npl > 0) pl%rot(i,:) = pack(vectemp(i,:), plmask(:)) @@ -632,37 +633,37 @@ module function netcdf_read_frame_system(self, nciu, param) result(ierr) end if ! if (param%ltides) then - ! call check( nf90_get_var(nciu%id, nciu%k2_varid, rtemp, start=[1, tslot]), "netcdf_read_frame_system nf90_getvar k2_varid" ) + ! call check( nf90_get_var(nc%id, nc%k2_varid, rtemp, start=[1, tslot]), "netcdf_read_frame_system nf90_getvar k2_varid" ) ! cb%k2 = rtemp(1) ! if (npl > 0) pl%k2(:) = pack(rtemp, plmask) - ! call check( nf90_get_var(nciu%id, nciu%Q_varid, rtemp, start=[1, tslot]), "netcdf_read_frame_system nf90_getvar Q_varid" ) + ! call check( nf90_get_var(nc%id, nc%Q_varid, rtemp, start=[1, tslot]), "netcdf_read_frame_system nf90_getvar Q_varid" ) ! cb%Q = rtemp(1) ! if (npl > 0) pl%Q(:) = pack(rtemp, plmask) ! end if - status = nf90_inq_varid(nciu%id, nciu%j2rp2_varname, nciu%j2rp2_varid) + status = nf90_inq_varid(nc%id, nc%j2rp2_varname, nc%j2rp2_varid) if (status == nf90_noerr) then - call check( nf90_get_var(nciu%id, nciu%j2rp2_varid, cb%j2rp2, start=[tslot]), "netcdf_read_frame_system nf90_getvar j2rp2_varid" ) + call check( nf90_get_var(nc%id, nc%j2rp2_varid, cb%j2rp2, start=[tslot]), "netcdf_read_frame_system nf90_getvar j2rp2_varid" ) else cb%j2rp2 = 0.0_DP end if - status = nf90_inq_varid(nciu%id, nciu%j4rp4_varname, nciu%j4rp4_varid) + status = nf90_inq_varid(nc%id, nc%j4rp4_varname, nc%j4rp4_varid) if (status == nf90_noerr) then - call check( nf90_get_var(nciu%id, nciu%j4rp4_varid, cb%j4rp4, start=[tslot]), "netcdf_read_frame_system nf90_getvar j4rp4_varid" ) + call check( nf90_get_var(nc%id, nc%j4rp4_varid, cb%j4rp4, start=[tslot]), "netcdf_read_frame_system nf90_getvar j4rp4_varid" ) else cb%j4rp4 = 0.0_DP end if - call self%read_particle_info(nciu, param, plmask, tpmask) + call self%read_particle_info(nc, param, plmask, tpmask) if (param%in_form == "EL") then call pl%el2xv(cb) call tp%el2xv(cb) end if ! if this is a GR-enabled run, check to see if we got the pseudovelocities in. Otherwise, we'll need to generate them. - if (param%lgr .and. .not.(nciu%lpseudo_vel_exists)) then + if (param%lgr .and. .not.(nc%lpseudo_vel_exists)) then call pl%set_mu(cb) call tp%set_mu(cb) call pl%v2pv(param) @@ -671,7 +672,7 @@ module function netcdf_read_frame_system(self, nciu, param) result(ierr) end associate - call nciu%close() + call nc%close() ierr = 0 return @@ -682,7 +683,7 @@ module function netcdf_read_frame_system(self, nciu, param) result(ierr) end function netcdf_read_frame_system - module subroutine netcdf_read_hdr_system(self, nciu, param) + module subroutine netcdf_read_hdr_system(self, nc, param) !! author: David A. Minton !! !! Reads header information (variables that change with time, but not particle id). @@ -691,7 +692,7 @@ module subroutine netcdf_read_hdr_system(self, nciu, param) implicit none ! Arguments class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object - class(netcdf_parameters), intent(inout) :: nciu !! Parameters used to for writing a NetCDF dataset to file + class(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, status, idmax @@ -700,15 +701,15 @@ module subroutine netcdf_read_hdr_system(self, nciu, param) tslot = param%ioutput - call check( nf90_inquire_dimension(nciu%id, nciu%id_dimid, len=idmax), "netcdf_read_hdr_system nf90_inquire_dimension id_dimid" ) - call check( nf90_get_var(nciu%id, nciu%time_varid, self%t, start=[tslot]), "netcdf_read_hdr_system nf90_getvar time_varid" ) + call check( nf90_inquire_dimension(nc%id, nc%id_dimid, len=idmax), "netcdf_read_hdr_system nf90_inquire_dimension id_dimid" ) + call check( nf90_get_var(nc%id, nc%time_varid, self%t, start=[tslot]), "netcdf_read_hdr_system nf90_getvar time_varid" ) allocate(gmtemp(idmax)) allocate(tpmask(idmax)) allocate(plmask(idmax)) allocate(plmmask(idmax)) - call check( nf90_get_var(nciu%id, nciu%Gmass_varid, gmtemp, start=[1,1], count=[idmax,1]), "netcdf_read_hdr_system nf90_getvar Gmass_varid" ) + call check( nf90_get_var(nc%id, nc%Gmass_varid, gmtemp, start=[1,1], count=[idmax,1]), "netcdf_read_hdr_system nf90_getvar Gmass_varid" ) plmask(:) = gmtemp(:) == gmtemp(:) tpmask(:) = .not. plmask(:) @@ -721,26 +722,26 @@ module subroutine netcdf_read_hdr_system(self, nciu, param) endwhere end select - status = nf90_inq_varid(nciu%id, nciu%npl_varname, nciu%npl_varid) + status = nf90_inq_varid(nc%id, nc%npl_varname, nc%npl_varid) if (status == nf90_noerr) then - call check( nf90_get_var(nciu%id, nciu%npl_varid, self%pl%nbody, start=[tslot]), "netcdf_read_hdr_system nf90_getvar npl_varid" ) + call check( nf90_get_var(nc%id, nc%npl_varid, self%pl%nbody, start=[tslot]), "netcdf_read_hdr_system nf90_getvar npl_varid" ) else self%pl%nbody = count(plmask(:)) end if - status = nf90_inq_varid(nciu%id, nciu%ntp_varname, nciu%ntp_varid) + status = nf90_inq_varid(nc%id, nc%ntp_varname, nc%ntp_varid) if (status == nf90_noerr) then - call check( nf90_get_var(nciu%id, nciu%ntp_varid, self%tp%nbody, start=[tslot]), "netcdf_read_hdr_system nf90_getvar ntp_varid" ) + call check( nf90_get_var(nc%id, nc%ntp_varid, self%tp%nbody, start=[tslot]), "netcdf_read_hdr_system nf90_getvar ntp_varid" ) else self%tp%nbody = count(tpmask(:)) end if if (param%integrator == SYMBA) then - status = nf90_inq_varid(nciu%id, nciu%nplm_varname, nciu%nplm_varid) + status = nf90_inq_varid(nc%id, nc%nplm_varname, nc%nplm_varid) select type(pl => self%pl) class is (symba_pl) if (status == nf90_noerr) then - call check( nf90_get_var(nciu%id, nciu%nplm_varid, pl%nplm, start=[tslot]), "netcdf_read_hdr_system nf90_getvar nplm_varid" ) + call check( nf90_get_var(nc%id, nc%nplm_varid, pl%nplm, start=[tslot]), "netcdf_read_hdr_system nf90_getvar nplm_varid" ) else pl%nplm = count(plmmask(:)) end if @@ -748,38 +749,38 @@ module subroutine netcdf_read_hdr_system(self, nciu, param) end if if (param%lenergy) then - status = nf90_inq_varid(nciu%id, nciu%ke_orb_varname, nciu%KE_orb_varid) - if (status == nf90_noerr) call check( nf90_get_var(nciu%id, nciu%KE_orb_varid, self%ke_orbit, start=[tslot]), "netcdf_read_hdr_system nf90_getvar KE_orb_varid" ) - status = nf90_inq_varid(nciu%id, nciu%ke_spin_varname, nciu%KE_spin_varid) - if (status == nf90_noerr) call check( nf90_get_var(nciu%id, nciu%KE_spin_varid, self%ke_spin, start=[tslot]), "netcdf_read_hdr_system nf90_getvar KE_spin_varid" ) - status = nf90_inq_varid(nciu%id, nciu%pe_varname, nciu%PE_varid) - if (status == nf90_noerr) call check( nf90_get_var(nciu%id, nciu%PE_varid, self%pe, start=[tslot]), "netcdf_read_hdr_system nf90_getvar PE_varid" ) - status = nf90_inq_varid(nciu%id, nciu%L_orb_varname, nciu%L_orb_varid) - if (status == nf90_noerr) call check( nf90_get_var(nciu%id, nciu%L_orb_varid, self%Lorbit(:), start=[1,tslot], count=[NDIM,1]), "netcdf_read_hdr_system nf90_getvar L_orb_varid" ) - status = nf90_inq_varid(nciu%id, nciu%L_spin_varname, nciu%L_spin_varid) - if (status == nf90_noerr) call check( nf90_get_var(nciu%id, nciu%L_spin_varid, self%Lspin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_read_hdr_system nf90_getvar L_spin_varid" ) - status = nf90_inq_varid(nciu%id, nciu%L_escape_varname, nciu%L_escape_varid) - if (status == nf90_noerr) call check( nf90_get_var(nciu%id, nciu%L_escape_varid, self%Lescape(:), start=[1, tslot], count=[NDIM,1]), "netcdf_read_hdr_system nf90_getvar L_escape_varid" ) - status = nf90_inq_varid(nciu%id, nciu%Ecollisions_varname, nciu%Ecollisions_varid) - if (status == nf90_noerr) call check( nf90_get_var(nciu%id, nciu%Ecollisions_varid, self%Ecollisions, start=[tslot]), "netcdf_read_hdr_system nf90_getvar Ecollisions_varid" ) - status = nf90_inq_varid(nciu%id, nciu%Euntracked_varname, nciu%Euntracked_varid) - if (status == nf90_noerr) call check( nf90_get_var(nciu%id, nciu%Euntracked_varid, self%Euntracked, start=[tslot]), "netcdf_read_hdr_system nf90_getvar Euntracked_varid" ) - status = nf90_inq_varid(nciu%id, nciu%GMescape_varname, nciu%GMescape_varid) - if (status == nf90_noerr) call check( nf90_get_var(nciu%id, nciu%GMescape_varid, self%GMescape, start=[tslot]), "netcdf_read_hdr_system nf90_getvar GMescape_varid" ) + status = nf90_inq_varid(nc%id, nc%ke_orb_varname, nc%KE_orb_varid) + if (status == nf90_noerr) call check( nf90_get_var(nc%id, nc%KE_orb_varid, self%ke_orbit, start=[tslot]), "netcdf_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 check( nf90_get_var(nc%id, nc%KE_spin_varid, self%ke_spin, start=[tslot]), "netcdf_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 check( nf90_get_var(nc%id, nc%PE_varid, self%pe, start=[tslot]), "netcdf_read_hdr_system nf90_getvar PE_varid" ) + status = nf90_inq_varid(nc%id, nc%L_orb_varname, nc%L_orb_varid) + if (status == nf90_noerr) call check( nf90_get_var(nc%id, nc%L_orb_varid, self%Lorbit(:), start=[1,tslot], count=[NDIM,1]), "netcdf_read_hdr_system nf90_getvar L_orb_varid" ) + status = nf90_inq_varid(nc%id, nc%L_spin_varname, nc%L_spin_varid) + if (status == nf90_noerr) call check( nf90_get_var(nc%id, nc%L_spin_varid, self%Lspin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_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 check( nf90_get_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1, tslot], count=[NDIM,1]), "netcdf_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 check( nf90_get_var(nc%id, nc%Ecollisions_varid, self%Ecollisions, start=[tslot]), "netcdf_read_hdr_system nf90_getvar Ecollisions_varid" ) + status = nf90_inq_varid(nc%id, nc%Euntracked_varname, nc%Euntracked_varid) + if (status == nf90_noerr) call check( nf90_get_var(nc%id, nc%Euntracked_varid, self%Euntracked, start=[tslot]), "netcdf_read_hdr_system nf90_getvar Euntracked_varid" ) + status = nf90_inq_varid(nc%id, nc%GMescape_varname, nc%GMescape_varid) + if (status == nf90_noerr) call check( nf90_get_var(nc%id, nc%GMescape_varid, self%GMescape, start=[tslot]), "netcdf_read_hdr_system nf90_getvar GMescape_varid" ) end if return end subroutine netcdf_read_hdr_system - module subroutine netcdf_read_particle_info_system(self, nciu, param, plmask, tpmask) + module subroutine netcdf_read_particle_info_system(self, nc, param, plmask, tpmask) !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton !! !! Reads particle information metadata from file implicit none ! Arguments class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object - class(netcdf_parameters), intent(inout) :: nciu !! Parameters used to identify a particular NetCDF dataset + class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters logical, dimension(:), intent(in) :: plmask !! Logical array indicating which index values belong to massive bodies logical, dimension(:), intent(in) :: tpmask !! Logical array indicating which index values belong to test particles @@ -819,12 +820,12 @@ module subroutine netcdf_read_particle_info_system(self, nciu, param, plmask, tp tpind(:) = pack([(i, i = 1, idmax)], tpmask(:)) end if - call check( nf90_get_var(nciu%id, nciu%id_varid, itemp), "netcdf_read_particle_info_system nf90_getvar id_varid" ) + call check( nf90_get_var(nc%id, nc%id_varid, itemp), "netcdf_read_particle_info_system nf90_getvar id_varid" ) cb%id = itemp(1) pl%id(:) = pack(itemp, plmask) tp%id(:) = pack(itemp, tpmask) - call check( nf90_get_var(nciu%id, nciu%name_varid, ctemp, count=[NAMELEN, idmax]), "netcdf_read_particle_info_system nf90_getvar name_varid" ) + call check( nf90_get_var(nc%id, nc%name_varid, ctemp, count=[NAMELEN, idmax]), "netcdf_read_particle_info_system nf90_getvar name_varid" ) call cb%info%set_value(name=ctemp(1)) do i = 1, npl call pl%info(i)%set_value(name=ctemp(plind(i))) @@ -833,7 +834,7 @@ module subroutine netcdf_read_particle_info_system(self, nciu, param, plmask, tp call tp%info(i)%set_value(name=ctemp(tpind(i))) end do - status = nf90_get_var(nciu%id, nciu%ptype_varid, ctemp, count=[NAMELEN, idmax]) + status = nf90_get_var(nc%id, nc%ptype_varid, ctemp, count=[NAMELEN, idmax]) if (status /= nf90_noerr) then ! Set default particle types call cb%info%set_value(particle_type=CB_TYPE_NAME) @@ -872,9 +873,9 @@ module subroutine netcdf_read_particle_info_system(self, nciu, param, plmask, tp if (param%lclose) then - status = nf90_inq_varid(nciu%id, nciu%origin_type_varname, nciu%origin_type_varid) + status = nf90_inq_varid(nc%id, nc%origin_type_varname, nc%origin_type_varid) if (status == nf90_noerr) then - call check( nf90_get_var(nciu%id, nciu%origin_type_varid, ctemp, count=[NAMELEN, idmax]), "netcdf_read_particle_info_system nf90_getvar origin_type_varid" ) + call check( nf90_get_var(nc%id, nc%origin_type_varid, ctemp, count=[NAMELEN, idmax]), "netcdf_read_particle_info_system nf90_getvar origin_type_varid" ) else ctemp = "Initial Conditions" end if @@ -887,9 +888,9 @@ module subroutine netcdf_read_particle_info_system(self, nciu, param, plmask, tp call tp%info(i)%set_value(origin_type=ctemp(tpind(i))) end do - status = nf90_inq_varid(nciu%id, nciu%origin_time_varname, nciu%origin_time_varid) + status = nf90_inq_varid(nc%id, nc%origin_time_varname, nc%origin_time_varid) if (status == nf90_noerr) then - call check( nf90_get_var(nciu%id, nciu%origin_time_varid, rtemp), "netcdf_read_particle_info_system nf90_getvar origin_time_varid" ) + call check( nf90_get_var(nc%id, nc%origin_time_varid, rtemp), "netcdf_read_particle_info_system nf90_getvar origin_time_varid" ) else rtemp = param%t0 end if @@ -902,11 +903,11 @@ module subroutine netcdf_read_particle_info_system(self, nciu, param, plmask, tp call tp%info(i)%set_value(origin_time=rtemp(tpind(i))) end do - status = nf90_inq_varid(nciu%id, nciu%origin_rh_varname, nciu%origin_rh_varid) + status = nf90_inq_varid(nc%id, nc%origin_rh_varname, nc%origin_rh_varid) if (status == nf90_noerr) then - call check( nf90_get_var(nciu%id, nciu%origin_rh_varid, vectemp(:,:)), "netcdf_read_particle_info_system nf90_getvar origin_rh_varid" ) + call check( nf90_get_var(nc%id, nc%origin_rh_varid, vectemp(:,:)), "netcdf_read_particle_info_system nf90_getvar origin_rh_varid" ) else if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then - call check( nf90_get_var(nciu%id, nciu%rh_varid, vectemp(:,:)), "netcdf_read_particle_info_system nf90_getvar rh_varid" ) + call check( nf90_get_var(nc%id, nc%rh_varid, vectemp(:,:)), "netcdf_read_particle_info_system nf90_getvar rh_varid" ) else vectemp(:,:) = 0._DP end if @@ -918,11 +919,11 @@ module subroutine netcdf_read_particle_info_system(self, nciu, param, plmask, tp call tp%info(i)%set_value(origin_rh=vectemp(:,tpind(i))) end do - status = nf90_inq_varid(nciu%id, nciu%origin_vh_varname, nciu%origin_vh_varid) + status = nf90_inq_varid(nc%id, nc%origin_vh_varname, nc%origin_vh_varid) if (status == nf90_noerr) then - call check( nf90_get_var(nciu%id, nciu%origin_vh_varid, vectemp(:,:)), "netcdf_read_particle_info_system nf90_getvar origin_vh_varid" ) + call check( nf90_get_var(nc%id, nc%origin_vh_varid, vectemp(:,:)), "netcdf_read_particle_info_system nf90_getvar origin_vh_varid" ) else if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then - call check( nf90_get_var(nciu%id, nciu%vh_varid, vectemp(:,:)), "netcdf_read_particle_info_system nf90_getvar vh_varid" ) + call check( nf90_get_var(nc%id, nc%vh_varid, vectemp(:,:)), "netcdf_read_particle_info_system nf90_getvar vh_varid" ) else vectemp(:,:) = 0._DP end if @@ -934,9 +935,9 @@ module subroutine netcdf_read_particle_info_system(self, nciu, param, plmask, tp call tp%info(i)%set_value(origin_vh=vectemp(:,tpind(i))) end do - status = nf90_inq_varid(nciu%id, nciu%collision_id_varname, nciu%collision_id_varid) + status = nf90_inq_varid(nc%id, nc%collision_id_varname, nc%collision_id_varid) if (status == nf90_noerr) then - call check( nf90_get_var(nciu%id, nciu%collision_id_varid, itemp), "netcdf_read_particle_info_system nf90_getvar collision_id_varid" ) + call check( nf90_get_var(nc%id, nc%collision_id_varid, itemp), "netcdf_read_particle_info_system nf90_getvar collision_id_varid" ) else itemp = 0.0_DP end if @@ -948,9 +949,9 @@ module subroutine netcdf_read_particle_info_system(self, nciu, param, plmask, tp call tp%info(i)%set_value(collision_id=itemp(tpind(i))) end do - status = nf90_inq_varid(nciu%id, nciu%discard_time_varname, nciu%discard_time_varid) + status = nf90_inq_varid(nc%id, nc%discard_time_varname, nc%discard_time_varid) if (status == nf90_noerr) then - call check( nf90_get_var(nciu%id, nciu%discard_time_varid, rtemp), "netcdf_read_particle_info_system nf90_getvar discard_time_varid" ) + call check( nf90_get_var(nc%id, nc%discard_time_varid, rtemp), "netcdf_read_particle_info_system nf90_getvar discard_time_varid" ) else rtemp = 0.0_DP end if @@ -963,9 +964,9 @@ module subroutine netcdf_read_particle_info_system(self, nciu, param, plmask, tp call tp%info(i)%set_value(discard_time=rtemp(tpind(i))) end do - status = nf90_inq_varid(nciu%id, nciu%discard_rh_varname, nciu%discard_rh_varid) + status = nf90_inq_varid(nc%id, nc%discard_rh_varname, nc%discard_rh_varid) if (status == nf90_noerr) then - call check( nf90_get_var(nciu%id, nciu%discard_rh_varid, vectemp(:,:)), "netcdf_read_particle_info_system nf90_getvar discard_rh_varid" ) + call check( nf90_get_var(nc%id, nc%discard_rh_varid, vectemp(:,:)), "netcdf_read_particle_info_system nf90_getvar discard_rh_varid" ) else vectemp(:,:) = 0.0_DP end if @@ -977,9 +978,9 @@ module subroutine netcdf_read_particle_info_system(self, nciu, param, plmask, tp call tp%info(i)%set_value(discard_rh=vectemp(:,tpind(i))) end do - status = nf90_inq_varid(nciu%id, nciu%discard_vh_varname, nciu%discard_vh_varid) + status = nf90_inq_varid(nc%id, nc%discard_vh_varname, nc%discard_vh_varid) if (status == nf90_noerr) then - call check( nf90_get_var(nciu%id, nciu%discard_vh_varid, vectemp(:,:)), "netcdf_read_particle_info_system nf90_getvar discard_vh_varid" ) + call check( nf90_get_var(nc%id, nc%discard_vh_varid, vectemp(:,:)), "netcdf_read_particle_info_system nf90_getvar discard_vh_varid" ) else vectemp(:,:) = 0.0_DP end if @@ -1013,7 +1014,7 @@ module subroutine netcdf_sync(self) end subroutine netcdf_sync - module subroutine netcdf_write_frame_base(self, nciu, param) + module subroutine netcdf_write_frame_base(self, nc, param) !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton !! !! Write a frame of output of either test particle or massive body data to the binary output file @@ -1021,7 +1022,7 @@ module subroutine netcdf_write_frame_base(self, nciu, param) implicit none ! Arguments class(swiftest_base), intent(in) :: self !! Swiftest particle object - class(netcdf_parameters), intent(inout) :: nciu !! Parameters used to identify a particular NetCDF dataset + class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i, j, tslot, idslot, old_mode @@ -1029,11 +1030,11 @@ module subroutine netcdf_write_frame_base(self, nciu, param) real(DP), dimension(NDIM) :: vh !! Temporary variable to store heliocentric velocity values when converting from pseudovelocity in GR-enabled runs real(DP) :: a, e, inc, omega, capom, capm, varpi, lam, f, cape, capf - call self%write_info(nciu, param) + call self%write_info(nc, param) tslot = param%ioutput - call check( nf90_set_fill(nciu%id, nf90_nofill, old_mode), "netcdf_write_frame_base nf90_set_fill" ) + call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "netcdf_write_frame_base nf90_set_fill" ) select type(self) class is (swiftest_body) associate(n => self%nbody) @@ -1049,13 +1050,13 @@ module subroutine netcdf_write_frame_base(self, nciu, param) if (param%lgr) call gr_pseudovel2vel(param, self%mu(j), self%rh(:, j), self%vh(:, j), vh(:)) if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then - call check( nf90_put_var(nciu%id, nciu%rh_varid, self%rh(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var rh_varid" ) + call check( nf90_put_var(nc%id, nc%rh_varid, self%rh(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var rh_varid" ) if (param%lgr) then !! Convert from pseudovelocity to heliocentric without replacing the current value of pseudovelocity - call check( nf90_put_var(nciu%id, nciu%vh_varid, vh(:), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var vh_varid" ) - call check( nf90_put_var(nciu%id, nciu%gr_pseudo_vh_varid, self%vh(:, j), start=[1,idslot, tslot],count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var gr_pseudo_vhx_varid" ) + call check( nf90_put_var(nc%id, nc%vh_varid, vh(:), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var vh_varid" ) + call check( nf90_put_var(nc%id, nc%gr_pseudo_vh_varid, self%vh(:, j), start=[1,idslot, tslot],count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var gr_pseudo_vhx_varid" ) else - call check( nf90_put_var(nciu%id, nciu%vh_varid, self%vh(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var vh_varid" ) + call check( nf90_put_var(nc%id, nc%vh_varid, self%vh(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var vh_varid" ) end if end if @@ -1069,36 +1070,36 @@ module subroutine netcdf_write_frame_base(self, nciu, param) self%vh(1,j), self%vh(2,j), self%vh(3,j), & a, e, inc, capom, omega, capm, varpi, lam, f, cape, capf) end if - call check( nf90_put_var(nciu%id, nciu%a_varid, a, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body a_varid" ) - call check( nf90_put_var(nciu%id, nciu%e_varid, e, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body e_varid" ) - call check( nf90_put_var(nciu%id, nciu%inc_varid, inc * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body inc_varid" ) - call check( nf90_put_var(nciu%id, nciu%capom_varid, capom * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body capom_varid" ) - call check( nf90_put_var(nciu%id, nciu%omega_varid, omega * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body omega_varid" ) - call check( nf90_put_var(nciu%id, nciu%capm_varid, capm * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body capm_varid" ) - call check( nf90_put_var(nciu%id, nciu%varpi_varid, varpi * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body varpi_varid" ) - call check( nf90_put_var(nciu%id, nciu%lam_varid, lam * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body lam_varid" ) - call check( nf90_put_var(nciu%id, nciu%f_varid, f * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body f_varid" ) + call check( nf90_put_var(nc%id, nc%a_varid, a, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body a_varid" ) + call check( nf90_put_var(nc%id, nc%e_varid, e, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body e_varid" ) + call check( nf90_put_var(nc%id, nc%inc_varid, inc * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body inc_varid" ) + call check( nf90_put_var(nc%id, nc%capom_varid, capom * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body capom_varid" ) + call check( nf90_put_var(nc%id, nc%omega_varid, omega * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body omega_varid" ) + call check( nf90_put_var(nc%id, nc%capm_varid, capm * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body capm_varid" ) + call check( nf90_put_var(nc%id, nc%varpi_varid, varpi * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body varpi_varid" ) + call check( nf90_put_var(nc%id, nc%lam_varid, lam * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body lam_varid" ) + call check( nf90_put_var(nc%id, nc%f_varid, f * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body f_varid" ) if (e < 1.0_DP) then - call check( nf90_put_var(nciu%id, nciu%cape_varid, cape * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body cape_varid" ) + call check( nf90_put_var(nc%id, nc%cape_varid, cape * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body cape_varid" ) else if (e > 1.0_DP) then - call check( nf90_put_var(nciu%id, nciu%cape_varid, capf * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body (capf) cape_varid" ) + call check( nf90_put_var(nc%id, nc%cape_varid, capf * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body (capf) cape_varid" ) end if end if select type(self) class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body - call check( nf90_put_var(nciu%id, nciu%Gmass_varid, self%Gmass(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body Gmass_varid" ) + call check( nf90_put_var(nc%id, nc%Gmass_varid, self%Gmass(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body Gmass_varid" ) if (param%lrhill_present) then - call check( nf90_put_var(nciu%id, nciu%rhill_varid, self%rhill(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body rhill_varid" ) + call check( nf90_put_var(nc%id, nc%rhill_varid, self%rhill(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body rhill_varid" ) end if - if (param%lclose) call check( nf90_put_var(nciu%id, nciu%radius_varid, self%radius(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body radius_varid" ) + if (param%lclose) call check( nf90_put_var(nc%id, nc%radius_varid, self%radius(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body radius_varid" ) if (param%lrotation) then - call check( nf90_put_var(nciu%id, nciu%Ip_varid, self%Ip(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var body Ip_varid" ) - call check( nf90_put_var(nciu%id, nciu%rot_varid, self%rot(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var body rotx_varid" ) + call check( nf90_put_var(nc%id, nc%Ip_varid, self%Ip(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var body Ip_varid" ) + call check( nf90_put_var(nc%id, nc%rot_varid, self%rot(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var body rotx_varid" ) end if ! if (param%ltides) then - ! call check( nf90_put_var(nciu%id, nciu%k2_varid, self%k2(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body k2_varid" ) - ! call check( nf90_put_var(nciu%id, nciu%Q_varid, self%Q(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body Q_varid" ) + ! call check( nf90_put_var(nc%id, nc%k2_varid, self%k2(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body k2_varid" ) + ! call check( nf90_put_var(nc%id, nc%Q_varid, self%Q(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body Q_varid" ) ! end if end select @@ -1106,55 +1107,55 @@ module subroutine netcdf_write_frame_base(self, nciu, param) end associate class is (swiftest_cb) idslot = self%id + 1 - call check( nf90_put_var(nciu%id, nciu%id_varid, self%id, start=[idslot]), "netcdf_write_frame_base nf90_put_var cb id_varid" ) + call check( nf90_put_var(nc%id, nc%id_varid, self%id, start=[idslot]), "netcdf_write_frame_base nf90_put_var cb id_varid" ) - call check( nf90_put_var(nciu%id, nciu%Gmass_varid, self%Gmass, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var cb Gmass_varid" ) - if (param%lclose) call check( nf90_put_var(nciu%id, nciu%radius_varid, self%radius, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var cb radius_varid" ) - call check( nf90_put_var(nciu%id, nciu%j2rp2_varid, self%j2rp2, start=[tslot]), "netcdf_write_frame_base nf90_put_var cb j2rp2_varid" ) - call check( nf90_put_var(nciu%id, nciu%j4rp4_varid, self%j4rp4, start=[tslot]), "netcdf_write_frame_base nf90_put_var cb j4rp4_varid" ) + call check( nf90_put_var(nc%id, nc%Gmass_varid, self%Gmass, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var cb Gmass_varid" ) + if (param%lclose) call check( nf90_put_var(nc%id, nc%radius_varid, self%radius, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var cb radius_varid" ) + call check( nf90_put_var(nc%id, nc%j2rp2_varid, self%j2rp2, start=[tslot]), "netcdf_write_frame_base nf90_put_var cb j2rp2_varid" ) + call check( nf90_put_var(nc%id, nc%j4rp4_varid, self%j4rp4, start=[tslot]), "netcdf_write_frame_base nf90_put_var cb j4rp4_varid" ) if (param%lrotation) then - call check( nf90_put_var(nciu%id, nciu%Ip_varid, self%Ip(:), start=[1, idslot, tslot], count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var cb Ip_varid" ) - call check( nf90_put_var(nciu%id, nciu%rot_varid, self%rot(:), start=[1, idslot, tslot], count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var cb rot_varid" ) + call check( nf90_put_var(nc%id, nc%Ip_varid, self%Ip(:), start=[1, idslot, tslot], count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var cb Ip_varid" ) + call check( nf90_put_var(nc%id, nc%rot_varid, self%rot(:), start=[1, idslot, tslot], count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var cb rot_varid" ) end if ! if (param%ltides) then - ! call check( nf90_put_var(nciu%id, nciu%k2_varid, self%k2, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var cb k2_varid" ) - ! call check( nf90_put_var(nciu%id, nciu%Q_varid, self%Q, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var cb Q_varid" ) + ! call check( nf90_put_var(nc%id, nc%k2_varid, self%k2, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var cb k2_varid" ) + ! call check( nf90_put_var(nc%id, nc%Q_varid, self%Q, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var cb Q_varid" ) ! end if end select - call check( nf90_set_fill(nciu%id, old_mode, old_mode), "netcdf_write_frame_base nf90_set_fill old_mode" ) + call check( nf90_set_fill(nc%id, old_mode, old_mode), "netcdf_write_frame_base nf90_set_fill old_mode" ) return end subroutine netcdf_write_frame_base - module subroutine netcdf_write_frame_system(self, nciu, param) + module subroutine netcdf_write_frame_system(self, nc, param) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! !! Write a frame (header plus records for each massive body and active test particle) to a output binary file implicit none ! Arguments class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(netcdf_parameters), intent(inout) :: nciu !! Parameters used to identify a particular NetCDF dataset + class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - call self%write_hdr(nciu, param) - call self%cb%write_frame(nciu, param) - call self%pl%write_frame(nciu, param) - call self%tp%write_frame(nciu, param) + call self%write_hdr(nc, param) + call self%cb%write_frame(nc, param) + call self%pl%write_frame(nc, param) + call self%tp%write_frame(nc, param) return end subroutine netcdf_write_frame_system - module subroutine netcdf_write_info_base(self, nciu, param) + module subroutine netcdf_write_info_base(self, nc, param) !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton !! !! Write all current particle to file implicit none ! Arguments class(swiftest_base), intent(in) :: self !! Swiftest particle object - class(netcdf_parameters), intent(inout) :: nciu !! Parameters used to identify a particular NetCDF dataset + class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i, j, idslot, old_mode @@ -1162,7 +1163,7 @@ module subroutine netcdf_write_info_base(self, nciu, param) character(len=NAMELEN) :: charstring ! This string of spaces of length NAMELEN is used to clear out any old data left behind inside the string variables - call check( nf90_set_fill(nciu%id, nf90_nofill, old_mode), "netcdf_write_info_base nf90_set_fill nf90_nofill" ) + call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "netcdf_write_info_base nf90_set_fill nf90_nofill" ) select type(self) class is (swiftest_body) @@ -1173,25 +1174,25 @@ module subroutine netcdf_write_info_base(self, nciu, param) do i = 1, n j = ind(i) idslot = self%id(j) + 1 - call check( nf90_put_var(nciu%id, nciu%id_varid, self%id(j), start=[idslot]), "netcdf_write_info_base nf90_put_var id_varid" ) + call check( nf90_put_var(nc%id, nc%id_varid, self%id(j), start=[idslot]), "netcdf_write_info_base nf90_put_var id_varid" ) charstring = trim(adjustl(self%info(j)%name)) - call check( nf90_put_var(nciu%id, nciu%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "netcdf_write_info_base nf90_put_var name_varid" ) + call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "netcdf_write_info_base nf90_put_var name_varid" ) charstring = trim(adjustl(self%info(j)%particle_type)) - call check( nf90_put_var(nciu%id, nciu%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "netcdf_write_info_base nf90_put_var particle_type_varid" ) + call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "netcdf_write_info_base nf90_put_var particle_type_varid" ) if (param%lclose) then charstring = trim(adjustl(self%info(j)%origin_type)) - call check( nf90_put_var(nciu%id, nciu%origin_type_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "netcdf_write_info_base nf90_put_var origin_type_varid" ) - call check( nf90_put_var(nciu%id, nciu%origin_time_varid, self%info(j)%origin_time, start=[idslot]), "netcdf_write_info_base nf90_put_var origin_time_varid" ) - call check( nf90_put_var(nciu%id, nciu%origin_rh_varid, self%info(j)%origin_rh(:), start=[1,idslot], count=[NDIM,1]), "netcdf_write_info_base nf90_put_var origin_rh_varid" ) - call check( nf90_put_var(nciu%id, nciu%origin_vh_varid, self%info(j)%origin_vh(:), start=[1,idslot], count=[NDIM,1]), "netcdf_write_info_base nf90_put_var origin_vh_varid" ) + call check( nf90_put_var(nc%id, nc%origin_type_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "netcdf_write_info_base nf90_put_var origin_type_varid" ) + call check( nf90_put_var(nc%id, nc%origin_time_varid, self%info(j)%origin_time, start=[idslot]), "netcdf_write_info_base nf90_put_var origin_time_varid" ) + call check( nf90_put_var(nc%id, nc%origin_rh_varid, self%info(j)%origin_rh(:), start=[1,idslot], count=[NDIM,1]), "netcdf_write_info_base nf90_put_var origin_rh_varid" ) + call check( nf90_put_var(nc%id, nc%origin_vh_varid, self%info(j)%origin_vh(:), start=[1,idslot], count=[NDIM,1]), "netcdf_write_info_base nf90_put_var origin_vh_varid" ) - call check( nf90_put_var(nciu%id, nciu%collision_id_varid, self%info(j)%collision_id, start=[idslot]), "netcdf_write_info_base nf90_put_var collision_id_varid" ) - call check( nf90_put_var(nciu%id, nciu%discard_time_varid, self%info(j)%discard_time, start=[idslot]), "netcdf_write_info_base nf90_put_var discard_time_varid" ) - call check( nf90_put_var(nciu%id, nciu%discard_rh_varid, self%info(j)%discard_rh(:), start=[1,idslot], count=[NDIM,1]), "netcdf_write_info_base nf90_put_var discard_rh_varid" ) - call check( nf90_put_var(nciu%id, nciu%discard_vh_varid, self%info(j)%discard_vh(:), start=[1,idslot], count=[NDIM,1]), "netcdf_write_info_base nf90_put_var discard_vh_varid" ) + call check( nf90_put_var(nc%id, nc%collision_id_varid, self%info(j)%collision_id, start=[idslot]), "netcdf_write_info_base nf90_put_var collision_id_varid" ) + call check( nf90_put_var(nc%id, nc%discard_time_varid, self%info(j)%discard_time, start=[idslot]), "netcdf_write_info_base nf90_put_var discard_time_varid" ) + call check( nf90_put_var(nc%id, nc%discard_rh_varid, self%info(j)%discard_rh(:), start=[1,idslot], count=[NDIM,1]), "netcdf_write_info_base nf90_put_var discard_rh_varid" ) + call check( nf90_put_var(nc%id, nc%discard_vh_varid, self%info(j)%discard_vh(:), start=[1,idslot], count=[NDIM,1]), "netcdf_write_info_base nf90_put_var discard_vh_varid" ) end if end do @@ -1199,36 +1200,36 @@ module subroutine netcdf_write_info_base(self, nciu, param) class is (swiftest_cb) idslot = self%id + 1 - call check( nf90_put_var(nciu%id, nciu%id_varid, self%id, start=[idslot]), "netcdf_write_info_base nf90_put_var cb id_varid" ) + call check( nf90_put_var(nc%id, nc%id_varid, self%id, start=[idslot]), "netcdf_write_info_base nf90_put_var cb id_varid" ) charstring = trim(adjustl(self%info%name)) - call check( nf90_put_var(nciu%id, nciu%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "netcdf_write_info_base nf90_put_var cb name_varid" ) + call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "netcdf_write_info_base nf90_put_var cb name_varid" ) charstring = trim(adjustl(self%info%particle_type)) - call check( nf90_put_var(nciu%id, nciu%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "netcdf_write_info_base nf90_put_var cb ptype_varid" ) + call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "netcdf_write_info_base nf90_put_var cb ptype_varid" ) if (param%lclose) then charstring = trim(adjustl(self%info%origin_type)) - call check( nf90_put_var(nciu%id, nciu%origin_type_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "netcdf_write_info_base nf90_put_var cb origin_type_varid" ) + call check( nf90_put_var(nc%id, nc%origin_type_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "netcdf_write_info_base nf90_put_var cb origin_type_varid" ) - call check( nf90_put_var(nciu%id, nciu%origin_time_varid, self%info%origin_time, start=[idslot]), "netcdf_write_info_base nf90_put_var cb origin_time_varid" ) - call check( nf90_put_var(nciu%id, nciu%origin_rh_varid, self%info%origin_rh(:), start=[1, idslot], count=[NDIM,1]), "netcdf_write_info_base nf90_put_var cb origin_rh_varid" ) - call check( nf90_put_var(nciu%id, nciu%origin_vh_varid, self%info%origin_vh(:), start=[1, idslot], count=[NDIM,1]), "netcdf_write_info_base nf90_put_var cb origin_vh_varid" ) + call check( nf90_put_var(nc%id, nc%origin_time_varid, self%info%origin_time, start=[idslot]), "netcdf_write_info_base nf90_put_var cb origin_time_varid" ) + call check( nf90_put_var(nc%id, nc%origin_rh_varid, self%info%origin_rh(:), start=[1, idslot], count=[NDIM,1]), "netcdf_write_info_base nf90_put_var cb origin_rh_varid" ) + call check( nf90_put_var(nc%id, nc%origin_vh_varid, self%info%origin_vh(:), start=[1, idslot], count=[NDIM,1]), "netcdf_write_info_base nf90_put_var cb origin_vh_varid" ) - call check( nf90_put_var(nciu%id, nciu%collision_id_varid, self%info%collision_id, start=[idslot]), "netcdf_write_info_base nf90_put_var cb collision_id_varid" ) - call check( nf90_put_var(nciu%id, nciu%discard_time_varid, self%info%discard_time, start=[idslot]), "netcdf_write_info_base nf90_put_var cb discard_time_varid" ) - call check( nf90_put_var(nciu%id, nciu%discard_rh_varid, self%info%discard_rh(:), start=[1, idslot], count=[NDIM,1]), "netcdf_write_info_base nf90_put_var cb discard_rh_varid" ) - call check( nf90_put_var(nciu%id, nciu%discard_vh_varid, self%info%discard_vh(:), start=[1, idslot], count=[NDIM,1]), "netcdf_write_info_base nf90_put_var cb discard_vh_varid" ) + call check( nf90_put_var(nc%id, nc%collision_id_varid, self%info%collision_id, start=[idslot]), "netcdf_write_info_base nf90_put_var cb collision_id_varid" ) + call check( nf90_put_var(nc%id, nc%discard_time_varid, self%info%discard_time, start=[idslot]), "netcdf_write_info_base nf90_put_var cb discard_time_varid" ) + call check( nf90_put_var(nc%id, nc%discard_rh_varid, self%info%discard_rh(:), start=[1, idslot], count=[NDIM,1]), "netcdf_write_info_base nf90_put_var cb discard_rh_varid" ) + call check( nf90_put_var(nc%id, nc%discard_vh_varid, self%info%discard_vh(:), start=[1, idslot], count=[NDIM,1]), "netcdf_write_info_base nf90_put_var cb discard_vh_varid" ) end if end select - call check( nf90_set_fill(nciu%id, old_mode, old_mode) ) + call check( nf90_set_fill(nc%id, old_mode, old_mode) ) return end subroutine netcdf_write_info_base - module subroutine netcdf_write_hdr_system(self, nciu, param) + module subroutine netcdf_write_hdr_system(self, nc, param) !! author: David A. Minton !! !! Writes header information (variables that change with time, but not particle id). @@ -1237,31 +1238,31 @@ module subroutine netcdf_write_hdr_system(self, nciu, param) implicit none ! Arguments class(swiftest_nbody_system), intent(in) :: self !! Swiftest nbody system object - class(netcdf_parameters), intent(inout) :: nciu !! Parameters used to for writing a NetCDF dataset to file + class(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 check( nf90_put_var(nciu%id, nciu%time_varid, self%t, start=[tslot]), "netcdf_write_hdr_system nf90_put_var time_varid" ) - call check( nf90_put_var(nciu%id, nciu%npl_varid, self%pl%nbody, start=[tslot]), "netcdf_write_hdr_system nf90_put_var npl_varid" ) - call check( nf90_put_var(nciu%id, nciu%ntp_varid, self%tp%nbody, start=[tslot]), "netcdf_write_hdr_system nf90_put_var ntp_varid" ) + call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "netcdf_write_hdr_system nf90_put_var time_varid" ) + call check( nf90_put_var(nc%id, nc%npl_varid, self%pl%nbody, start=[tslot]), "netcdf_write_hdr_system nf90_put_var npl_varid" ) + call check( nf90_put_var(nc%id, nc%ntp_varid, self%tp%nbody, start=[tslot]), "netcdf_write_hdr_system nf90_put_var ntp_varid" ) select type(pl => self%pl) class is (symba_pl) - call check( nf90_put_var(nciu%id, nciu%nplm_varid, pl%nplm, start=[tslot]), "netcdf_write_hdr_system nf90_put_var nplm_varid" ) + call check( nf90_put_var(nc%id, nc%nplm_varid, pl%nplm, start=[tslot]), "netcdf_write_hdr_system nf90_put_var nplm_varid" ) end select if (param%lenergy) then - call check( nf90_put_var(nciu%id, nciu%KE_orb_varid, self%ke_orbit, start=[tslot]), "netcdf_write_hdr_system nf90_put_var KE_orb_varid" ) - call check( nf90_put_var(nciu%id, nciu%KE_spin_varid, self%ke_spin, start=[tslot]), "netcdf_write_hdr_system nf90_put_var KE_spin_varid" ) - call check( nf90_put_var(nciu%id, nciu%PE_varid, self%pe, start=[tslot]), "netcdf_write_hdr_system nf90_put_var PE_varid" ) - call check( nf90_put_var(nciu%id, nciu%L_orb_varid, self%Lorbit(:), start=[1,tslot], count=[NDIM,1]), "netcdf_write_hdr_system nf90_put_var L_orb_varid" ) - call check( nf90_put_var(nciu%id, nciu%L_spin_varid, self%Lspin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_write_hdr_system nf90_put_var L_spin_varid" ) - call check( nf90_put_var(nciu%id, nciu%L_escape_varid, self%Lescape(:), start=[1,tslot], count=[NDIM,1]), "netcdf_write_hdr_system nf90_put_var L_escape_varid" ) - call check( nf90_put_var(nciu%id, nciu%Ecollisions_varid, self%Ecollisions, start=[tslot]), "netcdf_write_hdr_system nf90_put_var Ecollisions_varid" ) - call check( nf90_put_var(nciu%id, nciu%Euntracked_varid, self%Euntracked, start=[tslot]), "netcdf_write_hdr_system nf90_put_var Euntracked_varid" ) - call check( nf90_put_var(nciu%id, nciu%GMescape_varid, self%GMescape, start=[tslot]), "netcdf_write_hdr_system nf90_put_var GMescape_varid" ) + call check( nf90_put_var(nc%id, nc%KE_orb_varid, self%ke_orbit, start=[tslot]), "netcdf_write_hdr_system nf90_put_var KE_orb_varid" ) + call check( nf90_put_var(nc%id, nc%KE_spin_varid, self%ke_spin, start=[tslot]), "netcdf_write_hdr_system nf90_put_var KE_spin_varid" ) + call check( nf90_put_var(nc%id, nc%PE_varid, self%pe, start=[tslot]), "netcdf_write_hdr_system nf90_put_var PE_varid" ) + call check( nf90_put_var(nc%id, nc%L_orb_varid, self%Lorbit(:), start=[1,tslot], count=[NDIM,1]), "netcdf_write_hdr_system nf90_put_var L_orb_varid" ) + call check( nf90_put_var(nc%id, nc%L_spin_varid, self%Lspin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_write_hdr_system nf90_put_var L_spin_varid" ) + call check( nf90_put_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1,tslot], count=[NDIM,1]), "netcdf_write_hdr_system nf90_put_var L_escape_varid" ) + call check( nf90_put_var(nc%id, nc%Ecollisions_varid, self%Ecollisions, start=[tslot]), "netcdf_write_hdr_system nf90_put_var Ecollisions_varid" ) + call check( nf90_put_var(nc%id, nc%Euntracked_varid, self%Euntracked, start=[tslot]), "netcdf_write_hdr_system nf90_put_var Euntracked_varid" ) + call check( nf90_put_var(nc%id, nc%GMescape_varid, self%GMescape, start=[tslot]), "netcdf_write_hdr_system nf90_put_var GMescape_varid" ) end if return diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 4a3d98aed..ea6d3db23 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -68,8 +68,7 @@ module subroutine setup_construct_system(system, param) allocate(symba_pltpenc :: system%pltpenc_list) allocate(symba_plplenc :: system%plplenc_list) allocate(symba_plplenc :: system%plplcollision_list) - allocate(symba_pl :: system%snapshot%pl) - allocate(symba_tp :: system%snapshot%tp) + allocate(symba_encounter_storage :: system%encounter_history) end select case (RINGMOONS) write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled' @@ -93,7 +92,7 @@ module subroutine setup_finalize_system(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters associate(system => self) - call param%nciu%close() + call param%nc%close() end associate return diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index e58da2129..dd60f2e00 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -34,7 +34,7 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l lany_encounter = .false. if (self%nbody == 0) return - associate(pl => self, plplenc_list => system%plplenc_list, cb => system%cb, ienc_frame => system%ienc_frame) + associate(pl => self, plplenc_list => system%plplenc_list, cb => system%cb) npl = pl%nbody nplm = pl%nplm diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 9cfd8ba9a..0fceaa724 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -11,6 +11,173 @@ use swiftest contains + module subroutine symba_io_encounter_dump(self, param) + !! author: David A. Minton + !! + !! Dumps the time history of an encounter to file. + implicit none + ! Arguments + class(symba_encounter_storage(*)), intent(inout) :: self !! Encounter storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: i + + ! Most of this is just temporary test code just to get something working. Eventually this should get cleaned up. + + do i = 1, self%nframes + if (allocated(self%frame(i)%item)) then + select type(snapshot => self%frame(i)%item) + class is (symba_encounter_snapshot) + self%nc%ienc_frame = i + call snapshot%write_frame(self%nc,param) + end select + end if + end do + + + return + end subroutine symba_io_encounter_dump + + + module subroutine symba_io_encounter_initialize_output(self, param) + !! author: David A. Minton + !! + !! Initialize a NetCDF encounter file system. This is a simplified version of the main simulation output NetCDF file, but with fewer variables. + use, intrinsic :: ieee_arithmetic + use netcdf + implicit none + ! Arguments + class(symba_io_encounter_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: nvar, varid, vartype + real(DP) :: dfill + real(SP) :: sfill + logical :: fileExists + character(len=STRMAX) :: errmsg + integer(I4B) :: ndims + + + associate(nc => self) + dfill = ieee_value(dfill, IEEE_QUIET_NAN) + sfill = ieee_value(sfill, IEEE_QUIET_NAN) + + select case (param%out_type) + case("NETCDF_FLOAT") + self%out_type = NF90_FLOAT + case("NETCDF_DOUBLE") + self%out_type = NF90_DOUBLE + end select + + + ! Check if the file exists, and if it does, delete it + inquire(file=nc%enc_file, exist=fileExists) + if (fileExists) then + open(unit=LUN, file=nc%enc_file, status="old", err=667, iomsg=errmsg) + close(unit=LUN, status="delete") + end if + + call check( nf90_create(nc%enc_file, NF90_NETCDF4, nc%id), "symba_io_encounter_initialize_output nf90_create" ) + + ! Dimensions + call check( nf90_def_dim(nc%id, nc%time_dimname, NF90_UNLIMITED, nc%time_dimid), "symba_io_encounter_initialize_output nf90_def_dim time_dimid" ) ! Simulation time dimension + call check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "symba_io_encounter_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension + call check( nf90_def_dim(nc%id, nc%id_dimname, NF90_UNLIMITED, nc%id_dimid), "symba_io_encounter_initialize_output nf90_def_dim id_dimid" ) ! dimension to store particle id numbers + call check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "symba_io_encounter_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) + + ! Dimension coordinates + call check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, nc%time_dimid, nc%time_varid), "symba_io_encounter_initialize_output nf90_def_var time_varid" ) + call check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "symba_io_encounter_initialize_output nf90_def_var space_varid" ) + call check( nf90_def_var(nc%id, nc%id_dimname, NF90_INT, nc%id_dimid, nc%id_varid), "symba_io_encounter_initialize_output nf90_def_var id_varid" ) + + ! Variables + call check( nf90_def_var(nc%id, nc%name_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], nc%name_varid), "symba_io_encounter_initialize_output nf90_def_var name_varid" ) + call check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], nc%ptype_varid), "symba_io_encounter_initialize_output nf90_def_var ptype_varid" ) + call check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%rh_varid), "symba_io_encounter_initialize_output nf90_def_var rh_varid" ) + call check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%vh_varid), "symba_io_encounter_initialize_output nf90_def_var vh_varid" ) + call check( nf90_def_var(nc%id, nc%gmass_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%Gmass_varid), "symba_io_encounter_initialize_output nf90_def_var Gmass_varid" ) + if (param%lclose) then + call check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%radius_varid), "symba_io_encounter_initialize_output nf90_def_var radius_varid" ) + end if + if (param%lrotation) then + call check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%Ip_varid), "symba_io_encounter_initialize_output nf90_def_var Ip_varid" ) + call check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%rot_varid), "symba_io_encounter_initialize_output nf90_def_var rot_varid" ) + end if + + call check( nf90_inquire(nc%id, nVariables=nvar), "symba_io_encounter_initialize_output nf90_inquire nVariables" ) + do varid = 1, nvar + call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "symba_io_encounter_initialize_output nf90_inquire_variable" ) + select case(vartype) + case(NF90_INT) + call check( nf90_def_var_fill(nc%id, varid, 0, NF90_FILL_INT), "symba_io_encounter_initialize_output nf90_def_var_fill NF90_INT" ) + case(NF90_FLOAT) + call check( nf90_def_var_fill(nc%id, varid, 0, sfill), "symba_io_encounter_initialize_output nf90_def_var_fill NF90_FLOAT" ) + case(NF90_DOUBLE) + call check( nf90_def_var_fill(nc%id, varid, 0, dfill), "symba_io_encounter_initialize_output nf90_def_var_fill NF90_DOUBLE" ) + case(NF90_CHAR) + call check( nf90_def_var_fill(nc%id, varid, 0, 0), "symba_io_encounter_initialize_output nf90_def_var_fill NF90_CHAR" ) + end select + end do + + ! Take the file out of define mode + call check( nf90_enddef(nc%id), "symba_io_encounter_initialize_output nf90_enddef" ) + + ! Add in the space dimension coordinates + call check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "symba_io_encounter_initialize_output nf90_put_var space" ) + end associate + + return + + 667 continue + write(*,*) "Error creating encounter output file. " // trim(adjustl(errmsg)) + call util_exit(FAILURE) + end subroutine symba_io_encounter_initialize_output + + + module subroutine symba_io_encounter_write_frame(self, nc, param) + !! author: David A. Minton + !! + !! Write a frame of output of an encounter trajectory. + use netcdf + implicit none + ! Arguments + class(symba_encounter_snapshot), intent(in) :: self !! Swiftest encounter structure + class(symba_io_encounter_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: i, tslot, idslot, old_mode, n + character(len=NAMELEN) :: charstring + + tslot = self%tslot + call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "symba_io_encounter_write_frame nf90_set_fill" ) + select type(pl => self%pl) + class is (symba_pl) + n = size(pl%id(:)) + do i = 1, n + idslot = pl%id(i) + call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "symba_io_encounter_write_frame nf90_put_var time_varid" ) + call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[idslot]), "symba_io_encounter_write_frame nf90_put_var id_varid" ) + call check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var rh_varid" ) + call check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var vh_varid" ) + call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[idslot, tslot]), "symba_io_encounter_write_frame nf90_put_var body Gmass_varid" ) + if (param%lclose) call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[idslot, tslot]), "symba_io_encounter_write_frame nf90_put_var body radius_varid" ) + if (param%lrotation) then + call check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var body Ip_varid" ) + call check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1,idslot, tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var body rotx_varid" ) + end if + charstring = trim(adjustl(pl%info(i)%name)) + call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "symba_io_encounter_write_frame nf90_put_var name_varid" ) + charstring = trim(adjustl(pl%info(i)%particle_type)) + call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "symba_io_encounter_write_frame nf90_put_var particle_type_varid" ) + end do + end select + + call check( nf90_set_fill(nc%id, old_mode, old_mode) ) + + return + end subroutine symba_io_encounter_write_frame + + module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! @@ -66,6 +233,9 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms case ("ENCOUNTER_SAVE") call io_toupper(param_value) read(param_value, *) param%encounter_save + case ("FRAGMENTATION_SAVE") + call io_toupper(param_value) + read(param_value, *) param%fragmentation_save case("SEED") read(param_value, *) nseeds_from_file ! Because the number of seeds can vary between compilers/systems, we need to make sure we can handle cases in which the input file has a different @@ -116,12 +286,21 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms ! All reporting of collision information in SyMBA (including mergers) is now recorded in the Fraggle logfile call io_log_start(param, FRAGGLE_LOG_OUT, "Fraggle logfile") - if ((param%encounter_save /= "NONE") .and. (param%encounter_save /= "ALL") .and. (param%encounter_save /= "FRAGMENTATION")) then + if ((param%encounter_save /= "NONE") .and. (param%encounter_save /= "TRAJECTORY") .and. (param%encounter_save /= "CLOSEST")) then write(iomsg,*) 'Invalid encounter_save parameter: ',trim(adjustl(param%out_type)) - write(iomsg,*) 'Valid options are NONE, ALL, or FRAGMENTATION' + write(iomsg,*) 'Valid options are NONE, TRAJECTORY, or CLOSEST' + iostat = -1 + return + end if + + if ((param%fragmentation_save /= "NONE") .and. (param%fragmentation_save /= "TRAJECTORY") .and. (param%fragmentation_save /= "CLOSEST")) then + write(iomsg,*) 'Invalid fragmentation_save parameter: ',trim(adjustl(param%out_type)) + write(iomsg,*) 'Valid options are NONE, TRAJECTORY, or CLOSEST' iostat = -1 return end if + param%lencounter_save = (param%encounter_save == "TRAJECTORY") .or. (param%encounter_save == "CLOSEST") .or. & + (param%fragmentation_save == "TRAJECTORY") .or. (param%fragmentation_save == "CLOSEST") ! Call the base method (which also prints the contents to screen) call io_param_reader(param, unit, iotype, v_list, iostat, iomsg) @@ -176,6 +355,52 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms end subroutine symba_io_param_writer + module subroutine symba_io_start_encounter(self, param, t) + !! author: David A. Minton + !! + !! Initializes the new encounter and/or fragmentation history + implicit none + ! Arguments + class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current simulation time + + if (.not. allocated(self%encounter_history)) allocate(symba_encounter_storage :: self%encounter_history) + call self%encounter_history%reset() + + ! Empty out the time slot array for the next pass + self%encounter_history%tvals(:) = huge(1.0_DP) + + ! Take the snapshot at the start of the encounter + call self%snapshot(param, t) + + return + end subroutine symba_io_start_encounter + + + module subroutine symba_io_stop_encounter(self, param, t) + !! author: David A. Minton + !! + !! Saves the encounter and/or fragmentation data to file(s) + implicit none + ! Arguments + class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current simulation time + ! Internals + !character(STRMAX) + + ! Create and save the output file for this encounter + write(self%encounter_history%nc%enc_file, '("encounter_",I0.6,".nc")') param%iloop + call self%encounter_history%nc%initialize(param) + call self%encounter_history%dump(param) + call self%encounter_history%nc%close() + call self%encounter_history%reset() + + return + end subroutine symba_io_stop_encounter + + module subroutine symba_io_write_discard(self, param) !! author: David A. Minton !! @@ -187,12 +412,12 @@ module subroutine symba_io_write_discard(self, param) associate(pl => self%pl, npl => self%pl%nbody, pl_adds => self%pl_adds) - if (self%tp_discards%nbody > 0) call self%tp_discards%write_info(param%nciu, param) + if (self%tp_discards%nbody > 0) call self%tp_discards%write_info(param%nc, param) select type(pl_discards => self%pl_discards) class is (symba_merger) if (pl_discards%nbody == 0) return - call pl_discards%write_info(param%nciu, param) + call pl_discards%write_info(param%nc, param) end select end associate diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 9e9d2de41..e727ed9f3 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -37,8 +37,9 @@ module subroutine symba_step_system(self, param, t, dt) call self%reset(param) lencounter = pl%encounter_check(param, self, dt, 0) .or. tp%encounter_check(param, self, dt, 0) if (lencounter) then + if (param%lencounter_save) call self%start_encounter(param, t) call self%interp(param, t, dt) - !call self%encounter_history%dump(param) + if (param%lencounter_save) call self%stop_encounter(param, t+dt) else self%irec = -1 call helio_step_system(self, param, t, dt) @@ -180,72 +181,74 @@ recursive module subroutine symba_step_recur_system(self, param, t, ireci) logical :: lencounter, lplpl_collision, lpltp_collision associate(system => self, plplenc_list => self%plplenc_list, pltpenc_list => self%pltpenc_list) - select type(pl => self%pl) - class is (symba_pl) - select type(tp => self%tp) - class is (symba_tp) - system%irec = ireci - dtl = param%dt / (NTENC**ireci) - dth = 0.5_DP * dtl - IF (dtl / param%dt < VSMALL) THEN - write(*, *) "SWIFTEST Warning:" - write(*, *) " In symba_step_recur_system, local time step is too small" - write(*, *) " Roundoff error will be important!" - call util_exit(FAILURE) - END IF - irecp = ireci + 1 - if (ireci == 0) then - nloops = 1 - else - nloops = NTENC - end if - do j = 1, nloops - lencounter = plplenc_list%encounter_check(param, system, dtl, irecp) & - .or. pltpenc_list%encounter_check(param, system, dtl, irecp) - - call plplenc_list%kick(system, dth, irecp, 1) - call pltpenc_list%kick(system, dth, irecp, 1) - if (ireci /= 0) then - call plplenc_list%kick(system, dth, irecp, -1) - call pltpenc_list%kick(system, dth, irecp, -1) - end if - - if (param%lgr) then - call pl%gr_pos_kick(system, param, dth) - call tp%gr_pos_kick(system, param, dth) - end if - - call pl%drift(system, param, dtl) - call tp%drift(system, param, dtl) - - !call system% - - if (lencounter) call system%recursive_step(param, t+dth,irecp) + select type(param) + class is (symba_parameters) + select type(pl => self%pl) + class is (symba_pl) + select type(tp => self%tp) + class is (symba_tp) system%irec = ireci - - if (param%lgr) then - call pl%gr_pos_kick(system, param, dth) - call tp%gr_pos_kick(system, param, dth) - end if - - call plplenc_list%kick(system, dth, irecp, 1) - call pltpenc_list%kick(system, dth, irecp, 1) - if (ireci /= 0) then - call plplenc_list%kick(system, dth, irecp, -1) - call pltpenc_list%kick(system, dth, irecp, -1) - end if - - if (param%lclose) then - lplpl_collision = plplenc_list%collision_check(system, param, t+dtl, dtl, ireci) - lpltp_collision = pltpenc_list%collision_check(system, param, t+dtl, dtl, ireci) - - if (lplpl_collision) call plplenc_list%resolve_collision(system, param, t+dtl, dtl, ireci) - if (lpltp_collision) call pltpenc_list%resolve_collision(system, param, t+dtl, dtl, ireci) + dtl = param%dt / (NTENC**ireci) + dth = 0.5_DP * dtl + IF (dtl / param%dt < VSMALL) THEN + write(*, *) "SWIFTEST Warning:" + write(*, *) " In symba_step_recur_system, local time step is too small" + write(*, *) " Roundoff error will be important!" + call util_exit(FAILURE) + END IF + irecp = ireci + 1 + if (ireci == 0) then + nloops = 1 + else + nloops = NTENC end if - - call self%set_recur_levels(ireci) - - end do + do j = 1, nloops + lencounter = plplenc_list%encounter_check(param, system, dtl, irecp) & + .or. pltpenc_list%encounter_check(param, system, dtl, irecp) + + call plplenc_list%kick(system, dth, irecp, 1) + call pltpenc_list%kick(system, dth, irecp, 1) + if (ireci /= 0) then + call plplenc_list%kick(system, dth, irecp, -1) + call pltpenc_list%kick(system, dth, irecp, -1) + end if + + if (param%lgr) then + call pl%gr_pos_kick(system, param, dth) + call tp%gr_pos_kick(system, param, dth) + end if + + call pl%drift(system, param, dtl) + call tp%drift(system, param, dtl) + + if (lencounter) call system%recursive_step(param, t+dth,irecp) + system%irec = ireci + + if (param%lgr) then + call pl%gr_pos_kick(system, param, dth) + call tp%gr_pos_kick(system, param, dth) + end if + + call plplenc_list%kick(system, dth, irecp, 1) + call pltpenc_list%kick(system, dth, irecp, 1) + if (ireci /= 0) then + call plplenc_list%kick(system, dth, irecp, -1) + call pltpenc_list%kick(system, dth, irecp, -1) + end if + + if (param%lclose) then + lplpl_collision = plplenc_list%collision_check(system, param, t+dtl, dtl, ireci) + lpltp_collision = pltpenc_list%collision_check(system, param, t+dtl, dtl, ireci) + + if (lplpl_collision) call plplenc_list%resolve_collision(system, param, t+dtl, dtl, ireci) + if (lpltp_collision) call pltpenc_list%resolve_collision(system, param, t+dtl, dtl, ireci) + end if + if (param%lencounter_save) call system%snapshot(param, t+dtl) + + call self%set_recur_levels(ireci) + + end do + end select end select end select end associate @@ -278,8 +281,6 @@ module subroutine symba_step_reset_system(self, param) nenc_old = system%plplenc_list%nenc call system%plplenc_list%setup(0_I8B) call system%plplcollision_list%setup(0_I8B) - system%ienc_frame = 0 - if (allocated(system%encounter_history)) deallocate(system%encounter_history) if (npl > 0) then pl%lcollision(1:npl) = .false. call pl%reset_kinship([(i, i=1, npl)]) @@ -313,6 +314,7 @@ module subroutine symba_step_reset_system(self, param) tp%lfirst = param%lfirstkick pl%lfirst = param%lfirstkick + end associate end select end select diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 63ba3224b..20d295b9c 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -498,19 +498,29 @@ module subroutine symba_util_final_system(self) return end subroutine symba_util_final_system + module subroutine symba_util_final_encounter_snapshot(self) + !! author: David A. Minton + !! + !! Finalize the SyMBA encounter system snapshot object - deallocates all allocatables + implicit none + type(symba_encounter_snapshot), intent(inout) :: self !! SyMBA nbody system object + + call self%dealloc() - module subroutine symba_util_final_snapshot(self) + return + end subroutine symba_util_final_encounter_snapshot + + + module subroutine symba_util_final_encounter_storage(self) !! author: David A. Minton !! !! Finalize the SyMBA nbody system object - deallocates all allocatables implicit none ! Argument - type(symba_system_snapshot), intent(inout) :: self !! SyMBA nbody system object - - call self%dealloc() + type(symba_encounter_storage(*)), intent(inout) :: self !! SyMBA nbody system object return - end subroutine symba_util_final_snapshot + end subroutine symba_util_final_encounter_storage module subroutine symba_util_final_tp(self) @@ -694,7 +704,7 @@ module subroutine symba_util_rearray_pl(self, system, param) end where end select - call pl%write_info(param%nciu, param) + call pl%write_info(param%nc, param) deallocate(ldump_mask) ! Reindex the new list of bodies @@ -908,34 +918,39 @@ module subroutine symba_util_resize_storage(self, nnew) class(symba_nbody_system), intent(inout) :: self !! Swiftest encounter list integer(I4B), intent(in) :: nnew !! New size of list needed ! Internals - type(encounter_storage(nframes=:)), allocatable :: tmp - integer(I4B) :: i, nold + type(symba_encounter_storage(nframes=:)), allocatable :: tmp + integer(I4B) :: i, nold, nbig, iframe_old = 0 logical :: lmalloc lmalloc = allocated(self%encounter_history) if (lmalloc) then nold = self%encounter_history%nframes + iframe_old = self%encounter_history%iframe else nold = 0 end if if (nnew > nold) then - allocate(encounter_storage(8 * nnew) :: tmp) + nbig = nold + do while (nbig < nnew) + nbig = nbig * 2 + end do + allocate(symba_encounter_storage(nbig) :: tmp) if (lmalloc) then do i = 1, nold - if (allocated(self%encounter_history%frame(i)%item)) tmp%frame(i) = self%encounter_history%frame(i)%item + if (allocated(self%encounter_history%frame(i)%item)) call move_alloc(self%encounter_history%frame(i)%item, tmp%frame(i)%item) end do deallocate(self%encounter_history) end if call move_alloc(tmp,self%encounter_history) + self%encounter_history%iframe = iframe_old end if return end subroutine symba_util_resize_storage - module subroutine symba_util_resize_tp(self, nnew) !! author: David A. Minton !! @@ -1242,10 +1257,10 @@ module subroutine symba_util_spill_encounter_list(self, discards, lspill_list, l !! Note: Because the symba_plplenc currently does not contain any additional variable components, this method can recieve it as an input as well. implicit none ! Arguments - class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter list - class(encounter_list), intent(inout) :: discards !! Discarded object - logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards - logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list + class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter list + class(encounter_list), intent(inout) :: discards !! Discarded object + logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards + logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list associate(keeps => self) select type(discards) @@ -1294,23 +1309,112 @@ module subroutine symba_util_spill_tp(self, discards, lspill_list, ldestructive) end subroutine symba_util_spill_tp - module subroutine symba_util_take_system_snapshot(self, system, param, t) + module subroutine symba_util_take_encounter_snapshot(self, param, t) !! author: David A. Minton !! !! Takes a minimal snapshot of the state of the system during an encounter so that the trajectories - !! Can be played back through the encounter + !! can be played back through the encounter implicit none ! Internals - class(symba_system_snapshot), intent(inout) :: self !! SyMBA nbody system snapshot object - class(symba_nbody_system), intent(in) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! current time + class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! current time ! Arguments - logical, dimension(:), allocatable :: lmask + type(symba_encounter_snapshot), allocatable :: snapshot + integer(I4B) :: i, npl_snap, ntp_snap + + associate(npl => self%pl%nbody, ntp => self%tp%nbody) + + allocate(symba_encounter_snapshot :: snapshot) - !if (system%pl) + if (npl > 0) allocate(symba_pl :: snapshot%pl) + if (ntp > 0) allocate(symba_tp :: snapshot%tp) + if (npl + ntp == 0) return + npl_snap = npl + ntp_snap = ntp + + select type (pl => self%pl) + class is (symba_pl) + select type (tp => self%tp) + class is (symba_tp) + if (npl > 0) then + pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE .and. pl%levelg(1:npl) == self%irec + npl_snap = count(pl%lmask(1:npl)) + end if + if (ntp > 0) then + tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE .and. tp%levelg(1:ntp) == self%irec + ntp_snap = count(tp%lmask(1:ntp)) + end if + + ! Take snapshot of the currently encountering massive bodies + if (npl_snap > 0) then + allocate(snapshot%pl%id(npl_snap)) + allocate(snapshot%pl%info(npl_snap)) + allocate(snapshot%pl%Gmass(npl_snap)) + snapshot%pl%id(:) = pack(pl%id(1:npl), pl%lmask(1:npl)) + snapshot%pl%info(:) = pack(pl%info(1:npl), pl%lmask(1:npl)) + snapshot%pl%Gmass(:) = pack(pl%Gmass(1:npl), pl%lmask(1:npl)) + allocate(snapshot%pl%rh(NDIM,npl_snap)) + allocate(snapshot%pl%vh(NDIM,npl_snap)) + do i = 1, NDIM + snapshot%pl%rh(i,:) = pack(pl%rh(i,1:npl), pl%lmask(1:npl)) + snapshot%pl%vh(i,:) = pack(pl%vb(i,1:npl), pl%lmask(1:npl)) + end do + if (param%lclose) then + allocate(snapshot%pl%radius(npl_snap)) + snapshot%pl%radius(:) = pack(pl%radius(1:npl), pl%lmask(1:npl)) + end if + + if (param%lrotation) then + allocate(snapshot%pl%Ip(NDIM,npl_snap)) + allocate(snapshot%pl%rot(NDIM,npl_snap)) + do i = 1, NDIM + snapshot%pl%Ip(i,:) = pack(pl%Ip(i,1:npl), pl%lmask(1:npl)) + snapshot%pl%rot(i,:) = pack(pl%rot(i,1:npl), pl%lmask(1:npl)) + end do + end if + call snapshot%pl%sort("id", ascending=.true.) + end if + + ! Take snapshot of the currently encountering test particles + if (ntp_snap > 0) then + allocate(snapshot%tp%id(ntp_snap)) + allocate(snapshot%tp%info(ntp_snap)) + snapshot%tp%id(:) = pack(tp%id(1:ntp), tp%lmask(1:ntp)) + snapshot%tp%info(:) = pack(tp%info(1:ntp), tp%lmask(1:ntp)) + allocate(snapshot%tp%rh(NDIM,ntp_snap)) + allocate(snapshot%tp%vh(NDIM,ntp_snap)) + do i = 1, NDIM + snapshot%tp%rh(i,:) = pack(tp%rh(i,1:ntp), tp%lmask(1:ntp)) + snapshot%tp%vh(i,:) = pack(tp%vh(i,1:ntp), tp%lmask(1:ntp)) + end do + end if + + if (npl_snap + ntp_snap == 0) return + + snapshot%t = t + + ! Save the snapshot + self%encounter_history%iframe = self%encounter_history%iframe + 1 + call self%resize_storage(self%encounter_history%iframe) + + ! Find out which time slot this belongs in by searching for an existing slot + ! with the same value of time or the first available one + do i = 1, self%encounter_history%nframes + if (t <= self%encounter_history%tvals(i)) then + snapshot%tslot = i + self%encounter_history%tvals(i) = t + exit + end if + end do + + self%encounter_history%frame(self%encounter_history%iframe) = snapshot + + end select + end select + end associate return - end subroutine symba_util_take_system_snapshot + end subroutine symba_util_take_encounter_snapshot end submodule s_symba_util diff --git a/src/util/util_reset.f90 b/src/util/util_reset.f90 new file mode 100644 index 000000000..569846a68 --- /dev/null +++ b/src/util/util_reset.f90 @@ -0,0 +1,32 @@ +!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh +!! This file is part of Swiftest. +!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +!! You should have received a copy of the GNU General Public License along with Swiftest. +!! If not, see: https://www.gnu.org/licenses. + +submodule (swiftest_classes) s_util_reset + use swiftest +contains + + module subroutine util_reset_storage(self) + !! author: David A. Minton + !! + !! Resets a storage object by deallocating all items and resetting the frame counter to 0 + implicit none + ! Arguments + class(swiftest_storage(*)), intent(inout) :: self !! Swiftest storage object + ! Internals + integer(I4B) :: i + + do i = 1, self%nframes + if (allocated(self%frame(i)%item)) deallocate(self%frame(i)%item) + end do + self%iframe = 0 + + return + end subroutine util_reset_storage + +end submodule s_util_reset \ No newline at end of file