diff --git a/Makefile.Defines b/Makefile.Defines index 07126f842..820ad6d7d 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -65,8 +65,8 @@ GPAR = -fopenmp -ftree-parallelize-loops=4 GMEM = -fsanitize=undefined -fsanitize=address -fsanitize=leak GWARNINGS = -Wall -Warray-bounds -Wimplicit-interface -Wextra -Warray-temporaries -FFLAGS = $(IDEBUG) $(HEAPARR) -#FFLAGS = -init=snan,arrays -no-wrap-margin -O3 $(STRICTREAL) $(SIMDVEC) $(PAR) +#FFLAGS = $(IDEBUG) $(HEAPARR) +FFLAGS = -init=snan,arrays -no-wrap-margin -O3 $(STRICTREAL) $(SIMDVEC) $(PAR) FORTRAN = ifort #AR = xiar diff --git a/examples/symba_energy_momentum/collision_movie.py b/examples/symba_energy_momentum/collision_movie.py new file mode 100755 index 000000000..3fd3b6a86 --- /dev/null +++ b/examples/symba_energy_momentum/collision_movie.py @@ -0,0 +1,316 @@ +#!/usr/bin/env python3 +import swiftest +import numpy as np +import matplotlib.pyplot as plt +from matplotlib import animation +import matplotlib.collections as clt +from scipy.spatial.transform import Rotation as R + +xmin = -20.0 +xmax = 20.0 +ymin = -20.0 +ymax = 20.0 + +cases = ['supercat_head', 'supercat_off', 'disruption_head', 'disruption_off'] + +def scale_sim(ds): + + dsscale = ds.where(ds.id > 0, drop=True) # Remove the central body + + GMtot = dsscale['GMass'].sum(skipna=True, dim="id").isel(time=0) + rscale = ds['Radius'].sel(id=1, time=0) + dsscale['Radius'] /= rscale + + dsscale['radmarker'] = dsscale['Radius'].fillna(0) + + dsscale['px'] /= rscale + dsscale['py'] /= rscale + dsscale['pz'] /= rscale + + mpx = dsscale['GMass'] * dsscale['px'] + mpy = dsscale['GMass'] * dsscale['py'] + mpz = dsscale['GMass'] * dsscale['pz'] + xbsys = mpx.sum(skipna=True, dim="id") / GMtot + ybsys = mpy.sum(skipna=True, dim="id") / GMtot + zbsys = mpz.sum(skipna=True, dim="id") / GMtot + + mvx = dsscale['GMass'] * dsscale['vx'] + mvy = dsscale['GMass'] * dsscale['vy'] + mvz = dsscale['GMass'] * dsscale['vz'] + vxbsys = mvx.sum(skipna=True, dim="id") / GMtot + vybsys = mvy.sum(skipna=True, dim="id") / GMtot + vzbsys = mvz.sum(skipna=True, dim="id") / GMtot + + dsscale['pxb'] = dsscale['px'] - xbsys + dsscale['pyb'] = dsscale['py'] - ybsys + dsscale['pzb'] = dsscale['pz'] - zbsys + + dsscale['vxb'] = dsscale['vx'] - vxbsys + dsscale['vyb'] = dsscale['vy'] - vybsys + dsscale['vzb'] = dsscale['vz'] - vzbsys + + return dsscale + +class UpdatablePatchCollection(clt.PatchCollection): + def __init__(self, patches, *args, **kwargs): + self.patches = patches + clt.PatchCollection.__init__(self, patches, *args, **kwargs) + + def get_paths(self): + self.set_paths(self.patches) + return self._paths + +class AnimatedScatter(object): + """An animated scatter plot using matplotlib.animations.FuncAnimation.""" + + def __init__(self, ds, param): + + frame = 0 + nframes = ds['time'].size + self.ds = scale_sim(ds) + self.param = param + self.rot_angle = {} + + self.clist = {'Initial conditions' : 'xkcd:windows blue', + 'Disruption' : 'xkcd:baby poop', + 'Supercatastrophic' : 'xkcd:shocking pink', + 'Hit and run fragment' : 'xkcd:blue with a hint of purple', + 'Central body' : 'xkcd:almost black'} + + self.stream = self.data_stream(frame) + # Setup the figure and axes... + fig = plt.figure(figsize=(8,8), dpi=300) + plt.tight_layout(pad=0) + # set up the figure + self.ax = plt.Axes(fig, [0., 0., 1., 1.]) + self.ax.set_xlim(xmin, xmax) + self.ax.set_ylim(ymin, ymax) + self.ax.set_axis_off() + self.ax.set_aspect(1) + self.ax.get_xaxis().set_visible(False) + self.ax.get_yaxis().set_visible(False) + fig.add_axes(self.ax) + + # Then setup FuncAnimation. + self.ani = animation.FuncAnimation(fig, self.update, interval=1, frames=nframes, + init_func=self.setup_plot, blit=False) + self.ani.save(animfile, fps=60, dpi=300, + extra_args=['-vcodec', 'libx264']) + + def plot_pl_circles(self, pl, radmarker): + patches = [] + for i in range(pl.shape[0]): + s = plt.Circle((pl[i, 0], pl[i, 1]), radmarker[i]) + patches.append(s) + return patches + + def vec_props(self, c): + arrowprops = { + 'arrowstyle': '<|-', + 'mutation_scale': 20, + 'connectionstyle': 'arc3', + } + + arrow_args = { + 'xycoords': 'data', + 'textcoords': 'data', + 'arrowprops': arrowprops, + 'annotation_clip': True, + 'zorder': 100, + 'animated' : True + } + aarg = arrow_args.copy() + aprop = arrowprops.copy() + aprop['color'] = c + aarg['arrowprops'] = aprop + aarg['color'] = c + return aarg + + def plot_pl_vectors(self, pl, cval, r): + varrowend, varrowtip = self.velocity_vectors(pl, r) + arrows = [] + for i in range(pl.shape[0]): + aarg = self.vec_props(cval[i]) + a = self.ax.annotate("",xy=varrowend[i],xytext=varrowtip[i], **aarg) + arrows.append(a) + return arrows + + def plot_pl_spins(self, pl, id, cval, len): + sarrowend, sarrowtip = self.spin_arrows(pl, id, len) + arrows = [] + for i in range(pl.shape[0]): + aarg = self.vec_props(cval[i]) + aarg['arrowprops']['mutation_scale'] = 5 + aarg['arrowprops']['arrowstyle'] = "simple" + a = self.ax.annotate("",xy=sarrowend[i],xytext=sarrowtip[i], **aarg) + arrows.append(a) + return arrows + + def origin_to_color(self, origin): + cval = [] + for o in origin: + c = self.clist[o] + cval.append(c) + + return cval + + def velocity_vectors(self, pl, r): + px = pl[:, 0] + py = pl[:, 1] + vx = pl[:, 2] + vy = pl[:, 3] + vmag = np.sqrt(vx ** 2 + vy ** 2) + ux = np.zeros_like(vx) + uy = np.zeros_like(vx) + goodv = vmag > 0.0 + ux[goodv] = vx[goodv] / vmag[goodv] + uy[goodv] = vy[goodv] / vmag[goodv] + varrowend = [] + varrowtip = [] + for i in range(pl.shape[0]): + vend = (px[i], py[i]) + vtip = (px[i] + vx[i] * self.v_length, py[i] + vy[i] * self.v_length) + varrowend.append(vend) + varrowtip.append(vtip) + return varrowend, varrowtip + + def spin_arrows(self, pl, id, len): + px = pl[:, 0] + py = pl[:, 1] + sarrowend = [] + sarrowtip = [] + for i in range(pl.shape[0]): + endrel = np.array([0.0, len[i], 0.0]) + tiprel = np.array([0.0, -len[i], 0.0]) + r = R.from_rotvec(self.rot_angle[id[i]]) + endrel = r.apply(endrel) + tiprel = r.apply(tiprel) + send = (px[i] + endrel[0], py[i] + endrel[1]) + stip = (px[i] + tiprel[0], py[i] + tiprel[1]) + sarrowend.append(send) + sarrowtip.append(stip) + return sarrowend, sarrowtip + + def setup_plot(self): + # First frame + """Initial drawing of the scatter plot.""" + t, name, GMass, Radius, npl, pl, radmarker, origin = next(self.data_stream(0)) + + cval = self.origin_to_color(origin) + + # Scale markers to the size of the system + self.v_length = 0.50 # Length of arrow as fraction of velocity + + self.ax.margins(x=1, y=1) + self.ax.set_xlabel('x distance / ($R_1 + R_2$)', fontsize='16', labelpad=1) + self.ax.set_ylabel('y distance / ($R_1 + R_2$)', fontsize='16', labelpad=1) + + self.title = self.ax.text(0.50, 0.90, "", bbox={'facecolor': 'w', 'pad': 5}, transform=self.ax.transAxes, + ha="center", zorder=1000) + + self.title.set_text(titletext) + self.patches = self.plot_pl_circles(pl, radmarker) + + self.collection = UpdatablePatchCollection(self.patches, color=cval, alpha=0.5, zorder=50) + self.ax.add_collection(self.collection) + #self.varrows = self.plot_pl_vectors(pl, cval, radmarker) + self.sarrows = self.plot_pl_spins(pl, name, cval, radmarker) + + return self.collection, self.sarrows + + def update(self,frame): + """Update the scatter plot.""" + t, name, GMass, Radius, npl, pl, radmarker, origin = next(self.data_stream(frame)) + cval = self.origin_to_color(origin) + varrowend, varrowtip = self.velocity_vectors(pl, radmarker) + sarrowend, sarrowtip = self.spin_arrows(pl, name, radmarker) + for i, p in enumerate(self.patches): + p.set_center((pl[i, 0], pl[i,1])) + p.set_radius(radmarker[i]) + p.set_color(cval[i]) + #self.varrows[i].set_position(varrowtip[i]) + #self.varrows[i].xy = varrowend[i] + self.sarrows[i].set_position(sarrowtip[i]) + self.sarrows[i].xy = sarrowend[i] + + self.collection.set_paths(self.patches) + return self.collection, self.sarrows + + def data_stream(self, frame=0): + while True: + d = self.ds.isel(time=frame) + Radius = d['radmarker'].values + GMass = d['GMass'].values + x = d['pxb'].values + y = d['pyb'].values + vx = d['vxb'].values + vy = d['vyb'].values + name = d['id'].values + npl = d.id.count().values + id = d['id'].values + rotx = d['rot_x'].values + roty = d['rot_y'].values + rotz = d['rot_z'].values + + radmarker = d['radmarker'].values + origin = d['origin_type'].values + + t = self.ds.coords['time'].values[frame] + self.mask = np.logical_not(np.isnan(x)) + + x = np.nan_to_num(x, copy=False) + y = np.nan_to_num(y, copy=False) + vx = np.nan_to_num(vx, copy=False) + vy = np.nan_to_num(vy, copy=False) + radmarker = np.nan_to_num(radmarker, copy=False) + GMass = np.nan_to_num(GMass, copy=False) + Radius = np.nan_to_num(Radius, copy=False) + rotx = np.nan_to_num(rotx, copy=False) + roty = np.nan_to_num(roty, copy=False) + rotz = np.nan_to_num(rotz, copy=False) + rotvec = np.array([rotx, roty, rotz]) + self.rotvec = dict(zip(id, zip(*rotvec))) + + if frame == 0: + tmp = np.zeros_like(rotvec) + self.rot_angle = dict(zip(id, zip(*tmp))) + else: + t0 = self.ds.coords['time'].values[frame-1] + dt = t - t0 + idxactive = np.arange(id.size)[self.mask] + for i in id[idxactive]: + self.rot_angle[i] = self.rot_angle[i] + dt * np.array(self.rotvec[i]) + frame += 1 + yield t, name, GMass, Radius, npl, np.c_[x, y, vx, vy], radmarker, origin + +for case in cases: + if case == 'supercat_off': + animfile = 'movies/supercat_off_axis.mp4' + titletext = "Supercatastrophic - Off Axis" + paramfile = 'param.supercatastrophic_off_axis.in' + elif case == 'supercat_head': + animfile = 'movies/supercat_headon.mp4' + titletext = "Supercatastrophic - Head on" + paramfile = 'param.supercatastrophic_headon.in' + elif case == 'disruption_off': + animfile = 'movies/disruption_off_axis.mp4' + titletext = "Disruption - Off Axis" + paramfile = 'param.disruption_off_axis.in' + elif case == 'disruption_head': + animfile = 'movies/disruption_headon.mp4' + titletext = "Disruption- Head on" + paramfile = 'param.disruption_headon.in' + elif case == 'merger': + animfile = 'movies/merger.mp4' + titletext = "Merger" + paramfile = 'param.merger.in' + else: + print(f'{case} is an unknown case') + exit(-1) + sim = swiftest.Simulation(param_file=paramfile) + sim.bin2xr() + ds = sim.ds + print('Making animation') + anim = AnimatedScatter(ds,sim.param) + print('Animation finished') + plt.close(fig='all') diff --git a/examples/symba_energy_momentum/disruption_headon.in b/examples/symba_energy_momentum/disruption_headon.in index e1a5316bc..bc91bbdd0 100644 --- a/examples/symba_energy_momentum/disruption_headon.in +++ b/examples/symba_energy_momentum/disruption_headon.in @@ -1,11 +1,11 @@ 2 -2 1e-07 0.0009 +1 1e-07 0.0009 7e-06 1.0 -4.20E-05 0.0 0.00 6.28 0.0 0.4 0.4 0.4 !Ip 0.0 0.0 6.0e4 !rot -3 7e-10 0.0004 +2 7e-10 0.0004 3.25e-06 1.0 4.20E-05 0.0 0.00 -6.28 0.0 diff --git a/examples/symba_energy_momentum/disruption_off_axis.in b/examples/symba_energy_momentum/disruption_off_axis.in index b6bc29c26..792bb3a4a 100644 --- a/examples/symba_energy_momentum/disruption_off_axis.in +++ b/examples/symba_energy_momentum/disruption_off_axis.in @@ -1,11 +1,11 @@ 2 -2 1e-07 0.0009 +1 1e-07 0.0009 7e-06 1.0 -4.20E-05 0.0 0.00 6.28 0.0 0.4 0.4 0.4 !Ip 0.0 0.0 6.0e4 !rot -3 7e-10 0.0004 +2 7e-10 0.0004 3.25e-06 1.0 4.20E-05 0.0 -0.80 -6.28 0.0 diff --git a/examples/symba_energy_momentum/escape.in b/examples/symba_energy_momentum/escape.in index b8308af87..911cfce8e 100644 --- a/examples/symba_energy_momentum/escape.in +++ b/examples/symba_energy_momentum/escape.in @@ -1,11 +1,11 @@ 2 -2 1e-07 0.0009 +1 1e-07 0.0009 7e-05 99.9 0.0 0.0 100.00 10.00 0.0 0.4 0.4 0.4 !Ip 0.0 0.0 1000.0 !rot -3 1e-08 0.0004 +2 1e-08 0.0004 3.25e-05 1.0 4.20E-05 0.0 0.00 -6.28 0.0 diff --git a/examples/symba_energy_momentum/param.disruption_headon.in b/examples/symba_energy_momentum/param.disruption_headon.in index 6dbe1f788..4a291535e 100644 --- a/examples/symba_energy_momentum/param.disruption_headon.in +++ b/examples/symba_energy_momentum/param.disruption_headon.in @@ -7,7 +7,7 @@ TP_IN tp.in IN_TYPE ASCII ISTEP_OUT 1 ! output cadence every year BIN_OUT bin.disruption_headon.dat -PARTICLE_FILE particle.disruption_headon.dat +PARTICLE_OUT particle.disruption_headon.dat OUT_TYPE REAL8 ! double precision real output OUT_FORM XV ! osculating element output OUT_STAT REPLACE @@ -28,3 +28,4 @@ TU2S 3.1556925e7 DU2M 1.49598e11 ENERGY yes ROTATION yes +SEED 8 -223172604 -194186007 -2119403444 -114322815 -526658307 1075354356 2043693954 575062362 diff --git a/examples/symba_energy_momentum/param.disruption_off_axis.in b/examples/symba_energy_momentum/param.disruption_off_axis.in index 39303284e..0dfbae80a 100644 --- a/examples/symba_energy_momentum/param.disruption_off_axis.in +++ b/examples/symba_energy_momentum/param.disruption_off_axis.in @@ -7,7 +7,7 @@ TP_IN tp.in IN_TYPE ASCII ISTEP_OUT 1 ! output cadence every year BIN_OUT bin.disruption_off_axis.dat -PARTICLE_FILE particle.disruption_off_axis.dat +PARTICLE_OUT particle.disruption_off_axis.dat OUT_TYPE REAL8 ! double precision real output OUT_FORM XV ! osculating element output OUT_STAT REPLACE @@ -29,3 +29,4 @@ TU2S 3.1556925e7 DU2M 1.49598e11 ENERGY yes ROTATION yes +SEED 8 933097 -220886113 -118730874 233084005 32111237 -823335422 524551114 -61162322 diff --git a/examples/symba_energy_momentum/param.escape.in b/examples/symba_energy_momentum/param.escape.in index 99d572b75..90d118017 100644 --- a/examples/symba_energy_momentum/param.escape.in +++ b/examples/symba_energy_momentum/param.escape.in @@ -7,7 +7,7 @@ TP_IN tp.in IN_TYPE ASCII ISTEP_OUT 1 ! output cadence every year BIN_OUT bin.escape.dat -PARTICLE_FILE particle.escape.dat +PARTICLE_OUT particle.escape.dat OUT_TYPE REAL8 ! double precision real output OUT_FORM XV ! osculating element output OUT_STAT REPLACE @@ -30,3 +30,4 @@ TU2S 3.1556925e7 DU2M 1.49598e11 ENERGY yes ROTATION yes +SEED 8 -1109809 -120983313 -335849874 123308005 -625127 322235652 -3405804 -113111354 diff --git a/examples/symba_energy_momentum/param.sun.in b/examples/symba_energy_momentum/param.sun.in index 5e26e4cd3..a7748b19c 100644 --- a/examples/symba_energy_momentum/param.sun.in +++ b/examples/symba_energy_momentum/param.sun.in @@ -9,7 +9,7 @@ IN_TYPE ASCII ISTEP_OUT 1 ISTEP_DUMP 1 BIN_OUT bin.sun.dat -PARTICLE_FILE particle.sun.dat +PARTICLE_OUT particle.sun.dat OUT_TYPE REAL8 OUT_FORM XV ! osculating element output OUT_STAT REPLACE @@ -30,3 +30,4 @@ TU2S 3.1556925e7 DU2M 1.49598e11 ENERGY yes ROTATION yes +SEED 8 1230834 2346113 123409874 -123121105 -767545 -534058022 343309814 -12535638 diff --git a/examples/symba_energy_momentum/param.supercatastrophic_headon.in b/examples/symba_energy_momentum/param.supercatastrophic_headon.in index 3ba223ad9..e9b60e7da 100644 --- a/examples/symba_energy_momentum/param.supercatastrophic_headon.in +++ b/examples/symba_energy_momentum/param.supercatastrophic_headon.in @@ -7,7 +7,7 @@ TP_IN tp.in IN_TYPE ASCII ISTEP_OUT 1 ! output cadence every year BIN_OUT bin.supercatastrophic_headon.dat -PARTICLE_FILE particle.supercatastrophic_headon.dat +PARTICLE_OUT particle.supercatastrophic_headon.dat OUT_TYPE REAL8 ! double precision real output OUT_FORM XV ! osculating element output OUT_STAT REPLACE @@ -28,3 +28,4 @@ TU2S 3.1556925e7 DU2M 1.49598e11 ENERGY yes ROTATION yes +SEED 8 97 120384098 122231114 -1133345 112137 -239375422 120938114 -66674667 diff --git a/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in b/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in index 49b8b0dd7..0bf836be5 100644 --- a/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in +++ b/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in @@ -7,7 +7,7 @@ TP_IN tp.in IN_TYPE ASCII ISTEP_OUT 1 ! output cadence every year BIN_OUT bin.supercatastrophic_off_axis.dat -PARTICLE_FILE particle.supercatastrophic_off_axis.dat +PARTICLE_OUT particle.supercatastrophic_off_axis.dat OUT_TYPE REAL8 ! double precision real output OUT_FORM XV ! osculating element output OUT_STAT REPLACE @@ -28,3 +28,4 @@ TU2S 3.1556925e7 DU2M 1.49598e11 ENERGY yes ROTATION yes +SEED 8 92823097 -121212113 -736464874 352424135 34555257 -113243092 5640304 -394697 diff --git a/examples/symba_energy_momentum/sun.in b/examples/symba_energy_momentum/sun.in index 7117d93c3..2f3904e5d 100644 --- a/examples/symba_energy_momentum/sun.in +++ b/examples/symba_energy_momentum/sun.in @@ -1,11 +1,11 @@ 2 -2 2e-08 +1 2e-08 3e-04 5e-2 0.0 0.0 0.00 10.00 0.0 0.4 0.4 0.4 !Ip 100.0 100000.0 -2300.0 !rot -3 2e-08 +2 2e-08 3e-06 1.0 0.00E-05 0.0 0.00 6.28 0.0 diff --git a/examples/symba_energy_momentum/supercatastrophic_headon.in b/examples/symba_energy_momentum/supercatastrophic_headon.in index 7b420c9a0..6894837f9 100644 --- a/examples/symba_energy_momentum/supercatastrophic_headon.in +++ b/examples/symba_energy_momentum/supercatastrophic_headon.in @@ -1,11 +1,11 @@ 2 -2 1e-07 0.0009 +1 1e-07 0.0009 7e-06 1.0 -4.20E-05 0.0 0.00 6.28 0.0 0.4 0.4 0.4 !Ip 0.0 0.0 -6.0e4 !rot -3 1e-08 0.0004 +2 1e-08 0.0004 3.25e-06 1.0 4.20E-05 0.0 0.00 -6.28 0.0 diff --git a/examples/symba_energy_momentum/supercatastrophic_off_axis.in b/examples/symba_energy_momentum/supercatastrophic_off_axis.in index a464d037e..230ed071f 100644 --- a/examples/symba_energy_momentum/supercatastrophic_off_axis.in +++ b/examples/symba_energy_momentum/supercatastrophic_off_axis.in @@ -1,11 +1,11 @@ 2 -2 1e-07 0.0009 +1 1e-07 0.0009 7e-06 1.0 -4.20E-05 0.0 0.00 6.28 0.0 0.4 0.4 0.4 !Ip 0.0 0.0 -6.0e4 !rot -3 1e-08 0.0004 +2 1e-08 0.0004 3.25e-06 1.0 4.20E-05 0.0 1.00 -6.28 0.0 diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 81840ba05..0492a05f8 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -24,6 +24,7 @@ def real2float(realstr): """ return float(realstr.replace('d', 'E').replace('D', 'E')) + def read_swiftest_param(param_file_name, param): """ Reads in a Swiftest param.in file and saves it as a dictionary @@ -72,6 +73,8 @@ def read_swiftest_param(param_file_name, param): param['CHK_CLOSE'] = param['CHK_CLOSE'].upper() param['RHILL_PRESENT'] = param['RHILL_PRESENT'].upper() param['FRAGMENTATION'] = param['FRAGMENTATION'].upper() + if param['FRAGMENTATION'] == 'YES' and param['PARTICLE_OUT'] == '': + param['PARTICLE_OUT'] = 'particle.dat' param['ROTATION'] = param['ROTATION'].upper() param['TIDES'] = param['TIDES'].upper() param['ENERGY'] = param['ENERGY'].upper() @@ -82,6 +85,7 @@ def read_swiftest_param(param_file_name, param): print(f"{param_file_name} not found.") return param + def read_swifter_param(param_file_name): """ Reads in a Swifter param.in file and saves it as a dictionary @@ -165,6 +169,7 @@ def read_swifter_param(param_file_name): return param + def read_swift_param(param_file_name, startfile="swift.in"): """ Reads in a Swift param.in file and saves it as a dictionary @@ -251,6 +256,7 @@ def read_swift_param(param_file_name, startfile="swift.in"): return param + def write_swift_param(param, param_file_name): outfile = open(param_file_name, 'w') print(param['T0'], param['TSTOP'], param['DT'], file=outfile) @@ -262,6 +268,7 @@ def write_swift_param(param, param_file_name): outfile.close() return + def write_labeled_param(param, param_file_name): outfile = open(param_file_name, 'w') keylist = ['! VERSION', @@ -300,6 +307,7 @@ def write_labeled_param(param, param_file_name): outfile.close() return + def swifter_stream(f, param): """ Reads in a Swifter bin.dat file and returns a single frame of data as a datastream @@ -318,7 +326,7 @@ def swifter_stream(f, param): plid : int array IDs of massive bodies pvec : float array - (npl,N) - vector of N quantities or each particle (6 of XV/EL + Mass, Radius, etc) + (npl,N) - vector of N quantities or each particle (6 of XV/EL + GMass, Radius, etc) plab : string list Labels for the pvec data ntp : int @@ -376,7 +384,7 @@ def swifter_stream(f, param): tlab.append('omega') tlab.append('capm') plab = tlab.copy() - plab.append('Mass') + plab.append('GMass') plab.append('Radius') pvec = np.vstack([pvec, Mpl, Rpl]) @@ -401,11 +409,11 @@ def make_swiftest_labels(param): tlab.append('omega') tlab.append('capm') plab = tlab.copy() - plab.append('Mass') + plab.append('GMass') plab.append('Radius') if param['RHILL_PRESENT'] == 'YES': plab.append('Rhill') - clab = ['Mass', 'Radius', 'J_2', 'J_4'] + clab = ['GMass', 'Radius', 'J_2', 'J_4'] if param['ROTATION'] == 'YES': clab.append('Ip_x') clab.append('Ip_y') @@ -443,13 +451,13 @@ def swiftest_stream(f, param): cbid : int array ID of central body (always returns 0) cvec : float array - (npl,1) - vector of quantities for the massive body (Mass, Radius, J2, J4, etc) + (npl,1) - vector of quantities for the massive body (GMass, Radius, J2, J4, etc) npl : int Number of massive bodies plid : int array IDs of massive bodies pvec : float array - (npl,N) - vector of N quantities or each particle (6 of XV/EL + Mass, Radius, etc) + (npl,N) - vector of N quantities or each particle (6 of XV/EL + GMass, Radius, etc) plab : string list Labels for the pvec data ntp : int @@ -544,6 +552,7 @@ def swiftest_stream(f, param): npl, plid, pvec.T, plab, \ ntp, tpid, tvec.T, tlab + def swifter2xr(param): """ Converts a Swifter binary data file into an xarray DataSet. @@ -586,6 +595,7 @@ def swifter2xr(param): print(f"Successfully converted {ds.sizes['time']} output frames.") return ds + def swiftest2xr(param): """ Converts a Swiftest binary data file into an xarray DataSet. @@ -604,26 +614,29 @@ def swiftest2xr(param): cb = [] pl = [] tp = [] - with FortranFile(param['BIN_OUT'], 'r') as f: - for t, cbid, cvec, clab, \ - npl, plid, pvec, plab, \ - ntp, tpid, tvec, tlab in swiftest_stream(f, param): - # Prepare frames by adding an extra axis for the time coordinate - cbframe = np.expand_dims(cvec, axis=0) - plframe = np.expand_dims(pvec, axis=0) - tpframe = np.expand_dims(tvec, axis=0) + try: + with FortranFile(param['BIN_OUT'], 'r') as f: + for t, cbid, cvec, clab, \ + npl, plid, pvec, plab, \ + ntp, tpid, tvec, tlab in swiftest_stream(f, param): + # Prepare frames by adding an extra axis for the time coordinate + cbframe = np.expand_dims(cvec, axis=0) + plframe = np.expand_dims(pvec, axis=0) + tpframe = np.expand_dims(tvec, axis=0) + + # Create xarray DataArrays out of each body type + cbxr = xr.DataArray(cbframe, dims=dims, coords={'time': t, 'id': cbid, 'vec': clab}) + plxr = xr.DataArray(plframe, dims=dims, coords={'time': t, 'id': plid, 'vec': plab}) + tpxr = xr.DataArray(tpframe, dims=dims, coords={'time': t, 'id': tpid, 'vec': tlab}) + + cb.append(cbxr) + pl.append(plxr) + tp.append(tpxr) + sys.stdout.write('\r' + f"Reading in time {t[0]:.3e}") + sys.stdout.flush() + except IOError: + print(f"Error encountered reading in {param['BIN_OUT']}") - # Create xarray DataArrays out of each body type - cbxr = xr.DataArray(cbframe, dims=dims, coords={'time': t, 'id': cbid, 'vec': clab}) - plxr = xr.DataArray(plframe, dims=dims, coords={'time': t, 'id': plid, 'vec': plab}) - tpxr = xr.DataArray(tpframe, dims=dims, coords={'time': t, 'id': tpid, 'vec': tlab}) - - cb.append(cbxr) - pl.append(plxr) - tp.append(tpxr) - sys.stdout.write('\r' + f"Reading in time {t[0]:.3e}") - sys.stdout.flush() - cbda = xr.concat(cb, dim='time') plda = xr.concat(pl, dim='time') tpda = xr.concat(tp, dim='time') @@ -634,8 +647,75 @@ def swiftest2xr(param): print('\nCreating Dataset') ds = xr.combine_by_coords([cbds, plds, tpds]) print(f"Successfully converted {ds.sizes['time']} output frames.") + if param['PARTICLE_OUT'] != "": + ds = swiftest_particle_2xr(ds, param) + return ds + +def swiftest_particle_stream(f): + """ + Reads in a Swiftest particle.dat file and returns a single frame of particle data as a datastream + + Parameters + ---------- + f : file object + param : dict + + Yields + ------- + plid : int + ID of massive bodie + origin_type : string + The origin type for the body (Initial conditions, disruption, supercatastrophic, hit and run, etc) + origin_xh : float array + The origin heliocentric position vector + origin_vh : float array + The origin heliocentric velocity vector + """ + while True: # Loop until you read the end of file + try: + # Read multi-line header + plid = f.read_ints() # Try first part of the header + except: + break + origin_rec = f.read_record(np.dtype('a32'), np.dtype((' 0, drop=True) - pl = pl.where(np.invert(np.isnan(pl['Mass'])), drop=True).drop_vars(['J_2', 'J_4']) - tp = frame.where(np.isnan(frame['Mass']), drop=True).drop_vars(['Mass', 'Radius', 'J_2', 'J_4']) + pl = pl.where(np.invert(np.isnan(pl['GMass'])), drop=True).drop_vars(['J_2', 'J_4']) + tp = frame.where(np.isnan(frame['GMass']), drop=True).drop_vars(['GMass', 'Radius', 'J_2', 'J_4']) - GMSun = np.double(cb['Mass']) + GMSun = np.double(cb['GMass']) RSun = np.double(cb['Radius']) J2 = np.double(cb['J_2']) J4 = np.double(cb['J_4']) @@ -680,9 +760,9 @@ def swiftest_xr2infile(ds, param, framenum=-1): for i in pl.id: pli = pl.sel(id=i) if param['RHILL_PRESENT'] == 'YES': - print(i.values, pli['Mass'].values, pli['Rhill'].values, file=plfile) + print(i.values, pli['GMass'].values, pli['Rhill'].values, file=plfile) else: - print(i.values, pli['Mass'].values, file=plfile) + print(i.values, pli['GMass'].values, file=plfile) print(pli['Radius'].values, file=plfile) print(pli['px'].values, pli['py'].values, pli['pz'].values, file=plfile) print(pli['vx'].values, pli['vy'].values, pli['vz'].values, file=plfile) @@ -716,7 +796,7 @@ def swiftest_xr2infile(ds, param, framenum=-1): vx = pl['vx'].values vy = pl['vy'].values vz = pl['vz'].values - mass = pl['Mass'].values + Gmass = pl['GMass'].values radius = pl['Radius'].values plfile.write_record(npl) @@ -727,7 +807,7 @@ def swiftest_xr2infile(ds, param, framenum=-1): plfile.write_record(vx) plfile.write_record(vy) plfile.write_record(vz) - plfile.write_record(mass) + plfile.write_record(Gmass) if param['RHILL_PRESENT'] == 'YES': rhill = pl['Rhill'].values plfile.write_record(rhill) @@ -774,10 +854,10 @@ def swifter_xr2infile(ds, param, framenum=-1): frame = ds.isel(time=framenum) cb = frame.where(frame.id == 0, drop=True) pl = frame.where(frame.id > 0, drop=True) - pl = pl.where(np.invert(np.isnan(pl['Mass'])), drop=True).drop_vars(['J_2', 'J_4']) - tp = frame.where(np.isnan(frame['Mass']), drop=True).drop_vars(['Mass', 'Radius', 'J_2', 'J_4']) + pl = pl.where(np.invert(np.isnan(pl['GMass'])), drop=True).drop_vars(['J_2', 'J_4']) + tp = frame.where(np.isnan(frame['GMass']), drop=True).drop_vars(['GMass', 'Radius', 'J_2', 'J_4']) - GMSun = np.double(cb['Mass']) + GMSun = np.double(cb['GMass']) RSun = np.double(cb['Radius']) param['J2'] = np.double(cb['J_2']) param['J4'] = np.double(cb['J_4']) @@ -786,15 +866,15 @@ def swifter_xr2infile(ds, param, framenum=-1): # Swiftest Central body file plfile = open(param['PL_IN'], 'w') print(pl.id.count().values + 1, file=plfile) - print(cb.id.values[0], cb['Mass'].values[0], file=plfile) + print(cb.id.values[0], cb['GMass'].values[0], file=plfile) print('0.0 0.0 0.0', file=plfile) print('0.0 0.0 0.0', file=plfile) for i in pl.id: pli = pl.sel(id=i) if param['RHILL_PRESENT'] == "YES": - print(i.values, pli['Mass'].values, pli['Rhill'].values, file=plfile) + print(i.values, pli['GMass'].values, pli['Rhill'].values, file=plfile) else: - print(i.values, pli['Mass'].values, file=plfile) + print(i.values, pli['GMass'].values, file=plfile) if param['CHK_CLOSE'] == "YES": print(pli['Radius'].values, file=plfile) print(pli['px'].values, pli['py'].values, pli['pz'].values, file=plfile) @@ -1234,6 +1314,10 @@ def swifter2swiftest(swifter_param, plname="", tpname="", cbname="", conversion_ if swiftest_param['DISCARD_OUT'] == '': swiftest_param['DISCARD_OUT'] = "discard.out" + swiftest_param['PARTICLE_OUT'] = conversion_questions.get('PARTICLE_OUT', '') + if not swiftest_param['PARTICLE_OUT']: + swiftest_param['PARTICLE_OUT'] = input("PARTICLE_OUT: Particle info file name (Only used by SyMBA): []> ") + swiftest_param['! VERSION'] = "Swiftest parameter file converted from Swifter" return swiftest_param diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 05b6896b1..fc5075ab9 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -38,6 +38,7 @@ def __init__(self, codename="Swiftest", param_file=""): 'DU2M': constants.AU2M, 'EXTRA_FORCE': "NO", 'DISCARD_OUT': "discard.out", + 'PARTICLE_OUT' : "", 'BIG_DISCARD': "NO", 'CHK_CLOSE': "YES", 'RHILL_PRESENT': "YES", @@ -52,6 +53,7 @@ def __init__(self, codename="Swiftest", param_file=""): self.read_param(param_file, codename) return + def add(self, plname, date=date.today().isoformat(), idval=None): """ Adds a solar system body to an existing simulation DataSet. @@ -69,6 +71,7 @@ def add(self, plname, date=date.today().isoformat(), idval=None): self.ds = init_cond.solar_system_horizons(plname, idval, self.param, date, self.ds) return + def read_param(self, param_file, codename="Swiftest"): if codename == "Swiftest": self.param = io.read_swiftest_param(param_file, self.param) @@ -84,6 +87,7 @@ def read_param(self, param_file, codename="Swiftest"): self.codename = "Unknown" return + def write_param(self, param_file, param=None): if param is None: param = self.param @@ -97,6 +101,7 @@ def write_param(self, param_file, param=None): print('Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".') return + def convert(self, param_file, newcodename="Swiftest", plname="pl.swiftest.in", tpname="tp.swiftest.in", cbname="cb.swiftest.in", conversion_questions={}): """ Converts simulation input files from one code type to another (Swift, Swifter, or Swiftest). Returns the old parameter configuration. @@ -130,6 +135,7 @@ def convert(self, param_file, newcodename="Swiftest", plname="pl.swiftest.in", t print(f"Conversion from {self.codename} to {newcodename} is not supported.") return oldparam + def bin2xr(self): if self.codename == "Swiftest": self.ds = io.swiftest2xr(self.param) @@ -143,6 +149,7 @@ def bin2xr(self): print('Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".') return + def follow(self, codestyle="Swifter"): if self.ds is None: self.bin2xr() @@ -163,10 +170,13 @@ def follow(self, codestyle="Swifter"): ifol = None nskp = None fol = tool.follow_swift(self.ds, ifol=ifol, nskp=nskp) + else: + fol = None print('follow.out written') return fol + def save(self, param_file, framenum=-1, codename="Swiftest"): if codename == "Swiftest": io.swiftest_xr2infile(self.ds, self.param, framenum) diff --git a/python/swiftest/swiftest/tool.py b/python/swiftest/swiftest/tool.py index 741dc0f1b..a96610bc2 100644 --- a/python/swiftest/swiftest/tool.py +++ b/python/swiftest/swiftest/tool.py @@ -2,15 +2,17 @@ import numpy as np import os import glob -from pyslalib import slalib import xarray as xr """ Functions that recreate the Swift/Swifter tool programs """ -def sla_dranrm(angle): - func = np.vectorize(slalib.sla_dranrm) - return xr.apply_ufunc(func, angle) +def wrap_angle(angle): + while np.any(angle >= 2 * np.pi): + angle[angle >= 2 * np.pi] -= 2 * np.pi + while np.any(angle < 0.0): + angle[angle < 0.0] += 2 * np.pi + return angle def follow_swift(ds, ifol=None, nskp=None): """ @@ -36,11 +38,11 @@ def follow_swift(ds, ifol=None, nskp=None): ifol = int(intxt) print(f"Following particle {ifol}") if ifol < 0: # Negative numbers are planets - fol = ds.where(np.invert(np.isnan(ds['Mass'])), drop=True) + fol = ds.where(np.invert(np.isnan(ds['GMass'])), drop=True) fol = fol.where(np.invert(np.isnan(fol['a'])), drop=True) # Remove times where this body doesn't exist (but this also gets rid of the central body) fol = fol.isel(id = -ifol - 2) # Take 1 off for 0-indexed arrays in Python, and take 1 more off because the central body is gone elif ifol > 0: # Positive numbers are test particles - fol = ds.where(np.isnan(ds['Mass']), drop=True).drop_vars(['Mass', 'Radius']) + fol = ds.where(np.isnan(ds['GMass']), drop=True).drop_vars(['GMass', 'Radius']) fol = fol.where(np.invert(np.isnan(fol['a'])), drop=True) fol = fol.isel(id = ifol - 1) # Take 1 off for 0-indexed arrays in Python @@ -51,7 +53,7 @@ def follow_swift(ds, ifol=None, nskp=None): dr = 180.0 / np.pi fol['obar'] = fol['capom'] + fol['omega'] fol['obar'] = fol['obar'].fillna(0) - fol['obar'] = sla_dranrm(fol['obar']) + fol['obar'] = wrap_angle(fol['obar']) fol['obar'] = fol['obar'] * dr fol['inc'] = fol['inc'] * dr fol['capom'] = fol['capom'] * dr diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index 3d65a235d..9c8044c61 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -17,7 +17,7 @@ module subroutine discard_system(self, param) lpl_check = allocated(self%pl_discards) ltp_check = allocated(self%tp_discards) - associate(system => self, tp => self%tp, pl => self%pl, tp_discards => self%tp_discards, pl_discards => self%tp_discards) + associate(system => self, tp => self%tp, pl => self%pl, tp_discards => self%tp_discards, pl_discards => self%pl_discards) lpl_discards = .false. ltp_discards = .false. if (lpl_check) then @@ -31,8 +31,10 @@ module subroutine discard_system(self, param) end if if (lpl_discards .or. ltp_discards) call system%write_discard(param) + if (lpl_discards .and. param%lenergy) call self%conservation_report(param, lterminal=.true.) if (lpl_check) call pl_discards%setup(0,param) if (ltp_check) call tp_discards%setup(0,param) + end associate return diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index 0b7f2721a..73cef6240 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -25,7 +25,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, real(DP) :: mscale, dscale, vscale, tscale, Lscale, Escale ! Scale factors that reduce quantities to O(~1) in the collisional system real(DP) :: mtot real(DP), dimension(NDIM) :: xcom, vcom - integer(I4B) :: ii + integer(I4B) :: ii, npl_new logical, dimension(:), allocatable :: lexclude real(DP), dimension(NDIM, 2) :: rot, L_orb real(DP), dimension(:,:), allocatable :: x_frag, v_frag, v_r_unit, v_t_unit, v_h_unit @@ -42,11 +42,12 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, integer(I4B), parameter :: NFRAG_MIN = 7 !! The minimum allowable number of fragments (set to 6 because that's how many unknowns are needed in the tangential velocity calculation) real(DP) :: r_max_start, r_max_start_old, r_max, f_spin real(DP), parameter :: Ltol = 10 * epsilon(1.0_DP) - real(DP), parameter :: Etol = 1e-10_DP + real(DP), parameter :: Etol = 1e-9_DP integer(I4B), parameter :: MAXTRY = 3000 integer(I4B), parameter :: TANTRY = 3 logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes - class(swiftest_parameters), allocatable :: tmpparam + class(swiftest_nbody_system), allocatable :: tmpsys + class(swiftest_parameters), allocatable :: tmpparam if (nfrag < NFRAG_MIN) then write(*,*) "symba_frag_pos needs at least ",NFRAG_MIN," fragments, but only ",nfrag," were given." @@ -59,29 +60,24 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, call ieee_set_halting_mode(IEEE_ALL,fpe_quiet_modes) f_spin = 0.05_DP - mscale = 1.0_DP - dscale = 1.0_DP - vscale = 1.0_DP - tscale = 1.0_DP - Lscale = 1.0_DP - Escale = 1.0_DP - - associate(npl => system%pl%nbody, status => system%pl%status) - allocate(lexclude(npl)) - where (status(1:npl) == INACTIVE) ! Safety check in case one of the included bodies has been previously deactivated - lexclude(1:npl) = .true. - elsewhere - lexclude(1:npl) = .false. - end where - end associate allocate(x_frag, source=xb_frag) allocate(v_frag, source=vb_frag) + associate(pl => system%pl, npl => system%pl%nbody) + npl_new = npl + nfrag + allocate(lexclude(npl_new)) + lexclude(1:npl) = pl%status(1:npl) == INACTIVE + lexclude(npl+1:npl_new) = .true. + end associate + call set_scale_factors() call define_coordinate_system() + call construct_temporary_system() + + ! Calculate the initial energy of the system without the collisional family call calculate_system_energy(linclude_fragments=.false.) - + r_max_start = norm2(x(:,2) - x(:,1)) try = 1 lfailure = .false. @@ -92,7 +88,15 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, ke_avg_deficit = 0.0_DP subtry = 1 do + ! Initialize the fragments with 0 velocity and spin so we can divide up the balance between the tangential, radial, and spin components while conserving momentum + xb_frag(:,:) = 0.0_DP + vb_frag(:,:) = 0.0_DP + rot_frag(:,:) = 0.0_DP + v_t_mag(:) = 0.0_DP + v_r_mag(:) = 0.0_DP call set_fragment_position_vectors() + call calculate_system_energy(linclude_fragments=.true.) + ke_frag_budget = -dEtot - Qloss call set_fragment_tan_vel(lfailure) ke_avg_deficit = ke_avg_deficit - ke_radial subtry = subtry + 1 @@ -100,21 +104,23 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, !write(*,*) 'Trying new arrangement' end do ke_avg_deficit = ke_avg_deficit / subtry - if (lfailure) write(*,*) 'Failed to find tangential velocities' + !if (lfailure) write(*,*) 'Failed to find tangential velocities' if (.not.lfailure) then + call calculate_system_energy(linclude_fragments=.true.) + ke_radial = -dEtot - Qloss call set_fragment_radial_velocities(lfailure) - if (lfailure) write(*,*) 'Failed to find radial velocities' + ! if (lfailure) write(*,*) 'Failed to find radial velocities' if (.not.lfailure) then call calculate_system_energy(linclude_fragments=.true.) - !write(*,*) 'Qloss : ',Qloss - !write(*,*) '-dEtot: ',-dEtot - !write(*,*) 'delta : ',abs((dEtot + Qloss)) + ! write(*,*) 'Qloss : ',Qloss + ! write(*,*) '-dEtot: ',-dEtot + ! write(*,*) 'delta : ',abs((dEtot + Qloss)) if ((abs(dEtot + Qloss) > Etol) .or. (dEtot > 0.0_DP)) then - write(*,*) 'Failed due to high energy error: ',dEtot, abs(dEtot + Qloss) / Etol + !write(*,*) 'Failed due to high energy error: ',dEtot, abs(dEtot + Qloss) / Etol lfailure = .true. else if (abs(dLmag) / Lmag_before > Ltol) then - write(*,*) 'Failed due to high angular momentum error: ', dLmag / Lmag_before + !write(*,*) 'Failed due to high angular momentum error: ', dLmag / Lmag_before lfailure = .true. end if end if @@ -124,719 +130,773 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, call restructure_failed_fragments() try = try + 1 end do - write(*, "(' -------------------------------------------------------------------------------------')") - write(*, "(' Final diagnostic')") - write(*, "(' -------------------------------------------------------------------------------------')") - if (lfailure) then - write(*,*) "symba_frag_pos failed after: ",try," tries" - do ii = 1, nfrag - vb_frag(:, ii) = vcom(:) - end do - else - write(*,*) "symba_frag_pos succeeded after: ",try," tries" - write(*, "(' dL_tot should be very small' )") - write(*,fmtlabel) ' dL_tot |', dLmag / Lmag_before - write(*, "(' dE_tot should be negative and equal to Qloss' )") - write(*,fmtlabel) ' dE_tot |', dEtot / abs(Etot_before) - write(*,fmtlabel) ' Qloss |', -Qloss / abs(Etot_before) - write(*,fmtlabel) ' dE - Qloss |', (Etot_after - Etot_before + Qloss) / abs(Etot_before) - end if - write(*, "(' -------------------------------------------------------------------------------------')") - call restore_scale_factors() + call calculate_system_energy(linclude_fragments=.true.) + + ! write(*, "(' -------------------------------------------------------------------------------------')") + ! write(*, "(' Final diagnostic')") + ! write(*, "(' -------------------------------------------------------------------------------------')") + ! if (lfailure) then + ! write(*,*) "symba_frag_pos failed after: ",try," tries" + ! do ii = 1, nfrag + ! vb_frag(:, ii) = vcom(:) + ! end do + ! else + ! write(*,*) "symba_frag_pos succeeded after: ",try," tries" + ! write(*, "(' dL_tot should be very small' )") + ! write(*,fmtlabel) ' dL_tot |', dLmag / Lmag_before + ! write(*, "(' dE_tot should be negative and equal to Qloss' )") + ! write(*,fmtlabel) ' dE_tot |', dEtot / abs(Etot_before) + ! write(*,fmtlabel) ' Qloss |', -Qloss / abs(Etot_before) + ! write(*,fmtlabel) ' dE - Qloss |', (Etot_after - Etot_before + Qloss) / abs(Etot_before) + ! end if + ! write(*, "(' -------------------------------------------------------------------------------------')") + call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily return contains - ! Because of the complexity of this procedure, we have chosen to break it up into a series of nested subroutines. - - subroutine set_scale_factors() - !! author: David A. Minton - !! - !! Scales dimenional quantities to ~O(1) with respect to the collisional system. This scaling makes it easier for the non-linear minimization - !! to converge on a solution - implicit none - integer(I4B) :: i - - ! Find the center of mass of the collisional system - mtot = sum(mass(:)) - xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot - vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot - - ! Set scale factors - Escale = 0.5_DP * (mass(1) * dot_product(v(:,1), v(:,1)) + mass(2) * dot_product(v(:,2), v(:,2))) - dscale = sum(radius(:)) - mscale = mtot - vscale = sqrt(Escale / mscale) - tscale = dscale / vscale - Lscale = mscale * dscale * vscale - - xcom(:) = xcom(:) / dscale - vcom(:) = vcom(:) / vscale - - mtot = mtot / mscale - mass = mass / mscale - radius = radius / dscale - x = x / dscale - v = v / vscale - L_spin = L_spin / Lscale - do i = 1, 2 - rot(:,i) = L_spin(:,i) / (mass(i) * radius(i)**2 * Ip(3, i)) - end do + ! Because of the complexity of this procedure, we have chosen to break it up into a series of nested subroutines. + subroutine set_scale_factors() + !! author: David A. Minton + !! + !! Scales dimenional quantities to ~O(1) with respect to the collisional system. This scaling makes it easier for the non-linear minimization + !! to converge on a solution + implicit none + integer(I4B) :: i + + ! Find the center of mass of the collisional system + mtot = sum(mass(:)) + xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot + vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot + + ! Set scale factors + dscale = sum(radius(:)) + mscale = mtot + vscale = (mass(1) * norm2(v(:,1) - vcom(:)) + mass(2) * norm2(v(:,2) - vcom(:))) / mtot + tscale = dscale / vscale + Lscale = mscale * dscale * vscale + Escale = mscale * vscale**2 + + xcom(:) = xcom(:) / dscale + vcom(:) = vcom(:) / vscale + + mtot = mtot / mscale + mass = mass / mscale + radius = radius / dscale + x = x / dscale + v = v / vscale + L_spin = L_spin / Lscale + do i = 1, 2 + rot(:,i) = L_spin(:,i) / (mass(i) * radius(i)**2 * Ip(3, i)) + end do - m_frag = m_frag / mscale - rad_frag = rad_frag / dscale - Qloss = Qloss / Escale + m_frag = m_frag / mscale + rad_frag = rad_frag / dscale + Qloss = Qloss / Escale - return - end subroutine set_scale_factors - - subroutine restore_scale_factors() - !! author: David A. Minton - !! - !! Restores dimenional quantities back to the system units - implicit none - integer(I4B) :: i - - call ieee_set_halting_mode(IEEE_ALL,.false.) - ! Restore scale factors - xcom(:) = xcom(:) * dscale - vcom(:) = vcom(:) * vscale - - mtot = mtot * mscale - mass = mass * mscale - radius = radius * dscale - x = x * dscale - v = v * vscale - L_spin = L_spin * Lscale - do i = 1, 2 - rot(:,i) = L_spin(:,i) * (mass(i) * radius(i)**2 * Ip(3, i)) - end do + return + end subroutine set_scale_factors - m_frag = m_frag * mscale - rad_frag = rad_frag * dscale - rot_frag = rot_frag / tscale - x_frag = x_frag * dscale - v_frag = v_frag * vscale - Qloss = Qloss * Escale - do i = 1, nfrag - xb_frag(:, i) = x_frag(:, i) + xcom(:) - vb_frag(:, i) = v_frag(:, i) + vcom(:) - end do + subroutine restore_scale_factors() + !! author: David A. Minton + !! + !! Restores dimenional quantities back to the system units + implicit none + integer(I4B) :: i + + call ieee_set_halting_mode(IEEE_ALL,.false.) + ! Restore scale factors + xcom(:) = xcom(:) * dscale + vcom(:) = vcom(:) * vscale + + mtot = mtot * mscale + mass = mass * mscale + radius = radius * dscale + x = x * dscale + v = v * vscale + L_spin = L_spin * Lscale + do i = 1, 2 + rot(:,i) = L_spin(:,i) * (mass(i) * radius(i)**2 * Ip(3, i)) + end do - Etot_before = Etot_before * Escale - pe_before = pe_before * Escale - ke_spin_before = ke_spin_before * Escale - ke_orbit_before = ke_orbit_before * Escale - Ltot_before = Ltot_before * Lscale - Lmag_before = Lmag_before * Lscale - Etot_after = Etot_after * Escale - pe_after = pe_after * Escale - ke_spin_after = ke_spin_after * Escale - ke_orbit_after = ke_orbit_after * Escale - Ltot_after = Ltot_after * Lscale - Lmag_after = Lmag_after * Lscale - - mscale = 1.0_DP - dscale = 1.0_DP - vscale = 1.0_DP - tscale = 1.0_DP - Lscale = 1.0_DP - Escale = 1.0_DP + m_frag = m_frag * mscale + rad_frag = rad_frag * dscale + rot_frag = rot_frag / tscale + x_frag = x_frag * dscale + v_frag = v_frag * vscale + Qloss = Qloss * Escale - return - end subroutine restore_scale_factors - - subroutine define_coordinate_system() - !! author: David A. Minton - !! - !! Defines the collisional coordinate system, including the unit vectors of both the system and individual fragments. - implicit none - integer(I4B) :: i - real(DP), dimension(NDIM) :: x_cross_v, xc, vc, delta_r, delta_v - real(DP) :: r_col_norm, v_col_norm - - allocate(rmag(nfrag)) - allocate(rotmag(nfrag)) - allocate(v_r_mag(nfrag)) - allocate(v_t_mag(nfrag)) - allocate(v_r_unit(NDIM,nfrag)) - allocate(v_t_unit(NDIM,nfrag)) - allocate(v_h_unit(NDIM,nfrag)) - - rmag(:) = 0.0_DP - rotmag(:) = 0.0_DP - v_r_mag(:) = 0.0_DP - v_t_mag(:) = 0.0_DP - v_r_unit(:,:) = 0.0_DP - v_t_unit(:,:) = 0.0_DP - v_h_unit(:,:) = 0.0_DP - - L_orb(:, :) = 0.0_DP - ! Compute orbital angular momentum of pre-impact system - do i = 1, 2 - xc(:) = x(:, i) - xcom(:) - vc(:) = v(:, i) - vcom(:) - x_cross_v(:) = xc(:) .cross. vc(:) - L_orb(:, i) = mass(i) * x_cross_v(:) - end do + do i = 1, nfrag + xb_frag(:, i) = x_frag(:, i) + xcom(:) + vb_frag(:, i) = v_frag(:, i) + vcom(:) + end do - ! Compute orbital angular momentum of pre-impact system. This will be the normal vector to the collision fragment plane - L_frag_tot(:) = L_spin(:, 1) + L_spin(:, 2) + L_orb(:, 1) + L_orb(:, 2) + Etot_before = Etot_before * Escale + pe_before = pe_before * Escale + ke_spin_before = ke_spin_before * Escale + ke_orbit_before = ke_orbit_before * Escale + Ltot_before = Ltot_before * Lscale + Lmag_before = Lmag_before * Lscale + Etot_after = Etot_after * Escale + pe_after = pe_after * Escale + ke_spin_after = ke_spin_after * Escale + ke_orbit_after = ke_orbit_after * Escale + Ltot_after = Ltot_after * Lscale + Lmag_after = Lmag_after * Lscale + + dLmag = norm2(Ltot_after(:) - Ltot_before(:)) + dEtot = Etot_after - Etot_before + + call tmpsys%rescale(tmpparam, mscale**(-1), dscale**(-1), tscale**(-1)) + + mscale = 1.0_DP + dscale = 1.0_DP + vscale = 1.0_DP + tscale = 1.0_DP + Lscale = 1.0_DP + Escale = 1.0_DP + + return + end subroutine restore_scale_factors - delta_v(:) = v(:, 2) - v(:, 1) - v_col_norm = norm2(delta_v(:)) - delta_r(:) = x(:, 2) - x(:, 1) - r_col_norm = norm2(delta_r(:)) - ! We will initialize fragments on a plane defined by the pre-impact system, with the z-axis aligned with the angular momentum vector - ! and the y-axis aligned with the pre-impact distance vector. - y_col_unit(:) = delta_r(:) / r_col_norm - z_col_unit(:) = L_frag_tot(:) / norm2(L_frag_tot) - ! The cross product of the y- by z-axis will give us the x-axis - x_col_unit(:) = y_col_unit(:) .cross. z_col_unit(:) + subroutine define_coordinate_system() + !! author: David A. Minton + !! + !! Defines the collisional coordinate system, including the unit vectors of both the system and individual fragments. + implicit none + integer(I4B) :: i + real(DP), dimension(NDIM) :: x_cross_v, xc, vc, delta_r, delta_v + real(DP) :: r_col_norm, v_col_norm + + allocate(rmag(nfrag)) + allocate(rotmag(nfrag)) + allocate(v_r_mag(nfrag)) + allocate(v_t_mag(nfrag)) + allocate(v_r_unit(NDIM,nfrag)) + allocate(v_t_unit(NDIM,nfrag)) + allocate(v_h_unit(NDIM,nfrag)) + + rmag(:) = 0.0_DP + rotmag(:) = 0.0_DP + v_r_mag(:) = 0.0_DP + v_t_mag(:) = 0.0_DP + v_r_unit(:,:) = 0.0_DP + v_t_unit(:,:) = 0.0_DP + v_h_unit(:,:) = 0.0_DP + + L_orb(:, :) = 0.0_DP + ! Compute orbital angular momentum of pre-impact system + do i = 1, 2 + xc(:) = x(:, i) - xcom(:) + vc(:) = v(:, i) - vcom(:) + x_cross_v(:) = xc(:) .cross. vc(:) + L_orb(:, i) = mass(i) * x_cross_v(:) + end do + + ! Compute orbital angular momentum of pre-impact system. This will be the normal vector to the collision fragment plane + L_frag_tot(:) = L_spin(:, 1) + L_spin(:, 2) + L_orb(:, 1) + L_orb(:, 2) + + delta_v(:) = v(:, 2) - v(:, 1) + v_col_norm = norm2(delta_v(:)) + delta_r(:) = x(:, 2) - x(:, 1) + r_col_norm = norm2(delta_r(:)) + + ! We will initialize fragments on a plane defined by the pre-impact system, with the z-axis aligned with the angular momentum vector + ! and the y-axis aligned with the pre-impact distance vector. + y_col_unit(:) = delta_r(:) / r_col_norm + z_col_unit(:) = L_frag_tot(:) / norm2(L_frag_tot) + ! The cross product of the y- by z-axis will give us the x-axis + x_col_unit(:) = y_col_unit(:) .cross. z_col_unit(:) + + return + end subroutine define_coordinate_system + + + subroutine construct_temporary_system() + !! Author: David A. Minton + !! + !! Constructs a temporary internal system consisting of active bodies and additional fragments. This internal temporary system is used to calculate system energy with and without fragments + !! and optionally including fragments. + implicit none + ! Internals + logical, dimension(:), allocatable :: lexclude_tmp + + associate(pl => system%pl, npl => system%pl%nbody, cb => system%cb) + if (size(lexclude) /= npl + nfrag) then + allocate(lexclude_tmp(npl_new)) + lexclude_tmp(1:npl) = lexclude(1:npl) + call move_alloc(lexclude_tmp, lexclude) + end if + where (pl%status(1:npl) == INACTIVE) ! Safety check in case one of the included bodies has been previously deactivated + lexclude(1:npl) = .true. + elsewhere + lexclude(1:npl) = .false. + end where + lexclude(npl+1:npl_new) = .true. + if (allocated(tmpparam)) deallocate(tmpparam) + allocate(tmpparam, source=param) + call setup_construct_system(tmpsys, param) + call tmpsys%tp%setup(0, param) + deallocate(tmpsys%cb) + allocate(tmpsys%cb, source=cb) + call tmpsys%pl%setup(npl + nfrag, tmpparam) + call tmpsys%pl%fill(pl, .not.lexclude) + call tmpsys%rescale(tmpparam, mscale, dscale, tscale) + + end associate return - end subroutine define_coordinate_system - - subroutine calculate_system_energy(linclude_fragments) - !! Author: David A. Minton - !! - !! Calculates total system energy, including all bodies in the pl list that do not have a corresponding value of the lexclude array that is true - !! and optionally including fragments. - implicit none - ! Arguments - logical, intent(in) :: linclude_fragments - ! Internals - integer(I4B) :: i, npl_new, nplm - logical, dimension(:), allocatable :: ltmp - logical :: lk_plpl - class(swiftest_nbody_system), allocatable :: tmpsys - - ! Because we're making a copy of symba_pl with the excludes/fragments appended, we need to deallocate the - ! big k_plpl array and recreate it when we're done, otherwise we run the risk of blowing up the memory by - ! allocating two of these ginormous arrays simulteouously. This is not particularly efficient, but as this - ! subroutine should be called relatively infrequently, it shouldn't matter too much. - !if (allocated(pl%k_plpl)) deallocate(pl%k_plpl) - - ! Build the internal planet list out of the non-excluded bodies and optionally with fragments appended. This - ! will get passed to the energy calculation subroutine so that energy is computed exactly the same way is it - ! is in the main program. - associate(pl => system%pl, npl => system%pl%nbody, cb => system%cb) - lk_plpl = allocated(pl%k_plpl) - if (lk_plpl) deallocate(pl%k_plpl) - if (linclude_fragments) then ! Temporarily expand the planet list to feed it into symba_energy - lexclude(family(:)) = .true. - npl_new = npl + nfrag - else - npl_new = npl - end if - call setup_construct_system(tmpsys, param) - call tmpsys%tp%setup(0, param) - deallocate(tmpsys%cb) - allocate(tmpsys%cb, source=cb) - if (allocated(tmpparam)) deallocate(tmpparam) - allocate(tmpparam, source=param) - allocate(ltmp(npl_new)) - ltmp(:) = .false. - ltmp(1:npl) = .true. - call tmpsys%pl%setup(npl_new, param) - call tmpsys%pl%fill(pl, ltmp) - call tmpsys%rescale(tmpparam, mscale, dscale, tscale) - - if (linclude_fragments) then ! Append the fragments if they are included - ! Energy calculation requires the fragments to be in the system barcyentric frame + end subroutine construct_temporary_system + + + subroutine add_fragments_to_tmpsys() + !! Author: David A. Minton + !! + !! Adds fragments to the temporary system pl object + implicit none + ! Internals + integer(I4B) :: i + + associate(pl => system%pl, npl => system%pl%nbody) tmpsys%pl%mass(npl+1:npl_new) = m_frag(1:nfrag) tmpsys%pl%Gmass(npl+1:npl_new) = m_frag(1:nfrag) * tmpparam%GU tmpsys%pl%radius(npl+1:npl_new) = rad_frag(1:nfrag) - tmpsys%pl%xb(:,npl+1:npl_new) = xb_frag(:,1:nfrag) - tmpsys%pl%vb(:,npl+1:npl_new) = vb_frag(:,1:nfrag) - tmpsys%pl%status(npl+1:npl_new) = ACTIVE - if (param%lrotation) then + do concurrent (i = 1:nfrag) + tmpsys%pl%xb(:,npl+i) = xb_frag(:,i) + tmpsys%pl%vb(:,npl+i) = vb_frag(:,i) + tmpsys%pl%xh(:,npl+i) = xb_frag(:,i) - tmpsys%cb%xb(:) + tmpsys%pl%vh(:,npl+i) = vb_frag(:,i) - tmpsys%cb%vb(:) + end do + if (tmpparam%lrotation) then tmpsys%pl%Ip(:,npl+1:npl_new) = Ip_frag(:,1:nfrag) tmpsys%pl%rot(:,npl+1:npl_new) = rot_frag(:,1:nfrag) end if - call tmpsys%pl%b2h(tmpsys%cb) - ltmp(1:npl) = lexclude(1:npl) - ltmp(npl+1:npl_new) = .false. - call move_alloc(ltmp, lexclude) - end if + ! Disable the collisional family for subsequent energy calculations and coordinate shifts + lexclude(family(:)) = .true. + lexclude(npl+1:npl_new) = .false. + where(lexclude(:)) + tmpsys%pl%status(:) = INACTIVE + elsewhere + tmpsys%pl%status(:) = ACTIVE + end where + + end associate + + return + end subroutine add_fragments_to_tmpsys + + + subroutine calculate_system_energy(linclude_fragments) + !! Author: David A. Minton + !! + !! Calculates total system energy, including all bodies in the pl list that do not have a corresponding value of the lexclude array that is true + !! and optionally including fragments. + implicit none + ! Arguments + logical, intent(in) :: linclude_fragments + ! Internals + integer(I4B) :: i, nplm + logical, dimension(:), allocatable :: lexclude_tmp + logical :: lk_plpl + + ! Because we're making a copy of symba_pl with the excludes/fragments appended, we need to deallocate the + ! big k_plpl array and recreate it when we're done, otherwise we run the risk of blowing up the memory by + ! allocating two of these ginormous arrays simulteouously. This is not particularly efficient, but as this + ! subroutine should be called relatively infrequently, it shouldn't matter too much. + + ! Build the internal planet list out of the non-excluded bodies and optionally with fragments appended. This + ! will get passed to the energy calculation subroutine so that energy is computed exactly the same way is it + ! is in the main program. This will temporarily expand the planet list in a temporary system object called tmpsys to feed it into symba_energy + associate(pl => system%pl, npl => system%pl%nbody, cb => system%cb) + + where (lexclude(1:npl_new)) + tmpsys%pl%status(1:npl_new) = INACTIVE + elsewhere + tmpsys%pl%status(1:npl_new) = ACTIVE + end where - where (lexclude(1:npl_new)) - tmpsys%pl%status(1:npl_new) = INACTIVE - end where - - select type(plwksp => tmpsys%pl) - class is (symba_pl) - select type(param) - class is (symba_parameters) - plwksp%nplm = count(plwksp%Gmass > param%Gmtiny / mscale) + select type(plwksp => tmpsys%pl) + class is (symba_pl) + select type(param) + class is (symba_parameters) + plwksp%nplm = count(plwksp%Gmass > param%Gmtiny / mscale) + end select end select - end select - call tmpsys%pl%eucl_index() - call tmpsys%get_energy_and_momentum(param) - - ! Restore the big array - deallocate(tmpsys%pl%k_plpl) - select type(pl) - class is (symba_pl) - select type(param) - class is (symba_parameters) - nplm = count(pl%Gmass > param%Gmtiny) + + lk_plpl = allocated(pl%k_plpl) + if (lk_plpl) deallocate(pl%k_plpl) + + call tmpsys%pl%eucl_index() + + call tmpsys%get_energy_and_momentum(param) + + ! Restore the big array + deallocate(tmpsys%pl%k_plpl) + select type(pl) + class is (symba_pl) + select type(param) + class is (symba_parameters) + pl%nplm = count(pl%Gmass > param%Gmtiny) + end select end select - end select - if (lk_plpl) call pl%eucl_index() - - ! Calculate the current fragment energy and momentum balances - if (linclude_fragments) then - Ltot_after(:) = tmpsys%Lorbit(:) + tmpsys%Lspin(:) - Lmag_after = norm2(Ltot_after(:)) - ke_orbit_after = tmpsys%ke_orbit - ke_spin_after = tmpsys%ke_spin - pe_after = tmpsys%pe - Etot_after = ke_orbit_after + ke_spin_after + pe_after - dEtot = Etot_after - Etot_before - dLmag = norm2(Ltot_after(:) - Ltot_before(:)) - else - Ltot_before(:) = tmpsys%Lorbit(:) + tmpsys%Lspin(:) - Lmag_before = norm2(Ltot_before(:)) - ke_orbit_before = tmpsys%ke_orbit - ke_spin_before = tmpsys%ke_spin - pe_before = tmpsys%pe - Etot_before = ke_orbit_before + ke_spin_before + pe_before - end if - end associate - return - end subroutine calculate_system_energy - - subroutine shift_vector_to_origin(m_frag, vec_frag) - !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Adjusts the position or velocity of the fragments as needed to align them with the origin - implicit none - ! Arguments - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:,:), intent(inout) :: vec_frag !! Fragment positions or velocities in the center of mass frame - - ! Internals - real(DP), dimension(NDIM) :: mvec_frag, COM_offset - integer(I4B) :: i - - mvec_frag(:) = 0.0_DP - - do i = 1, nfrag - mvec_frag = mvec_frag(:) + vec_frag(:,i) * m_frag(i) - end do - COM_offset(:) = -mvec_frag(:) / mtot - do i = 1, nfrag - vec_frag(:, i) = vec_frag(:, i) + COM_offset(:) - end do - return - end subroutine shift_vector_to_origin - - subroutine set_fragment_position_vectors() - !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Initializes the orbits of the fragments around the center of mass. The fragments are initially placed on a plane defined by the - !! pre-impact angular momentum. They are distributed on an ellipse surrounding the center of mass. - !! The initial positions do not conserve energy or momentum, so these need to be adjusted later. - - implicit none - real(DP) :: dis, rad - real(DP), dimension(NDIM) :: L_sigma - logical, dimension(:), allocatable :: loverlap - integer(I4B) :: i, j - - allocate(loverlap(nfrag)) - - ! Place the fragments into a region that is big enough that we should usually not have overlapping bodies - ! An overlapping bodies will collide in the next time step, so it's not a major problem if they do (it just slows the run down) - r_max = r_max_start - rad = sum(radius(:)) - - ! We will treat the first two fragments of the list as special cases. They get initialized the maximum distances apart along the original impactor distance vector. - ! This is done because in a regular disruption, the first body is the largest, the second the second largest, and the rest are smaller equal-mass fragments. - - call random_number(x_frag(:,3:nfrag)) - loverlap(:) = .true. - do while (any(loverlap(3:nfrag))) - x_frag(:, 1) = x(:, 1) - xcom(:) - x_frag(:, 2) = x(:, 2) - xcom(:) - r_max = r_max + 0.1_DP * rad - do i = 3, nfrag - if (loverlap(i)) then - call random_number(x_frag(:,i)) - x_frag(:, i) = 2 * (x_frag(:, i) - 0.5_DP) * r_max + if (lk_plpl) call pl%eucl_index() + + ! Calculate the current fragment energy and momentum balances + if (linclude_fragments) then + Ltot_after(:) = tmpsys%Lorbit(:) + tmpsys%Lspin(:) + Lmag_after = norm2(Ltot_after(:)) + ke_orbit_after = tmpsys%ke_orbit + ke_spin_after = tmpsys%ke_spin + pe_after = tmpsys%pe + Etot_after = ke_orbit_after + ke_spin_after + pe_after + dEtot = Etot_after - Etot_before + dLmag = norm2(Ltot_after(:) - Ltot_before(:)) + else + Ltot_before(:) = tmpsys%Lorbit(:) + tmpsys%Lspin(:) + Lmag_before = norm2(Ltot_before(:)) + ke_orbit_before = tmpsys%ke_orbit + ke_spin_before = tmpsys%ke_spin + pe_before = tmpsys%pe + Etot_before = ke_orbit_before + ke_spin_before + pe_before end if + end associate + + return + end subroutine calculate_system_energy + + + subroutine shift_vector_to_origin(m_frag, vec_frag) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Adjusts the position or velocity of the fragments as needed to align them with the origin + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses + real(DP), dimension(:,:), intent(inout) :: vec_frag !! Fragment positions or velocities in the center of mass frame + + ! Internals + real(DP), dimension(NDIM) :: mvec_frag, COM_offset + integer(I4B) :: i + + mvec_frag(:) = 0.0_DP + + do i = 1, nfrag + mvec_frag = mvec_frag(:) + vec_frag(:,i) * m_frag(i) + end do + COM_offset(:) = -mvec_frag(:) / mtot + do i = 1, nfrag + vec_frag(:, i) = vec_frag(:, i) + COM_offset(:) end do - loverlap(:) = .false. - do j = 1, nfrag - do i = j + 1, nfrag - dis = norm2(x_frag(:,j) - x_frag(:,i)) - loverlap(i) = loverlap(i) .or. (dis <= (rad_frag(i) + rad_frag(j))) + + return + end subroutine shift_vector_to_origin + + + subroutine set_fragment_position_vectors() + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Initializes the orbits of the fragments around the center of mass. The fragments are initially placed on a plane defined by the + !! pre-impact angular momentum. They are distributed on an ellipse surrounding the center of mass. + !! The initial positions do not conserve energy or momentum, so these need to be adjusted later. + + implicit none + real(DP) :: dis, rad + real(DP), dimension(NDIM) :: L_sigma + logical, dimension(:), allocatable :: loverlap + integer(I4B) :: i, j + + allocate(loverlap(nfrag)) + + ! Place the fragments into a region that is big enough that we should usually not have overlapping bodies + ! An overlapping bodies will collide in the next time step, so it's not a major problem if they do (it just slows the run down) + r_max = r_max_start + rad = sum(radius(:)) + + ! We will treat the first two fragments of the list as special cases. They get initialized the maximum distances apart along the original impactor distance vector. + ! This is done because in a regular disruption, the first body is the largest, the second the second largest, and the rest are smaller equal-mass fragments. + + call random_number(x_frag(:,3:nfrag)) + loverlap(:) = .true. + do while (any(loverlap(3:nfrag))) + x_frag(:, 1) = x(:, 1) - xcom(:) + x_frag(:, 2) = x(:, 2) - xcom(:) + r_max = r_max + 0.1_DP * rad + do i = 3, nfrag + if (loverlap(i)) then + call random_number(x_frag(:,i)) + x_frag(:, i) = 2 * (x_frag(:, i) - 0.5_DP) * r_max + end if + end do + loverlap(:) = .false. + do j = 1, nfrag + do i = j + 1, nfrag + dis = norm2(x_frag(:,j) - x_frag(:,i)) + loverlap(i) = loverlap(i) .or. (dis <= (rad_frag(i) + rad_frag(j))) + end do end do end do - end do - call shift_vector_to_origin(m_frag, x_frag) - - do i = 1, nfrag - rmag(i) = norm2(x_frag(:, i)) - v_r_unit(:, i) = x_frag(:, i) / rmag(i) - call random_number(L_sigma(:)) ! Randomize the tangential velocity direction. This helps to ensure that the tangential velocity doesn't completely line up with the angular momentum vector, - ! otherwise we can get an ill-conditioned system - v_h_unit(:, i) = z_col_unit(:) + 2e-1_DP * (L_sigma(:) - 0.5_DP) - v_h_unit(:, i) = v_h_unit(:, i) / norm2(v_h_unit(:, i)) - v_t_unit(:, i) = v_h_unit(:, i) .cross. v_r_unit(:, i) - xb_frag(:,i) = x_frag(:,i) + xcom(:) - end do + call shift_vector_to_origin(m_frag, x_frag) + + do i = 1, nfrag + rmag(i) = norm2(x_frag(:, i)) + v_r_unit(:, i) = x_frag(:, i) / rmag(i) + call random_number(L_sigma(:)) ! Randomize the tangential velocity direction. This helps to ensure that the tangential velocity doesn't completely line up with the angular momentum vector, + ! otherwise we can get an ill-conditioned system + v_h_unit(:, i) = z_col_unit(:) + 2e-1_DP * (L_sigma(:) - 0.5_DP) + v_h_unit(:, i) = v_h_unit(:, i) / norm2(v_h_unit(:, i)) + v_t_unit(:, i) = v_h_unit(:, i) .cross. v_r_unit(:, i) + xb_frag(:,i) = x_frag(:,i) + xcom(:) + end do + + call add_fragments_to_tmpsys() + + xcom(:) = 0.0_DP + do i = 1, nfrag + xcom(:) = xcom(:) + m_frag(i) * xb_frag(:,i) + end do + xcom(:) = xcom(:) / mtot - return - end subroutine set_fragment_position_vectors - - subroutine set_fragment_tan_vel(lerr) - !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Adjusts the tangential velocities and spins of a collection of fragments such that they conserve angular momentum without blowing the fragment kinetic energy budget. - !! This procedure works in several stages, with a goal to solve the angular and linear momentum constraints on the fragments, while still leaving a positive balance of - !! our fragment kinetic energy (ke_frag_budget) that we can put into the radial velocity distribution. - !! - !! The first thing we'll try to do is solve for the tangential velocities of the first 6 fragments, using angular and linear momentum as constraints and an initial - !! tangential velocity distribution for the remaining bodies (if there are any) that distributes their angular momentum equally between them. - !! If that doesn't work and we blow our kinetic energy budget, we will attempt to find a tangential velocity distribution that minimizes the kinetic energy while - !! conserving momentum. - !! - !! A failure will trigger a restructuring of the fragments so we will try new values of the radial position distribution. - implicit none - ! Arguments - logical, intent(out) :: lerr - ! Internals - integer(I4B) :: i - real(DP), parameter :: TOL = 1e-4_DP - real(DP), dimension(:), allocatable :: v_t_initial - type(lambda_obj) :: spinfunc - type(lambda_obj_err) :: objective_function - real(DP), dimension(NDIM) :: L_frag_spin, L_remainder, Li, rot_L, rot_ke - - ! Initialize the fragments with 0 velocity and spin so we can divide up the balance between the tangential, radial, and spin components while conserving momentum - lerr = .false. - vb_frag(:,:) = 0.0_DP - rot_frag(:,:) = 0.0_DP - v_t_mag(:) = 0.0_DP - v_r_mag(:) = 0.0_DP - - call calculate_system_energy(linclude_fragments=.true.) - ke_frag_budget = -dEtot - Qloss - !write(*,*) '***************************************************' - !write(*,*) 'Original dis : ',norm2(x(:,2) - x(:,1)) - !write(*,*) 'r_max : ',r_max - !write(*,*) 'f_spin : ',f_spin - !write(*,*) '***************************************************' - !write(*,*) 'Energy balance so far: ' - !write(*,*) 'ke_frag_budget : ',ke_frag_budget - !write(*,*) 'ke_orbit_before: ',ke_orbit_before - !write(*,*) 'ke_orbit_after : ',ke_orbit_after - !write(*,*) 'ke_spin_before : ',ke_spin_before - !write(*,*) 'ke_spin_after : ',ke_spin_after - !write(*,*) 'pe_before : ',pe_before - !write(*,*) 'pe_after : ',pe_after - !write(*,*) 'Qloss : ',Qloss - !write(*,*) '***************************************************' - if (ke_frag_budget < 0.0_DP) then - write(*,*) 'Negative ke_frag_budget: ',ke_frag_budget - r_max_start = r_max_start / 2 - lerr = .true. return - end if + end subroutine set_fragment_position_vectors - allocate(v_t_initial, mold=v_t_mag) - L_frag_spin(:) = 0.0_DP - ke_frag_spin = 0.0_DP - ! Start the first two bodies with the same rotation as the original two impactors, then distribute the remaining angular momentum among the rest - do i = 1, 2 - rot_frag(:, i) = rot(:, i) - L_frag_spin(:) = L_frag_spin(:) + m_frag(i) * rad_frag(i)**2 * Ip_frag(3, i) * rot_frag(:, i) - end do - L_frag_orb(:) = L_frag_tot(:) - L_frag_spin(:) - L_frag_spin(:) = 0.0_DP - do i = 1, nfrag - ! Convert a fraction (f_spin) of either the remaining angular momentum or kinetic energy budget into spin, whichever gives the smaller rotation so as not to blow any budgets - rot_ke(:) = sqrt(2 * f_spin * ke_frag_budget / (nfrag * m_frag(i) * rad_frag(i)**2 * Ip_frag(3, i))) * L_frag_orb(:) / norm2(L_frag_orb(:)) - rot_L(:) = f_spin * L_frag_orb(:) / (nfrag * m_frag(i) * rad_frag(i)**2 * Ip_frag(3, i)) - if (norm2(rot_ke) < norm2(rot_L)) then - rot_frag(:,i) = rot_frag(:, i) + rot_ke(:) - else - rot_frag(:, i) = rot_frag(:, i) + rot_L(:) + subroutine set_fragment_tan_vel(lerr) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Adjusts the tangential velocities and spins of a collection of fragments such that they conserve angular momentum without blowing the fragment kinetic energy budget. + !! This procedure works in several stages, with a goal to solve the angular and linear momentum constraints on the fragments, while still leaving a positive balance of + !! our fragment kinetic energy (ke_frag_budget) that we can put into the radial velocity distribution. + !! + !! The first thing we'll try to do is solve for the tangential velocities of the first 6 fragments, using angular and linear momentum as constraints and an initial + !! tangential velocity distribution for the remaining bodies (if there are any) that distributes their angular momentum equally between them. + !! If that doesn't work and we blow our kinetic energy budget, we will attempt to find a tangential velocity distribution that minimizes the kinetic energy while + !! conserving momentum. + !! + !! A failure will trigger a restructuring of the fragments so we will try new values of the radial position distribution. + implicit none + ! Arguments + logical, intent(out) :: lerr + ! Internals + integer(I4B) :: i + real(DP), parameter :: TOL = 1e-4_DP + real(DP), dimension(:), allocatable :: v_t_initial + real(DP), dimension(nfrag) :: kefrag + type(lambda_obj) :: spinfunc + type(lambda_obj_err) :: objective_function + real(DP), dimension(NDIM) :: L_frag_spin, L_remainder, Li, rot_L, rot_ke + + lerr = .false. + + if (ke_frag_budget < 0.0_DP) then + write(*,*) 'Negative ke_frag_budget: ',ke_frag_budget + r_max_start = r_max_start / 2 + lerr = .true. + return end if - L_frag_spin(:) = L_frag_spin(:) + m_frag(i) * rad_frag(i)**2 * Ip_frag(3, i) * rot_frag(:, i) - ke_frag_spin = ke_frag_spin + m_frag(i) * Ip_frag(3, i) * rad_frag(i)**2 * dot_product(rot_frag(:, i), rot_frag(:, i)) - end do - ke_frag_spin = 0.5_DP * ke_frag_spin - ! Convert a fraction of the pre-impact angular momentum into fragment spin angular momentum - L_frag_orb(:) = L_frag_tot(:) - L_frag_spin(:) - L_remainder(:) = L_frag_orb(:) - ! Next we will solve for the tangential component of the velocities that both conserves linear momentum and uses the remaining angular momentum not used in spin. - ! This will be done using a linear solver that solves for the tangential velocities of the first 6 fragments, constrained by the linear and angular momentum vectors, - ! which is embedded in a non-linear minimizer that will adjust the tangential velocities of the remaining i>6 fragments to minimize kinetic energy for a given momentum solution - ! The initial conditions fed to the minimizer for the fragments will be the remaining angular momentum distributed between the fragments. - do i = 1, nfrag - v_t_initial(i) = norm2(L_remainder(:)) / ((nfrag - i + 1) * m_frag(i) * norm2(x_frag(:,i))) - Li(:) = m_frag(i) * x_frag(:,i) .cross. v_t_initial(i) * v_t_unit(:, i) - L_remainder(:) = L_remainder(:) - Li(:) - end do - ! Find the local kinetic energy minimum for the system that conserves linear and angular momentum - objective_function = lambda_obj(tangential_objective_function, lerr) - v_t_mag(7:nfrag) = util_minimize_bfgs(objective_function, nfrag-6, v_t_initial(7:nfrag), TOL, lerr) - ! Now that the KE-minimized values of the i>6 fragments are found, calculate the momentum-conserving solution for tangential velociteis - v_t_initial(7:nfrag) = v_t_mag(7:nfrag) - v_t_mag(1:nfrag) = solve_fragment_tan_vel(v_t_mag_input=v_t_initial(7:nfrag), lerr=lerr) - - ! Perform one final shift of the radial velocity vectors to align with the center of mass of the collisional system (the origin) - vb_frag(:,1:nfrag) = vmag_to_vb(v_r_mag(1:nfrag), v_r_unit(:,1:nfrag), v_t_mag(1:nfrag), v_t_unit(:,1:nfrag), m_frag(1:nfrag), vcom(:)) - ! Now do a kinetic energy budget check to make sure we are still within the budget. - ke_frag_orbit = 0.0_DP - do i = 1, nfrag - v_frag(:, i) = vb_frag(:, i) - vcom(:) - ke_frag_orbit = ke_frag_orbit + m_frag(i) * dot_product(vb_frag(:, i), vb_frag(:, i)) - end do - ke_frag_orbit = 0.5_DP * ke_frag_orbit - ke_radial = ke_frag_budget - ke_frag_orbit - ke_frag_spin + allocate(v_t_initial, mold=v_t_mag) - ! If we are over the energy budget, flag this as a failure so we can try again - lerr = (ke_radial < 0.0_DP) - !write(*,*) 'Tangential' - !write(*,*) 'ke_frag_budget: ',ke_frag_budget - !write(*,*) 'ke_frag_orbit : ',ke_frag_orbit - !write(*,*) 'ke_frag_spin : ',ke_frag_spin - !write(*,*) 'ke_radial : ',ke_radial + L_frag_spin(:) = 0.0_DP + ke_frag_spin = 0.0_DP + ! Start the first two bodies with the same rotation as the original two impactors, then distribute the remaining angular momentum among the rest + do i = 1, 2 + rot_frag(:, i) = rot(:, i) + L_frag_spin(:) = L_frag_spin(:) + m_frag(i) * rad_frag(i)**2 * Ip_frag(3, i) * rot_frag(:, i) + end do + L_frag_orb(:) = L_frag_tot(:) - L_frag_spin(:) + L_frag_spin(:) = 0.0_DP + do i = 1, nfrag + ! Convert a fraction (f_spin) of either the remaining angular momentum or kinetic energy budget into spin, whichever gives the smaller rotation so as not to blow any budgets + rot_ke(:) = sqrt(2 * f_spin * ke_frag_budget / (nfrag * m_frag(i) * rad_frag(i)**2 * Ip_frag(3, i))) * L_frag_orb(:) / norm2(L_frag_orb(:)) + rot_L(:) = f_spin * L_frag_orb(:) / (nfrag * m_frag(i) * rad_frag(i)**2 * Ip_frag(3, i)) + if (norm2(rot_ke) < norm2(rot_L)) then + rot_frag(:,i) = rot_frag(:, i) + rot_ke(:) + else + rot_frag(:, i) = rot_frag(:, i) + rot_L(:) + end if + L_frag_spin(:) = L_frag_spin(:) + m_frag(i) * rad_frag(i)**2 * Ip_frag(3, i) * rot_frag(:, i) + ke_frag_spin = ke_frag_spin + m_frag(i) * Ip_frag(3, i) * rad_frag(i)**2 * dot_product(rot_frag(:, i), rot_frag(:, i)) + end do + ke_frag_spin = 0.5_DP * ke_frag_spin + ! Convert a fraction of the pre-impact angular momentum into fragment spin angular momentum + L_frag_orb(:) = L_frag_tot(:) - L_frag_spin(:) + L_remainder(:) = L_frag_orb(:) + ! Next we will solve for the tangential component of the velocities that both conserves linear momentum and uses the remaining angular momentum not used in spin. + ! This will be done using a linear solver that solves for the tangential velocities of the first 6 fragments, constrained by the linear and angular momentum vectors, + ! which is embedded in a non-linear minimizer that will adjust the tangential velocities of the remaining i>6 fragments to minimize kinetic energy for a given momentum solution + ! The initial conditions fed to the minimizer for the fragments will be the remaining angular momentum distributed between the fragments. + do i = 1, nfrag + v_t_initial(i) = norm2(L_remainder(:)) / ((nfrag - i + 1) * m_frag(i) * norm2(x_frag(:,i))) + Li(:) = m_frag(i) * x_frag(:,i) .cross. v_t_initial(i) * v_t_unit(:, i) + L_remainder(:) = L_remainder(:) - Li(:) + end do - return + ! Find the local kinetic energy minimum for the system that conserves linear and angular momentum + objective_function = lambda_obj(tangential_objective_function, lerr) + v_t_mag(7:nfrag) = util_minimize_bfgs(objective_function, nfrag-6, v_t_initial(7:nfrag), TOL, lerr) + ! Now that the KE-minimized values of the i>6 fragments are found, calculate the momentum-conserving solution for tangential velociteis + v_t_initial(7:nfrag) = v_t_mag(7:nfrag) + v_t_mag(1:nfrag) = solve_fragment_tan_vel(v_t_mag_input=v_t_initial(7:nfrag), lerr=lerr) + + ! Perform one final shift of the radial velocity vectors to align with the center of mass of the collisional system (the origin) + vb_frag(:,1:nfrag) = vmag_to_vb(v_r_mag(1:nfrag), v_r_unit(:,1:nfrag), v_t_mag(1:nfrag), v_t_unit(:,1:nfrag), m_frag(1:nfrag), vcom(:)) + call add_fragments_to_tmpsys() + + ! Now do a kinetic energy budget check to make sure we are still within the budget. + kefrag = 0.0_DP + do concurrent(i = 1:nfrag) + v_frag(:, i) = vb_frag(:, i) - vcom(:) + kefrag(i) = m_frag(i) * dot_product(vb_frag(:, i), vb_frag(:, i)) + end do + ke_frag_orbit = 0.5_DP * sum(kefrag(:)) + ke_radial = ke_frag_budget - ke_frag_orbit - ke_frag_spin + + ! If we are over the energy budget, flag this as a failure so we can try again + lerr = (ke_radial < 0.0_DP) + ! write(*,*) 'Tangential' + ! write(*,*) 'Failure? ',lerr + ! write(*,*) 'ke_frag_budget: ',ke_frag_budget + ! write(*,*) 'ke_frag_spin : ',ke_frag_spin + ! write(*,*) 'ke_tangential : ',ke_frag_orbit + ! write(*,*) 'ke_remainder : ',ke_radial - end subroutine set_fragment_tan_vel - - function tangential_objective_function(v_t_mag_input, lerr) result(fval) - !! Author: David A. Minton - !! - !! Objective function for evaluating how close our fragment velocities get to minimizing KE error from our required value - implicit none - ! Arguments - real(DP), dimension(:), intent(in) :: v_t_mag_input !! Unknown tangential component of velocity vector set previously by angular momentum constraint - logical, intent(out) :: lerr !! Error flag - ! Result - real(DP) :: fval - ! Internals - integer(I4B) :: i - real(DP), dimension(:,:), allocatable :: v_shift - real(DP), dimension(:), allocatable :: v_t_new - real(DP) :: keo - - lerr = .false. - - allocate(v_shift(NDIM, nfrag)) - allocate(v_t_new(nfrag)) - - v_t_new(:) = solve_fragment_tan_vel(v_t_mag_input=v_t_mag_input(:), lerr=lerr) - v_shift(:,:) = vmag_to_vb(v_r_mag, v_r_unit, v_t_new, v_t_unit, m_frag, vcom) - - keo = 0.0_DP - do i = 1, nfrag - keo = keo + m_frag(i) * dot_product(v_shift(:, i), v_shift(:, i)) - end do - keo = 0.5_DP * keo - fval = keo - lerr = .false. - return - end function tangential_objective_function - - function solve_fragment_tan_vel(lerr, v_t_mag_input) result(v_t_mag_output) - !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Adjusts the positions, velocities, and spins of a collection of fragments such that they conserve angular momentum - implicit none - ! Arguments - logical, intent(out) :: lerr !! Error flag - real(DP), dimension(:), optional, intent(in) :: v_t_mag_input !! Unknown tangential velocities for fragments 7:nfrag - ! Internals - integer(I4B) :: i - ! Result - real(DP), dimension(:), allocatable :: v_t_mag_output - - real(DP), dimension(2 * NDIM, 2 * NDIM) :: A ! LHS of linear equation used to solve for momentum constraint in Gauss elimination code - real(DP), dimension(2 * NDIM) :: b ! RHS of linear equation used to solve for momentum constraint in Gauss elimination code - real(DP), dimension(NDIM) :: L_lin_others, L_orb_others, L, vtmp - - v_frag(:,:) = 0.0_DP - lerr = .false. - - ! We have 6 constraint equations (2 vector constraints in 3 dimensions each) - ! The first 3 are that the linear momentum of the fragments is zero with respect to the collisional barycenter - ! The second 3 are that the sum of the angular momentum of the fragments is conserved from the pre-impact state - L_lin_others(:) = 0.0_DP - L_orb_others(:) = 0.0_DP - do i = 1, nfrag - if (i <= 2 * NDIM) then ! The tangential velocities of the first set of bodies will be the unknowns we will solve for to satisfy the constraints - A(1:3, i) = m_frag(i) * v_t_unit(:, i) - L(:) = v_r_unit(:, i) .cross. v_t_unit(:, i) - A(4:6, i) = m_frag(i) * rmag(i) * L(:) - else if (present(v_t_mag_input)) then - vtmp(:) = v_t_mag_input(i - 6) * v_t_unit(:, i) - L_lin_others(:) = L_lin_others(:) + m_frag(i) * vtmp(:) - L(:) = x_frag(:, i) .cross. vtmp(:) - L_orb_others(:) = L_orb_others(:) + m_frag(i) * L(:) - end if - end do - b(1:3) = -L_lin_others(:) - b(4:6) = L_frag_orb(:) - L_orb_others(:) - allocate(v_t_mag_output(nfrag)) - v_t_mag_output(1:6) = util_solve_linear_system(A, b, 6, lerr) - if (present(v_t_mag_input)) v_t_mag_output(7:nfrag) = v_t_mag_input(:) - - return - end function solve_fragment_tan_vel - - subroutine set_fragment_radial_velocities(lerr) - !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! - !! Adjust the fragment velocities to set the fragment orbital kinetic energy. This will minimize the difference between the fragment kinetic energy and the energy budget - implicit none - ! Arguments - logical, intent(out) :: lerr - ! Internals - real(DP), parameter :: TOL = 1e-10_DP - integer(I4B) :: i, j - real(DP), dimension(:), allocatable :: v_r_initial, v_r_sigma - real(DP), dimension(:,:), allocatable :: v_r - type(lambda_obj) :: objective_function - - ! Set the "target" ke_orbit_after (the value of the orbital kinetic energy that the fragments ought to have) - - allocate(v_r_initial, source=v_r_mag) - ! Initialize radial velocity magnitudes with a random value that is approximately 10% of that found by distributing the kinetic energy equally - allocate(v_r_sigma, source=v_r_mag) - call random_number(v_r_sigma(1:nfrag)) - v_r_sigma(1:nfrag) = sqrt(1.0_DP + 2 * (v_r_sigma(1:nfrag) - 0.5_DP) * 1e-4_DP) - v_r_initial(1:nfrag) = v_r_sigma(1:nfrag) * sqrt(abs(ke_radial) / (2 * m_frag(1:nfrag) * nfrag)) - - ! Initialize the lambda function using a structure constructor that calls the init method - ! Minimize the ke objective function using the BFGS optimizer - objective_function = lambda_obj(radial_objective_function) - v_r_mag = util_minimize_bfgs(objective_function, nfrag, v_r_initial, TOL, lerr) - ! Shift the radial velocity vectors to align with the center of mass of the collisional system (the origin) - vb_frag(:,1:nfrag) = vmag_to_vb(v_r_mag(1:nfrag), v_r_unit(:,1:nfrag), v_t_mag(1:nfrag), v_t_unit(:,1:nfrag), m_frag(1:nfrag), vcom(:)) - do i = 1, nfrag - v_frag(:, i) = vb_frag(:, i) - vcom(:) - end do - ke_frag_orbit = 0.0_DP - do i = 1, nfrag - ke_frag_orbit = ke_frag_orbit + m_frag(i) * dot_product(vb_frag(:, i), vb_frag(:, i)) - end do - ke_frag_orbit = 0.5_DP * ke_frag_orbit - !write(*,*) 'Radial' - !write(*,*) 'Failure? ',lerr - !write(*,*) 'ke_frag_budget: ',ke_frag_budget - !write(*,*) 'ke_frag_orbit : ',ke_frag_orbit - !write(*,*) 'ke_frag_spin : ',ke_frag_spin - !write(*,*) 'ke_remainder : ',ke_frag_budget - (ke_frag_orbit + ke_frag_spin) - lerr = .false. + return + end subroutine set_fragment_tan_vel - return - end subroutine set_fragment_radial_velocities - - function radial_objective_function(v_r_mag_input) result(fval) - !! Author: David A. Minton - !! - !! Objective function for evaluating how close our fragment velocities get to minimizing KE error from our required value - implicit none - ! Arguments - real(DP), dimension(:), intent(in) :: v_r_mag_input !! Unknown radial component of fragment velocity vector - ! Result - real(DP) :: fval !! The objective function result, which is the square of the difference between the calculated fragment kinetic energy and our target - !! Minimizing this brings us closer to our objective - ! Internals - integer(I4B) :: i - real(DP), dimension(:,:), allocatable :: v_shift - - allocate(v_shift, mold=vb_frag) - v_shift(:,:) = vmag_to_vb(v_r_mag_input, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom) - fval = 2 * ke_frag_budget - do i = 1, nfrag - fval = fval - m_frag(i) * (Ip_frag(3, i) * rad_frag(i)**2 * dot_product(rot_frag(:, i), rot_frag(:, i)) + dot_product(v_shift(:, i), v_shift(:, i))) - end do - ! The following ensures that fval = 0 is a local minimum, which is what the BFGS method is searching for - fval = (fval / (2 * ke_radial))**2 - return + function tangential_objective_function(v_t_mag_input, lerr) result(fval) + !! Author: David A. Minton + !! + !! Objective function for evaluating how close our fragment velocities get to minimizing KE error from our required value + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: v_t_mag_input !! Unknown tangential component of velocity vector set previously by angular momentum constraint + logical, intent(out) :: lerr !! Error flag + ! Result + real(DP) :: fval + ! Internals + integer(I4B) :: i + real(DP), dimension(NDIM,nfrag) :: v_shift + real(DP), dimension(nfrag) :: v_t_new, kearr + real(DP) :: keo - end function radial_objective_function - - function vmag_to_vb(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom) result(vb) - !! Author: David A. Minton - !! - !! Converts radial and tangential velocity magnitudes into barycentric velocity - implicit none - ! Arguments - real(DP), dimension(:), intent(in) :: v_r_mag !! Unknown radial component of fragment velocity vector - real(DP), dimension(:), intent(in) :: v_t_mag !! Tangential component of velocity vector set previously by angular momentum constraint - real(DP), dimension(:,:), intent(in) :: v_r_unit, v_t_unit !! Radial and tangential unit vectors for each fragment - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:), intent(in) :: vcom !! Barycentric velocity of collisional system center of mass - ! Result - real(DP), dimension(:,:), allocatable :: vb - ! Internals - integer(I4B) :: i - - allocate(vb, mold=v_r_unit) - ! Make sure the velocity magnitude stays positive - do i = 1, nfrag - vb(:,i) = abs(v_r_mag(i)) * v_r_unit(:, i) - end do - ! In order to keep satisfying the kinetic energy constraint, we must shift the origin of the radial component of the velocities to the center of mass - call shift_vector_to_origin(m_frag, vb) - - do i = 1, nfrag - vb(:, i) = vb(:, i) + v_t_mag(i) * v_t_unit(:, i) + vcom(:) - end do + lerr = .false. - end function vmag_to_vb - - subroutine restructure_failed_fragments() - !! Author: David A. Minton - !! - !! We failed to find a set of positions and velocities that satisfy all the constraints, and so we will alter the fragments and try again. - implicit none - integer(I4B) :: i - real(DP), dimension(:), allocatable :: m_frag_new, rad_frag_new - real(DP), dimension(:,:), allocatable :: xb_frag_new, vb_frag_new, Ip_frag_new, rot_frag_new - real(DP) :: delta_r, delta_r_max - real(DP), parameter :: ke_avg_deficit_target = 0.0_DP - - ! Introduce a bit of noise in the radius determination so we don't just flip flop between similar failed positions - call random_number(delta_r_max) - delta_r_max = sum(radius(:)) * (1.0_DP + 2e-1_DP * (delta_r_max - 0.5_DP)) - if (try > 2) then - ! Linearly interpolate the last two failed solution ke deficits to find a new distance value to try - delta_r = (r_max_start - r_max_start_old) * (ke_avg_deficit_target - ke_avg_deficit_old) / (ke_avg_deficit - ke_avg_deficit_old) - if (abs(delta_r) > delta_r_max) delta_r = sign(delta_r_max, delta_r) - else - delta_r = delta_r_max - end if - r_max_start_old = r_max_start - r_max_start = r_max_start + delta_r ! The larger lever arm can help if the problem is in the angular momentum step - if (f_spin > epsilon(1.0_DP)) then - f_spin = f_spin / 2 - else - f_spin = 0.0_DP - end if - end subroutine restructure_failed_fragments + v_t_new(:) = solve_fragment_tan_vel(v_t_mag_input=v_t_mag_input(:), lerr=lerr) + v_shift(:,:) = vmag_to_vb(v_r_mag, v_r_unit, v_t_new, v_t_unit, m_frag, vcom) + + kearr = 0.0_DP + do concurrent(i = 1:nfrag) + kearr(i) = m_frag(i) * dot_product(v_shift(:, i), v_shift(:, i)) + end do + keo = 0.5_DP * sum(kearr(:)) + fval = keo + lerr = .false. + + return + end function tangential_objective_function + function solve_fragment_tan_vel(lerr, v_t_mag_input) result(v_t_mag_output) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Adjusts the positions, velocities, and spins of a collection of fragments such that they conserve angular momentum + implicit none + ! Arguments + logical, intent(out) :: lerr !! Error flag + real(DP), dimension(:), optional, intent(in) :: v_t_mag_input !! Unknown tangential velocities for fragments 7:nfrag + ! Internals + integer(I4B) :: i + ! Result + real(DP), dimension(:), allocatable :: v_t_mag_output + + real(DP), dimension(2 * NDIM, 2 * NDIM) :: A ! LHS of linear equation used to solve for momentum constraint in Gauss elimination code + real(DP), dimension(2 * NDIM) :: b ! RHS of linear equation used to solve for momentum constraint in Gauss elimination code + real(DP), dimension(NDIM) :: L_lin_others, L_orb_others, L, vtmp + + v_frag(:,:) = 0.0_DP + lerr = .false. + + ! We have 6 constraint equations (2 vector constraints in 3 dimensions each) + ! The first 3 are that the linear momentum of the fragments is zero with respect to the collisional barycenter + ! The second 3 are that the sum of the angular momentum of the fragments is conserved from the pre-impact state + L_lin_others(:) = 0.0_DP + L_orb_others(:) = 0.0_DP + do i = 1, nfrag + if (i <= 2 * NDIM) then ! The tangential velocities of the first set of bodies will be the unknowns we will solve for to satisfy the constraints + A(1:3, i) = m_frag(i) * v_t_unit(:, i) + L(:) = v_r_unit(:, i) .cross. v_t_unit(:, i) + A(4:6, i) = m_frag(i) * rmag(i) * L(:) + else if (present(v_t_mag_input)) then + vtmp(:) = v_t_mag_input(i - 6) * v_t_unit(:, i) + L_lin_others(:) = L_lin_others(:) + m_frag(i) * vtmp(:) + L(:) = x_frag(:, i) .cross. vtmp(:) + L_orb_others(:) = L_orb_others(:) + m_frag(i) * L(:) + end if + end do + b(1:3) = -L_lin_others(:) + b(4:6) = L_frag_orb(:) - L_orb_others(:) + allocate(v_t_mag_output(nfrag)) + v_t_mag_output(1:6) = util_solve_linear_system(A, b, 6, lerr) + if (present(v_t_mag_input)) v_t_mag_output(7:nfrag) = v_t_mag_input(:) + + return + end function solve_fragment_tan_vel + + + subroutine set_fragment_radial_velocities(lerr) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! + !! Adjust the fragment velocities to set the fragment orbital kinetic energy. This will minimize the difference between the fragment kinetic energy and the energy budget + implicit none + ! Arguments + logical, intent(out) :: lerr + ! Internals + real(DP), parameter :: TOL = 1e-10_DP + integer(I4B) :: i, j + real(DP), dimension(:), allocatable :: v_r_initial, v_r_sigma + real(DP), dimension(:,:), allocatable :: v_r + real(DP), dimension(nfrag) :: kearr, kespinarr + type(lambda_obj) :: objective_function + + ! Set the "target" ke_orbit_after (the value of the orbital kinetic energy that the fragments ought to have) + + allocate(v_r_initial, source=v_r_mag) + ! Initialize radial velocity magnitudes with a random value that is approximately 10% of that found by distributing the kinetic energy equally + allocate(v_r_sigma, source=v_r_mag) + call random_number(v_r_sigma(1:nfrag)) + v_r_sigma(1:nfrag) = sqrt(1.0_DP + 2 * (v_r_sigma(1:nfrag) - 0.5_DP) * 1e-4_DP) + v_r_initial(1:nfrag) = v_r_sigma(1:nfrag) * sqrt(abs(2 * ke_radial) / (m_frag(1:nfrag) * nfrag)) + + ! Initialize the lambda function using a structure constructor that calls the init method + ! Minimize the ke objective function using the BFGS optimizer + objective_function = lambda_obj(radial_objective_function) + v_r_mag = util_minimize_bfgs(objective_function, nfrag, v_r_initial, TOL, lerr) + ! Shift the radial velocity vectors to align with the center of mass of the collisional system (the origin) + vb_frag(:,1:nfrag) = vmag_to_vb(v_r_mag(1:nfrag), v_r_unit(:,1:nfrag), v_t_mag(1:nfrag), v_t_unit(:,1:nfrag), m_frag(1:nfrag), vcom(:)) + do i = 1, nfrag + v_frag(:, i) = vb_frag(:, i) - vcom(:) + end do + call add_fragments_to_tmpsys() + + do concurrent(i = 1:nfrag) + kearr(i) = m_frag(i) * dot_product(vb_frag(:, i), vb_frag(:, i)) + kespinarr(i) = m_frag(i) * Ip_frag(3, i) * rad_frag(i)**2 * dot_product(rot_frag(:,i), rot_frag(:,i)) + end do + ke_frag_orbit = 0.5_DP * sum(kearr(:)) + ke_frag_spin = 0.5_DP * sum(kespinarr(:)) + ! write(*,*) 'Radial' + ! write(*,*) 'Failure? ',lerr + ! write(*,*) 'ke_frag_budget: ',ke_frag_budget + ! write(*,*) 'ke_frag_spin : ',ke_frag_spin + ! write(*,*) 'ke_frag_orbit : ',ke_frag_orbit + ! write(*,*) 'ke_remainder : ',ke_frag_budget - (ke_frag_orbit + ke_frag_spin) + lerr = .false. + + return + end subroutine set_fragment_radial_velocities + + + function radial_objective_function(v_r_mag_input) result(fval) + !! Author: David A. Minton + !! + !! Objective function for evaluating how close our fragment velocities get to minimizing KE error from our required value + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: v_r_mag_input !! Unknown radial component of fragment velocity vector + ! Result + real(DP) :: fval !! The objective function result, which is the square of the difference between the calculated fragment kinetic energy and our target + !! Minimizing this brings us closer to our objective + ! Internals + integer(I4B) :: i + real(DP), dimension(:,:), allocatable :: v_shift + real(DP), dimension(nfrag) :: kearr + real(DP) :: keo + + allocate(v_shift, mold=vb_frag) + v_shift(:,:) = vmag_to_vb(v_r_mag_input, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom) + do concurrent(i = 1:nfrag) + kearr(i) = m_frag(i) * (Ip_frag(3, i) * rad_frag(i)**2 * dot_product(rot_frag(:, i), rot_frag(:, i)) + dot_product(v_shift(:, i), v_shift(:, i))) + end do + keo = 2 * ke_frag_budget - sum(kearr(:)) + ! The following ensures that fval = 0 is a local minimum, which is what the BFGS method is searching for + fval = (keo / (2 * ke_radial))**2 + + return + end function radial_objective_function + + + function vmag_to_vb(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom) result(vb) + !! Author: David A. Minton + !! + !! Converts radial and tangential velocity magnitudes into barycentric velocity + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: v_r_mag !! Unknown radial component of fragment velocity vector + real(DP), dimension(:), intent(in) :: v_t_mag !! Tangential component of velocity vector set previously by angular momentum constraint + real(DP), dimension(:,:), intent(in) :: v_r_unit, v_t_unit !! Radial and tangential unit vectors for each fragment + real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses + real(DP), dimension(:), intent(in) :: vcom !! Barycentric velocity of collisional system center of mass + ! Result + real(DP), dimension(:,:), allocatable :: vb + ! Internals + integer(I4B) :: i + + allocate(vb, mold=v_r_unit) + ! Make sure the velocity magnitude stays positive + do i = 1, nfrag + vb(:,i) = abs(v_r_mag(i)) * v_r_unit(:, i) + end do + ! In order to keep satisfying the kinetic energy constraint, we must shift the origin of the radial component of the velocities to the center of mass + call shift_vector_to_origin(m_frag, vb) + + do i = 1, nfrag + vb(:, i) = vb(:, i) + v_t_mag(i) * v_t_unit(:, i) + vcom(:) + end do + + return + end function vmag_to_vb + + + subroutine restructure_failed_fragments() + !! Author: David A. Minton + !! + !! We failed to find a set of positions and velocities that satisfy all the constraints, and so we will alter the fragments and try again. + implicit none + integer(I4B) :: i + real(DP), dimension(:), allocatable :: m_frag_new, rad_frag_new + real(DP), dimension(:,:), allocatable :: xb_frag_new, vb_frag_new, Ip_frag_new, rot_frag_new + real(DP) :: delta_r, delta_r_max + real(DP), parameter :: ke_avg_deficit_target = 0.0_DP + + ! Introduce a bit of noise in the radius determination so we don't just flip flop between similar failed positions + call random_number(delta_r_max) + delta_r_max = sum(radius(:)) * (1.0_DP + 2e-1_DP * (delta_r_max - 0.5_DP)) + if (try > 2) then + ! Linearly interpolate the last two failed solution ke deficits to find a new distance value to try + delta_r = (r_max_start - r_max_start_old) * (ke_avg_deficit_target - ke_avg_deficit_old) / (ke_avg_deficit - ke_avg_deficit_old) + if (abs(delta_r) > delta_r_max) delta_r = sign(delta_r_max, delta_r) + else + delta_r = delta_r_max + end if + r_max_start_old = r_max_start + r_max_start = r_max_start + delta_r ! The larger lever arm can help if the problem is in the angular momentum step + if (f_spin > epsilon(1.0_DP)) then + f_spin = f_spin / 2 + else + f_spin = 0.0_DP + end if + + return + end subroutine restructure_failed_fragments end subroutine fragmentation_initialize diff --git a/src/io/io.f90 b/src/io/io.f90 index ebfe2b1ef..03fdc2e17 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -31,14 +31,15 @@ module subroutine io_conservation_report(self, param, lterminal) Euntracked => self%Euntracked, Eorbit_orig => param%Eorbit_orig, Mtot_orig => param%Mtot_orig, & Ltot_orig => param%Ltot_orig(:), Lmag_orig => param%Lmag_orig, Lorbit_orig => param%Lorbit_orig(:), Lspin_orig => param%Lspin_orig(:), & lfirst => param%lfirstenergy) - if (lfirst) then - if (param%out_stat == "OLD") then - open(unit = EGYIU, file = ENERGY_FILE, form = "formatted", status = "old", action = "write", position = "append") - else - open(unit = EGYIU, file = ENERGY_FILE, form = "formatted", status = "replace", action = "write") + if (param%energy_out /= "") then + if (lfirst .and. (param%out_stat /= "OLD")) then + open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "replace", action = "write") + else + open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "old", action = "write", position = "append") write(EGYIU,EGYHEADER) end if end if + call pl%h2b(cb) call system%get_energy_and_momentum(param) ke_orbit_now = system%ke_orbit ke_spin_now = system%ke_spin @@ -58,8 +59,10 @@ module subroutine io_conservation_report(self, param, lterminal) lfirst = .false. end if - write(EGYIU,EGYFMT) param%t, Eorbit_now, Ecollisions, Ltot_now, Mtot_now - flush(EGYIU) + if (param%energy_out /= "") then + write(EGYIU,EGYFMT) param%t, Eorbit_now, Ecollisions, Ltot_now, Mtot_now + close(EGYIU) + end if if (.not.lfirst .and. lterminal) then Lmag_now = norm2(Ltot_now) Lerror = norm2(Ltot_now - Ltot_orig) / Lmag_orig @@ -421,6 +424,8 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) param%enc_out = param_value case ("DISCARD_OUT") param%discard_out = param_value + case ("ENERGY_OUT") + param%energy_out = param_value case ("EXTRA_FORCE") call io_toupper(param_value) if (param_value == "YES" .or. param_value == 'T') param%lextra_force = .true. @@ -493,7 +498,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) read(param_value, *) param%Ecollisions case("EUNTRACKED") read(param_value, *) param%Euntracked - case ("NPLMAX", "NTPMAX", "GMTINY", "PARTICLE_FILE", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters + case ("NPLMAX", "NTPMAX", "GMTINY", "PARTICLE_OUT", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters case default write(iomsg,*) "Unknown parameter -> ",param_name iostat = -1 @@ -535,6 +540,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) iostat = -1 return end if + param%lrestart = (param%out_stat == "APPEND") if (param%outfile /= "") then if ((param%out_type /= REAL4_TYPE) .and. (param%out_type /= REAL8_TYPE) .and. & (param%out_type /= SWIFTER_REAL4_TYPE) .and. (param%out_type /= SWIFTER_REAL8_TYPE)) then @@ -552,6 +558,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) iostat = -1 return end if + end if if (param%qmin > 0.0_DP) then if ((param%qmin_coord /= "HELIO") .and. (param%qmin_coord /= "BARY")) then diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index 4f7255237..6ffb7ba1b 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -163,7 +163,7 @@ module subroutine rmvs_util_append_pl(self, source, lsource_mask) implicit none class(rmvs_pl), intent(inout) :: self !! RMVS massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine rmvs_util_append_pl module subroutine rmvs_util_append_tp(self, source, lsource_mask) @@ -171,7 +171,7 @@ module subroutine rmvs_util_append_tp(self, source, lsource_mask) implicit none class(rmvs_tp), intent(inout) :: self !! RMVS test particle object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine rmvs_util_append_tp module subroutine rmvs_util_fill_pl(self, inserts, lfill_list) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 4d0e98704..fa5ec8b97 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -45,7 +45,7 @@ module swiftest_classes real(QP) :: DU2M = -1.0_QP !! Converts distance unit to centimeters real(DP) :: GU = -1.0_DP !! Universal gravitational constant in the system units real(DP) :: inv_c2 = -1.0_DP !! Inverse speed of light squared in the system units - character(STRMAX) :: ennergy_out = "" !! Name of output energy and momentum report file + character(STRMAX) :: energy_out = "" !! Name of output energy and momentum report file ! Logical flags to turn on or off various features of the code logical :: lrhill_present = .false. !! Hill radii are given as an input rather than calculated by the code (can be used to inflate close encounter regions manually) @@ -64,13 +64,14 @@ module swiftest_classes real(DP), dimension(NDIM) :: Ltot_orig = 0.0_DP !! Initial total angular momentum vector real(DP), dimension(NDIM) :: Lorbit_orig = 0.0_DP !! Initial orbital angular momentum real(DP), dimension(NDIM) :: Lspin_orig = 0.0_DP !! Initial spin angular momentum vector - real(DP), dimension(NDIM) :: Ltot = 0.0_DP !! System angular momentum vector - real(DP), dimension(NDIM) :: Lescape = 0.0_DP !! Angular momentum of bodies that escaped the system (used for bookeeping) - real(DP) :: Mescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping) - real(DP) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions - real(DP) :: Euntracked = 0.0_DP !! Energy gained from system due to escaped bodies + real(DP), dimension(NDIM) :: Ltot = 0.0_DP !! System angular momentum vector + real(DP), dimension(NDIM) :: Lescape = 0.0_DP !! Angular momentum of bodies that escaped the system (used for bookeeping) + real(DP) :: Mescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping) + real(DP) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions + real(DP) :: Euntracked = 0.0_DP !! Energy gained from system due to escaped bodies logical :: lfirstenergy = .true. !! This is the first time computing energe logical :: lfirstkick = .true. !! Initiate the first kick in a symplectic step + logical :: lrestart = .false. !! Indicates whether or not this is a restarted run ! Future features not implemented or in development logical :: lgr = .false. !! Turn on GR @@ -832,39 +833,44 @@ end subroutine user_kick_getacch_body end interface interface util_append - module subroutine util_append_arr_char_string(arr, source, lsource_mask) + module subroutine util_append_arr_char_string(arr, source, nold, nsrc, lsource_mask) implicit none character(len=STRMAX), dimension(:), allocatable, intent(inout) :: arr !! Destination array character(len=STRMAX), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_arr_char_string - module subroutine util_append_arr_DP(arr, source, lsource_mask) + module subroutine util_append_arr_DP(arr, source, nold, nsrc, lsource_mask) implicit none real(DP), dimension(:), allocatable, intent(inout) :: arr !! Destination array real(DP), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_arr_DP - module subroutine util_append_arr_DPvec(arr, source, lsource_mask) + module subroutine util_append_arr_DPvec(arr, source, nold, nsrc, lsource_mask) implicit none real(DP), dimension(:,:), allocatable, intent(inout) :: arr !! Destination array real(DP), dimension(:,:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_arr_DPvec - module subroutine util_append_arr_I4B(arr, source, lsource_mask) + module subroutine util_append_arr_I4B(arr, source, nold, nsrc, lsource_mask) implicit none integer(I4B), dimension(:), allocatable, intent(inout) :: arr !! Destination array integer(I4B), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_arr_I4B - module subroutine util_append_arr_logical(arr, source, lsource_mask) + module subroutine util_append_arr_logical(arr, source, nold, nsrc, lsource_mask) implicit none logical, dimension(:), allocatable, intent(inout) :: arr !! Destination array logical, dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_arr_logical end interface @@ -872,22 +878,22 @@ end subroutine util_append_arr_logical module subroutine util_append_body(self, source, lsource_mask) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest body object - class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + class(swiftest_body), intent(in) :: source !! Source object to append + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_body module subroutine util_append_pl(self, source, lsource_mask) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_pl module subroutine util_append_tp(self, source, lsource_mask) implicit none class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_tp module subroutine util_coord_b2h_pl(self, cb) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index dfe1c0326..4628202f8 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -18,7 +18,7 @@ module symba_classes integer(I4B), parameter :: PARTICLEUNIT = 44 !! File unit number for the binary particle info output file type, extends(swiftest_parameters) :: symba_parameters - character(STRMAX) :: particle_file = PARTICLE_OUTFILE !! Name of output particle information file + character(STRMAX) :: particle_out = PARTICLE_OUTFILE !! Name of output particle information file real(DP) :: GMTINY = -1.0_DP !! Smallest mass that is fully gravitating integer(I4B), dimension(:), allocatable :: seed !! Random seeds logical :: lfragmentation = .false. !! Do fragmentation modeling instead of simple merger. @@ -27,33 +27,17 @@ module symba_classes procedure :: writer => symba_io_param_writer end type symba_parameters - !******************************************************************************************************************************** - ! symba_cb class definitions and method interfaces - !******************************************************************************************************************************* - !> SyMBA central body particle class - type, extends(helio_cb) :: symba_cb - real(DP) :: M0 = 0.0_DP !! Initial mass of the central body - real(DP) :: dM = 0.0_DP !! Change in mass of the central body - real(DP) :: R0 = 0.0_DP !! Initial radius of the central body - real(DP) :: dR = 0.0_DP !! Change in the radius of the central body - contains - end type symba_cb - !******************************************************************************************************************************** ! symba_particle_info class definitions and method interfaces !******************************************************************************************************************************* !> Class definition for the particle origin information object. This object is used to track time, location, and collisional regime !> of fragments produced in collisional events. - type, extends(swiftest_base) :: symba_particle_info + type :: symba_particle_info + sequence character(len=32) :: origin_type !! String containing a description of the origin of the particle (e.g. Initial Conditions, Supercatastrophic, Disruption, etc.) real(DP) :: origin_time !! The time of the particle's formation real(DP), dimension(NDIM) :: origin_xh !! The heliocentric distance vector at the time of the particle's formation real(DP), dimension(NDIM) :: origin_vh !! The heliocentric velocity vector at the time of the particle's formation - contains - procedure :: dump => symba_io_dump_particle_info !! I/O routine for dumping particle info to file - procedure :: initialize => symba_io_initialize_particle_info !! I/O routine for reading in particle info data - procedure :: read_frame => symba_io_read_frame_info !! I/O routine for reading in a single frame of particle info - procedure :: write_frame => symba_io_write_frame_info !! I/O routine for writing out a single frame of particle info end type symba_particle_info !******************************************************************************************************************************** @@ -66,6 +50,19 @@ module symba_classes integer(I4B), dimension(:), allocatable :: child !! Index of children particles end type symba_kinship + !******************************************************************************************************************************** + ! symba_cb class definitions and method interfaces + !******************************************************************************************************************************* + !> SyMBA central body particle class + type, extends(helio_cb) :: symba_cb + real(DP) :: M0 = 0.0_DP !! Initial mass of the central body + real(DP) :: dM = 0.0_DP !! Change in mass of the central body + real(DP) :: R0 = 0.0_DP !! Initial radius of the central body + real(DP) :: dR = 0.0_DP !! Change in the radius of the central body + type(symba_particle_info) :: info + contains + end type symba_cb + !******************************************************************************************************************************** ! symba_pl class definitions and method interfaces !******************************************************************************************************************************* @@ -118,6 +115,7 @@ module symba_classes integer(I4B), dimension(:), allocatable :: nplenc !! number of encounters with planets this time step integer(I4B), dimension(:), allocatable :: levelg !! level at which this particle should be moved integer(I4B), dimension(:), allocatable :: levelm !! deepest encounter level achieved this time step + type(symba_particle_info), dimension(:), allocatable :: info contains procedure :: drift => symba_drift_tp !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level procedure :: encounter_check => symba_encounter_check_tp !! Checks if any test particles are undergoing a close encounter with a massive body @@ -327,12 +325,13 @@ module subroutine symba_io_write_discard(self, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine symba_io_write_discard - module subroutine symba_io_dump_particle_info(self, param, msg) - use swiftest_classes, only : swiftest_parameters + module subroutine symba_io_dump_particle_info(system, param, lincludecb, tpidx, plidx) implicit none - class(symba_particle_info), intent(inout) :: self !! Swiftest base object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - character(*), optional, intent(in) :: msg !! Message to display with dump operation + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA extensions + logical, optional, intent(in) :: lincludecb !! Set to true to include the central body (default is false) + integer(I4B), dimension(:), optional, intent(in) :: tpidx !! Array of test particle indices to append to the particle file + integer(I4B), dimension(:), optional, intent(in) :: plidx !! Array of massive body indices to append to the particle file end subroutine symba_io_dump_particle_info module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) @@ -356,23 +355,12 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms integer, intent(out) :: iostat !! IO status code character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 end subroutine symba_io_param_writer - - module subroutine symba_io_initialize_particle_info(self, param) - use swiftest_classes, only : swiftest_parameters - implicit none - class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine symba_io_initialize_particle_info - module subroutine symba_io_read_frame_info(self, iu, param, form, ierr) - use swiftest_classes, only : swiftest_parameters + module subroutine symba_io_read_particle(system, param) implicit none - class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - character(*), intent(in) :: form !! Input format code ("XV" or "EL") - integer(I4B), intent(out) :: ierr !! Error code - end subroutine symba_io_read_frame_info + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system file + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA extensions + end subroutine symba_io_read_particle module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters @@ -402,14 +390,12 @@ module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn) integer(I4B), intent(in) :: irec !! Current recursion level integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration end subroutine symba_kick_pltpenc - - module subroutine symba_io_write_frame_info(self, iu, param) - use swiftest_classes, only : swiftest_parameters + + module subroutine symba_setup_initialize_particle_info(system, param) implicit none - class(symba_particle_info), intent(in) :: self !! SyMBA particle info object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine symba_io_write_frame_info + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA extensions + end subroutine symba_setup_initialize_particle_info module subroutine symba_setup_initialize_system(self, param) use swiftest_classes, only : swiftest_parameters @@ -488,18 +474,20 @@ end subroutine symba_step_reset_system end interface interface util_append - module subroutine symba_util_append_arr_info(arr, source, lsource_mask) + module subroutine symba_util_append_arr_info(arr, source, nold, nsrc, lsource_mask) implicit none type(symba_particle_info), dimension(:), allocatable, intent(inout) :: arr !! Destination array type(symba_particle_info), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine symba_util_append_arr_info - module subroutine symba_util_append_arr_kin(arr, source, lsource_mask) + module subroutine symba_util_append_arr_kin(arr, source, nold, nsrc, lsource_mask) implicit none type(symba_kinship), dimension(:), allocatable, intent(inout) :: arr !! Destination array type(symba_kinship), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine symba_util_append_arr_kin end interface @@ -509,7 +497,7 @@ module subroutine symba_util_append_merger(self, source, lsource_mask) implicit none class(symba_merger), intent(inout) :: self !! SyMBA massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine symba_util_append_merger module subroutine symba_util_append_pl(self, source, lsource_mask) @@ -517,7 +505,7 @@ module subroutine symba_util_append_pl(self, source, lsource_mask) implicit none class(symba_pl), intent(inout) :: self !! SyMBA massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine symba_util_append_pl module subroutine symba_util_append_tp(self, source, lsource_mask) @@ -525,7 +513,7 @@ module subroutine symba_util_append_tp(self, source, lsource_mask) implicit none class(symba_tp), intent(inout) :: self !! SyMBA test particle object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine symba_util_append_tp end interface diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90 index a79f52bca..6d5b26394 100644 --- a/src/modules/whm_classes.f90 +++ b/src/modules/whm_classes.f90 @@ -232,7 +232,7 @@ module subroutine whm_util_append_pl(self, source, lsource_mask) implicit none class(whm_pl), intent(inout) :: self !! WHM massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine whm_util_append_pl module subroutine whm_util_spill_pl(self, discards, lspill_list, ldestructive) diff --git a/src/rmvs/rmvs_util.f90 b/src/rmvs/rmvs_util.f90 index 9f9cf0037..ee9ce6932 100644 --- a/src/rmvs/rmvs_util.f90 +++ b/src/rmvs/rmvs_util.f90 @@ -11,21 +11,23 @@ module subroutine rmvs_util_append_pl(self, source, lsource_mask) !! Arguments class(rmvs_pl), intent(inout) :: self !! RMVS massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to select type(source) class is (rmvs_pl) - call whm_util_append_pl(self, source, lsource_mask) + associate(nold => self%nbody, nsrc => source%nbody) + call util_append(self%nenc, source%nenc, nold, nsrc, lsource_mask) + call util_append(self%tpenc1P, source%tpenc1P, nold, nsrc, lsource_mask) + call util_append(self%plind, source%plind, nold, nsrc, lsource_mask) - call util_append(self%nenc, source%nenc, lsource_mask) - call util_append(self%tpenc1P, source%tpenc1P, lsource_mask) - call util_append(self%plind, source%plind, lsource_mask) + ! The following are not implemented as RMVS doesn't make use of fill operations on pl type + ! So they are here as a placeholder in case someone wants to extend the RMVS class for some reason + !call util_append(self%outer, source%outer, nold, nsrc, lsource_mask) + !call util_append(self%inner, source%inner, nold, nsrc, lsource_mask) + !call util_append(self%planetocentric, source%planetocentric, nold, nsrc, lsource_mask) - ! The following are not implemented as RMVS doesn't make use of fill operations on pl type - ! So they are here as a placeholder in case someone wants to extend the RMVS class for some reason - !call util_append(self%outer, source%outer, lsource_mask) - !call util_append(self%inner, source%inner, lsource_mask) - !call util_append(self%planetocentric, source%planetocentric, lsource_mask) + call whm_util_append_pl(self, source, lsource_mask) + end associate class default write(*,*) "Invalid object passed to the append method. Source must be of class rmvs_pl or its descendents!" call util_exit(FAILURE) @@ -44,15 +46,17 @@ module subroutine rmvs_util_append_tp(self, source, lsource_mask) !! Arguments class(rmvs_tp), intent(inout) :: self !! RMVS test particle object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to select type(source) class is (rmvs_tp) - call util_append_tp(self, source, lsource_mask) ! Note: whm_tp does not have its own append method, so we skip back to the base class + associate(nold => self%nbody, nsrc => source%nbody) + call util_append(self%lperi, source%lperi, nold, nsrc, lsource_mask) + call util_append(self%plperP, source%plperP, nold, nsrc, lsource_mask) + call util_append(self%plencP, source%plencP, nold, nsrc, lsource_mask) - call util_append(self%lperi, source%lperi, lsource_mask) - call util_append(self%plperP, source%plperP, lsource_mask) - call util_append(self%plencP, source%plencP, lsource_mask) + call util_append_tp(self, source, lsource_mask) ! Note: whm_tp does not have its own append method, so we skip back to the base class + end associate class default write(*,*) "Invalid object passed to the append method. Source must be of class rmvs_tp or its descendents!" call util_exit(FAILURE) @@ -139,8 +143,6 @@ module subroutine rmvs_util_resize_pl(self, nnew) class(rmvs_pl), intent(inout) :: self !! RMVS massive body object integer(I4B), intent(in) :: nnew !! New size neded - call whm_util_resize_pl(self, nnew) - call util_resize(self%nenc, nnew) call util_resize(self%tpenc1P, nnew) call util_resize(self%plind, nnew) @@ -151,6 +153,7 @@ module subroutine rmvs_util_resize_pl(self, nnew) !call util_resize(self%inner, nnew) !call util_resize(self%planetocentric, nnew) + call whm_util_resize_pl(self, nnew) return end subroutine rmvs_util_resize_pl @@ -164,13 +167,13 @@ module subroutine rmvs_util_resize_tp(self, nnew) class(rmvs_tp), intent(inout) :: self !! RMVS test particle object integer(I4B), intent(in) :: nnew !! New size neded - call util_resize_tp(self, nnew) - call util_resize(self%lperi, nnew) call util_resize(self%plperP, nnew) call util_resize(self%plencP, nnew) call util_resize(self%xheliocentric, nnew) + call util_resize_tp(self, nnew) + return end subroutine rmvs_util_resize_tp diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 5f6d3926a..253fb2700 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -271,7 +271,6 @@ subroutine symba_discard_peri_pl(pl, system, param) pl%lfirst = lfirst_orig return - end subroutine symba_discard_peri_pl @@ -285,6 +284,8 @@ module subroutine symba_discard_pl(self, system, param) class(symba_pl), intent(inout) :: self !! SyMBA test particle object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + real(DP) :: Eorbit_before, Eorbit_after select type(system) class is (symba_nbody_system) @@ -309,8 +310,17 @@ module subroutine symba_discard_pl(self, system, param) end if if (any(pl%ldiscard(:))) then + if (param%lenergy) then + call system%get_energy_and_momentum(param) + Eorbit_before = system%te + end if call symba_discard_nonplpl_conservation(self, system, param) call pl%rearray(system, param) + if (param%lenergy) then + call system%get_energy_and_momentum(param) + Eorbit_after = system%te + system%Ecollisions = Eorbit_after - Eorbit_before + end if end if end associate diff --git a/src/symba/symba_fragmentation.f90 b/src/symba/symba_fragmentation.f90 index a27d31e12..b36c54e9a 100644 --- a/src/symba/symba_fragmentation.f90 +++ b/src/symba/symba_fragmentation.f90 @@ -28,124 +28,63 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v, real(DP), dimension(2) :: vol real(DP), dimension(:, :), allocatable :: vb_frag, xb_frag, rot_frag, Ip_frag real(DP), dimension(:), allocatable :: m_frag, rad_frag + integer(I4B), dimension(:), allocatable :: id_frag logical :: lfailure - logical, dimension(system%pl%nbody) :: lmask - class(symba_pl), allocatable :: plnew - select type(pl => system%pl) - class is (symba_pl) - select type(pl_discards => system%pl_discards) - class is (symba_merger) - associate(pl_adds => system%pl_adds, cb => system%cb) - ! Collisional fragments will be uniformly distributed around the pre-impact barycenter - nfrag = NFRAG_DISRUPT - allocate(m_frag(nfrag)) - allocate(rad_frag(nfrag)) - allocate(xb_frag(NDIM, nfrag)) - allocate(vb_frag(NDIM, nfrag)) - allocate(rot_frag(NDIM, nfrag)) - allocate(Ip_frag(NDIM, nfrag)) - - mtot = sum(mass(:)) - xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot - vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot - - ! Get mass weighted mean of Ip and average density - Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mtot - vol(:) = 4._DP / 3._DP * PI * radius(:)**3 - avg_dens = mtot / sum(vol(:)) - - ! Distribute the mass among fragments, with a branch to check for the size of the second largest fragment - m_frag(1) = mass_res(1) - if (mass_res(2) > mass_res(1) / 3._DP) then - m_frag(2) = mass_res(2) - istart = 3 - else - istart = 2 - end if - ! Distribute remaining mass among the remaining bodies - do i = istart, nfrag - m_frag(i) = (mtot - sum(m_frag(1:istart - 1))) / (nfrag - istart + 1) - end do - - ! Distribute any residual mass if there is any and set the radius - m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) - rad_frag(:) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) - - do i = 1, nfrag - Ip_frag(:, i) = Ip_new(:) - end do - - call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & - nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) - - if (lfailure) then - write(*,*) 'No fragment solution found, so treat as a pure hit-and-run' - status = ACTIVE - nfrag = 0 - else - ! Populate the list of new bodies - write(*,'("Generating ",I2.0," fragments")') nfrag - status = DISRUPTION - - ! Add the family bodies to the subtraction list - nfamily = size(family(:)) - lmask(:) = .false. - lmask(family(:)) = .true. - pl%status(family(:)) = MERGED - nstart = pl_discards%nbody + 1 - nend = pl_discards%nbody + nfamily - call pl_discards%append(pl, lmask) - ! Record how many bodies were subtracted in this event - pl_discards%ncomp(nstart:nend) = nfamily - - allocate(plnew, mold=pl) - call plnew%setup(nfrag, param) - - plnew%id(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)] - system%maxid = system%maxid + nfrag - plnew%status(:) = ACTIVE - plnew%lcollision(:) = .false. - plnew%ldiscard(:) = .false. - plnew%xb(:,:) = xb_frag(:, :) - plnew%vb(:,:) = vb_frag(:, :) - do i = 1, nfrag - plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) - plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) - end do - plnew%mass(:) = m_frag(:) - plnew%Gmass(:) = param%GU * m_frag(:) - plnew%density(:) = avg_dens - plnew%radius(:) = rad_frag(:) - plnew%info(:)%origin_type = "Disruption" - plnew%info(:)%origin_time = param%t - do i = 1, nfrag - plnew%info(i)%origin_xh(:) = plnew%xh(:,i) - plnew%info(i)%origin_vh(:) = plnew%vh(:,i) - end do - if (param%lrotation) then - plnew%Ip(:,:) = Ip_frag(:,:) - plnew%rot(:,:) = rot_frag(:,:) - end if - if (param%ltides) then - ibiggest = maxloc(pl%Gmass(family(:)), dim=1) - plnew%Q = pl%Q(ibiggest) - plnew%k2 = pl%k2(ibiggest) - plnew%tlag = pl%tlag(ibiggest) - end if - - ! Append the new merged body to the list and record how many we made - nstart = pl_adds%nbody + 1 - nend = pl_adds%nbody + plnew%nbody - call pl_adds%append(plnew) - pl_adds%ncomp(nstart:nend) = plnew%nbody - - call plnew%setup(0, param) - deallocate(plnew) - end if - end associate - end select - end select + ! Collisional fragments will be uniformly distributed around the pre-impact barycenter + nfrag = NFRAG_DISRUPT + allocate(m_frag(nfrag)) + allocate(rad_frag(nfrag)) + allocate(xb_frag(NDIM, nfrag)) + allocate(vb_frag(NDIM, nfrag)) + allocate(rot_frag(NDIM, nfrag)) + allocate(Ip_frag(NDIM, nfrag)) + allocate(id_frag(nfrag)) + + mtot = sum(mass(:)) + xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot + vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot + + ! Get mass weighted mean of Ip and average density + Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mtot + vol(:) = 4._DP / 3._DP * PI * radius(:)**3 + avg_dens = mtot / sum(vol(:)) + + ! Distribute the mass among fragments, with a branch to check for the size of the second largest fragment + m_frag(1) = mass_res(1) + if (mass_res(2) > mass_res(1) / 3._DP) then + m_frag(2) = mass_res(2) + istart = 3 + else + istart = 2 + end if + ! Distribute remaining mass among the remaining bodies + do i = istart, nfrag + m_frag(i) = (mtot - sum(m_frag(1:istart - 1))) / (nfrag - istart + 1) + end do + + ! Distribute any residual mass if there is any and set the radius + m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) + rad_frag(:) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) + id_frag(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)] + + do i = 1, nfrag + Ip_frag(:, i) = Ip_new(:) + end do + + call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) + + if (lfailure) then + write(*,*) 'No fragment solution found, so treat as a pure hit-and-run' + status = ACTIVE + nfrag = 0 + else + ! Populate the list of new bodies + write(*,'("Generating ",I2.0," fragments")') nfrag + status = DISRUPTION + call symba_fragmentation_mergeaddsub(system, param, family, id_frag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, status) + end if return end function symba_fragmentation_casedisruption @@ -168,7 +107,7 @@ module function symba_fragmentation_casehitandrun(system, param, family, x, v, m ! Result integer(I4B) :: status !! Status flag assigned to this outcome ! Internals - integer(I4B) :: i, nfrag, jproj, jtarg, idstart, ibiggest, nfamily, nstart, nend + integer(I4B) :: i, nfrag, jproj, jtarg, idstart, ibiggest, nfamily real(DP) :: mtot, avg_dens real(DP), dimension(NDIM) :: xcom, vcom real(DP), dimension(2) :: vol @@ -177,136 +116,73 @@ module function symba_fragmentation_casehitandrun(system, param, family, x, v, m integer(I4B), dimension(:), allocatable :: id_frag logical :: lpure logical, dimension(system%pl%nbody) :: lmask - class(symba_pl), allocatable :: plnew - select type(pl => system%pl) - class is (symba_pl) - select type(pl_discards => system%pl_discards) - class is (symba_merger) - associate(pl_adds => system%pl_adds, cb => system%cb) - mtot = sum(mass(:)) - xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot - vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot - lpure = .false. - - ! The largest body will stay untouched - if (mass(1) > mass(2)) then - jtarg = 1 - jproj = 2 - else - jtarg = 2 - jproj = 1 - end if - - if (mass_res(2) > 0.9_DP * mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched - write(*,*) 'Pure hit and run. No new fragments generated.' - nfrag = 0 - lpure = .true. - else ! Imperfect hit and run, so we'll keep the largest body and destroy the other - nfrag = NFRAG_DISRUPT - 1 - lpure = .false. - allocate(m_frag(nfrag)) - allocate(id_frag(nfrag)) - allocate(rad_frag(nfrag)) - allocate(xb_frag(NDIM, nfrag)) - allocate(vb_frag(NDIM, nfrag)) - allocate(rot_frag(NDIM, nfrag)) - allocate(Ip_frag(NDIM, nfrag)) - m_frag(1) = mass(jtarg) - ibiggest = maxloc(pl%Gmass(family(:)), dim=1) - id_frag(1) = pl%id(ibiggest) - rad_frag(1) = radius(jtarg) - xb_frag(:, 1) = x(:, jtarg) - vb_frag(:, 1) = v(:, jtarg) - Ip_frag(:,1) = Ip(:, jtarg) - - ! Get mass weighted mean of Ip and average density - vol(:) = 4._DP / 3._DP * pi * radius(:)**3 - avg_dens = mass(jproj) / vol(jproj) - m_frag(2:nfrag) = (mtot - m_frag(1)) / (nfrag - 1) - rad_frag(2:nfrag) = (3 * m_frag(2:nfrag) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) - m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) - - do i = 1, nfrag - Ip_frag(:, i) = Ip(:, jproj) - end do - - ! Put the fragments on the circle surrounding the center of mass of the system - call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & - nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lpure) - if (lpure) then - write(*,*) 'Should have been a pure hit and run instead' - nfrag = 0 - else - write(*,'("Generating ",I2.0," fragments")') nfrag - end if - end if - if (lpure) then - status = ACTIVE - else - status = HIT_AND_RUN - - ! Add the family bodies to the subtraction list - nfamily = size(family(:)) - lmask(:) = .false. - lmask(family(:)) = .true. - pl%status(family(:)) = MERGED - nstart = pl_discards%nbody + 1 - nend = pl_discards%nbody + nfamily - call pl_discards%append(pl, lmask) - ! Record how many bodies were subtracted in this event - pl_discards%ncomp(nstart:nend) = nfamily - - allocate(plnew, mold=pl) - call plnew%setup(nfrag, param) - - plnew%id(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)] - system%maxid = system%maxid + nfrag - plnew%status(:) = ACTIVE - plnew%lcollision(:) = .false. - plnew%ldiscard(:) = .false. - plnew%xb(:,:) = xb_frag(:, :) - plnew%vb(:,:) = vb_frag(:, :) - do i = 1, nfrag - plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) - plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) - end do - plnew%mass(:) = m_frag(:) - plnew%Gmass(:) = param%GU * m_frag(:) - plnew%density(:) = avg_dens - plnew%radius(:) = rad_frag(:) - plnew%info(:)%origin_type = "Hit and run fragment" - plnew%info(:)%origin_time = param%t - do i = 1, nfrag - plnew%info(i)%origin_xh(:) = plnew%xh(:,i) - plnew%info(i)%origin_vh(:) = plnew%vh(:,i) - end do - if (param%lrotation) then - plnew%Ip(:,:) = Ip_frag(:,:) - plnew%rot(:,:) = rot_frag(:,:) - end if - if (param%ltides) then - ibiggest = maxloc(pl%Gmass(family(:)), dim=1) - plnew%Q = pl%Q(ibiggest) - plnew%k2 = pl%k2(ibiggest) - plnew%tlag = pl%tlag(ibiggest) - end if - - ! Append the new merged body to the list and record how many we made - nstart = pl_adds%nbody + 1 - nend = pl_adds%nbody + plnew%nbody - call pl_adds%append(plnew) - pl_adds%ncomp(nstart:nend) = plnew%nbody - - call plnew%setup(0, param) - deallocate(plnew) - - end if - end associate - end select - end select + mtot = sum(mass(:)) + xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot + vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot + lpure = .false. + + ! The largest body will stay untouched + if (mass(1) > mass(2)) then + jtarg = 1 + jproj = 2 + else + jtarg = 2 + jproj = 1 + end if + + if (mass_res(2) > 0.9_DP * mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched + write(*,*) 'Pure hit and run. No new fragments generated.' + nfrag = 0 + lpure = .true. + else ! Imperfect hit and run, so we'll keep the largest body and destroy the other + nfrag = NFRAG_DISRUPT - 1 + lpure = .false. + allocate(m_frag(nfrag)) + allocate(id_frag(nfrag)) + allocate(rad_frag(nfrag)) + allocate(xb_frag(NDIM, nfrag)) + allocate(vb_frag(NDIM, nfrag)) + allocate(rot_frag(NDIM, nfrag)) + allocate(Ip_frag(NDIM, nfrag)) + m_frag(1) = mass(jtarg) + ibiggest = maxloc(system%pl%Gmass(family(:)), dim=1) + id_frag(1) = system%pl%id(ibiggest) + rad_frag(1) = radius(jtarg) + xb_frag(:, 1) = x(:, jtarg) + vb_frag(:, 1) = v(:, jtarg) + Ip_frag(:,1) = Ip(:, jtarg) + + ! Get mass weighted mean of Ip and average density + vol(:) = 4._DP / 3._DP * pi * radius(:)**3 + avg_dens = mass(jproj) / vol(jproj) + m_frag(2:nfrag) = (mtot - m_frag(1)) / (nfrag - 1) + rad_frag(2:nfrag) = (3 * m_frag(2:nfrag) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) + m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) + id_frag(2:nfrag) = [(i, i = system%maxid + 1, system%maxid + nfrag - 1)] + + do i = 1, nfrag + Ip_frag(:, i) = Ip(:, jproj) + end do + + ! Put the fragments on the circle surrounding the center of mass of the system + call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lpure) + if (lpure) then + write(*,*) 'Should have been a pure hit and run instead' + nfrag = 0 + else + write(*,'("Generating ",I2.0," fragments")') nfrag + end if + end if + if (lpure) then + status = ACTIVE + else + status = HIT_AND_RUN + call symba_fragmentation_mergeaddsub(system, param, family, id_frag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, status) + end if - return + return end function symba_fragmentation_casehitandrun @@ -328,120 +204,75 @@ module function symba_fragmentation_casemerge(system, param, family, x, v, mass, ! Result integer(I4B) :: status !! Status flag assigned to this outcome ! Internals - integer(I4B) :: i, j, ibiggest, nfamily, nstart, nend - real(DP) :: mass_new, radius_new, volume_new, pe - real(DP), dimension(NDIM) :: xcom, vcom, xc, vc, xcrossv + integer(I4B) :: i, j, ibiggest, nfamily + real(DP) :: volume_new, pe + real(DP), dimension(NDIM) :: xc, vc, xcrossv real(DP), dimension(2) :: vol real(DP), dimension(NDIM) :: L_orb_old, L_spin_old - real(DP), dimension(NDIM) :: L_spin_new, rot_new, Ip_new + real(DP), dimension(NDIM) :: L_spin_new logical, dimension(system%pl%nbody) :: lmask - class(symba_pl), allocatable :: plnew + real(DP), dimension(NDIM, 1) :: vb_frag, xb_frag, rot_frag, Ip_frag + real(DP), dimension(1) :: m_frag, rad_frag + integer(I4B), dimension(1) :: id_frag select type(pl => system%pl) class is (symba_pl) - select type(pl_discards => system%pl_discards) - class is (symba_merger) - associate(pl_adds => system%pl_adds, cb => system%cb) - status = MERGED - write(*, '("Merging bodies ",99(I8,",",:))') pl%id(family(:)) - mass_new = sum(mass(:)) - - ! Merged body is created at the barycenter of the original bodies - xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mass_new - vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mass_new - - ! Get mass weighted mean of Ip and - vol(:) = 4._DP / 3._DP * PI * radius(:)**3 - volume_new = sum(vol(:)) - radius_new = (3 * volume_new / (4 * PI))**(1._DP / 3._DP) - - L_orb_old(:) = 0.0_DP - - ! Compute orbital angular momentum of pre-impact system - do i = 1, 2 - xc(:) = x(:, i) - xcom(:) - vc(:) = v(:, i) - vcom(:) - xcrossv(:) = xc(:) .cross. vc(:) - L_orb_old(:) = L_orb_old(:) + mass(i) * xcrossv(:) - end do - - if (param%lrotation) then - Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mass_new - L_spin_old(:) = L_spin(:,1) + L_spin(:,2) + write(*, '("Merging bodies ",99(I8,",",:))') pl%id(family(:)) - ! Conserve angular momentum by putting pre-impact orbital momentum into spin of the new body - L_spin_new(:) = L_orb_old(:) + L_spin_old(:) - - ! Assume prinicpal axis rotation on 3rd Ip axis - rot_new(:) = L_spin_new(:) / (Ip_new(3) * mass_new * radius_new**2) - else ! If spin is not enabled, we will consider the lost pre-collision angular momentum as "escaped" and add it to our bookkeeping variable - system%Lescape(:) = system%Lescape(:) + L_orb_old(:) - end if - - ! Keep track of the component of potential energy due to the pre-impact family for book-keeping - nfamily = size(family(:)) - pe = 0.0_DP - do j = 1, nfamily - do i = j + 1, nfamily - pe = pe - pl%mass(i) * pl%mass(j) / norm2(pl%xb(:, i) - pl%xb(:, j)) - end do - end do - system%Ecollisions = system%Ecollisions + pe - system%Euntracked = system%Euntracked - pe - - ! Add the family bodies to the subtraction list - lmask(:) = .false. - lmask(family(:)) = .true. - pl%status(family(:)) = MERGED - nstart = pl_discards%nbody + 1 - nend = pl_discards%nbody + nfamily - call pl_discards%append(pl, lmask) - ! Record how many bodies were subtracted in this event - pl_discards%ncomp(nstart:nend) = nfamily + ibiggest = maxloc(pl%Gmass(family(:)), dim=1) + id_frag(1) = pl%id(family(ibiggest)) - ! Create the new merged body - allocate(plnew, mold=pl) - call plnew%setup(1, param) + m_frag(1) = sum(mass(:)) + + ! Merged body is created at the barycenter of the original bodies + xb_frag(:,1) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / m_frag(1) + vb_frag(:,1) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / m_frag(1) + + ! Get mass weighted mean of Ip and + vol(:) = 4._DP / 3._DP * PI * radius(:)**3 + volume_new = sum(vol(:)) + rad_frag(1) = (3 * volume_new / (4 * PI))**(1._DP / 3._DP) - ! The merged body's name will be that of the largest of the two parents - ibiggest = maxloc(pl%Gmass(family(:)), dim=1) - plnew%id(1) = pl%id(family(ibiggest)) - plnew%status(1) = ACTIVE - plnew%lcollision = .false. - plnew%ldiscard = .false. - plnew%xb(:,1) = xcom(:) - plnew%vb(:,1) = vcom(:) - plnew%xh(:,1) = xcom(:) - cb%xb(:) - plnew%vh(:,1) = vcom(:) - cb%vb(:) - plnew%mass(1) = mass_new - plnew%Gmass(1) = param%GU * mass_new - plnew%density(1) = mass_new / volume_new - plnew%radius(1) = radius_new - plnew%info(1) = pl%info(family(ibiggest)) - if (param%lrotation) then - plnew%Ip(:,1) = Ip_new(:) - plnew%rot(:,1) = rot_new(:) - end if - if (param%ltides) then - plnew%Q = pl%Q(ibiggest) - plnew%k2 = pl%k2(ibiggest) - plnew%tlag = pl%tlag(ibiggest) - end if + L_orb_old(:) = 0.0_DP - ! Append the new merged body to the list and record how many we made - nstart = pl_adds%nbody + 1 - nend = pl_adds%nbody + plnew%nbody - call pl_adds%append(plnew) - pl_adds%ncomp(nstart:nend) = plnew%nbody + ! Compute orbital angular momentum of pre-impact system + do i = 1, 2 + xc(:) = x(:, i) - xb_frag(:,1) + vc(:) = v(:, i) - vb_frag(:,1) + xcrossv(:) = xc(:) .cross. vc(:) + L_orb_old(:) = L_orb_old(:) + mass(i) * xcrossv(:) + end do + + if (param%lrotation) then + Ip_frag(:,1) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / m_frag(1) + L_spin_old(:) = L_spin(:,1) + L_spin(:,2) - call plnew%setup(0, param) - deallocate(plnew) - end associate - end select + ! Conserve angular momentum by putting pre-impact orbital momentum into spin of the new body + L_spin_new(:) = L_orb_old(:) + L_spin_old(:) + + ! Assume prinicpal axis rotation on 3rd Ip axis + rot_frag(:,1) = L_spin_new(:) / (Ip_frag(3,1) * m_frag(1) * rad_frag(1)**2) + else ! If spin is not enabled, we will consider the lost pre-collision angular momentum as "escaped" and add it to our bookkeeping variable + system%Lescape(:) = system%Lescape(:) + L_orb_old(:) + end if + + ! Keep track of the component of potential energy due to the pre-impact family for book-keeping + nfamily = size(family(:)) + pe = 0.0_DP + do j = 1, nfamily + do i = j + 1, nfamily + pe = pe - pl%mass(i) * pl%mass(j) / norm2(pl%xb(:, i) - pl%xb(:, j)) + end do + end do + system%Ecollisions = system%Ecollisions + pe + system%Euntracked = system%Euntracked - pe + + status = MERGED + call symba_fragmentation_mergeaddsub(system, param, family, id_frag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, status) + end select return - end function symba_fragmentation_casemerge @@ -469,122 +300,171 @@ module function symba_fragmentation_casesupercatastrophic(system, param, family, real(DP), dimension(NDIM) :: Ip_new real(DP), dimension(:, :), allocatable :: vb_frag, xb_frag, rot_frag, Ip_frag real(DP), dimension(:), allocatable :: m_frag, rad_frag + integer(I4B), dimension(:), allocatable :: id_frag logical :: lfailure logical, dimension(system%pl%nbody) :: lmask + + ! Collisional fragments will be uniformly distributed around the pre-impact barycenter + nfrag = NFRAG_SUPERCAT + allocate(m_frag(nfrag)) + allocate(rad_frag(nfrag)) + allocate(id_frag(nfrag)) + allocate(xb_frag(NDIM, nfrag)) + allocate(vb_frag(NDIM, nfrag)) + allocate(rot_frag(NDIM, nfrag)) + allocate(Ip_frag(NDIM, nfrag)) + + mtot = sum(mass(:)) + xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot + vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot + + ! Get mass weighted mean of Ip and average density + Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mtot + vol(:) = 4._DP / 3._DP * pi * radius(:)**3 + avg_dens = mtot / sum(vol(:)) + + ! If we are adding the first and largest fragment (lr), check to see if its mass is SMALLER than an equal distribution of + ! mass between all fragments. If so, we will just distribute the mass equally between the fragments + min_frag_mass = mtot / nfrag + if (mass_res(1) < min_frag_mass) then + m_frag(:) = min_frag_mass + else + m_frag(1) = mass_res(1) + m_frag(2:nfrag) = (mtot - mass_res(1)) / (nfrag - 1) + end if + ! Distribute any residual mass if there is any and set the radius + m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) + rad_frag(:) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) + id_frag(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)] + + do i = 1, nfrag + Ip_frag(:, i) = Ip_new(:) + end do + + call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) + + if (lfailure) then + write(*,*) 'No fragment solution found, so treat as a pure hit-and-run' + status = ACTIVE + nfrag = 0 + else + ! Populate the list of new bodies + write(*,'("Generating ",I2.0," fragments")') nfrag + status = SUPERCATASTROPHIC + call symba_fragmentation_mergeaddsub(system, param, family, id_frag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, status) + end if + + return + end function symba_fragmentation_casesupercatastrophic + + + subroutine symba_fragmentation_mergeaddsub(system, param, family, id_frag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, status) + !! author: David A. Minton + !! + !! Fills the pl_discards and pl_adds with removed and added bodies + !! + implicit none + ! Arguments + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision + integer(I4B), dimension(:), intent(in) :: id_frag !! List of fragment ids + real(DP), dimension(:), intent(in) :: m_frag, rad_frag !! Distribution of fragment mass and radii + real(DP), dimension(:,:), intent(in) :: Ip_frag !! Fragment rotational inertia vectors + real(DP), dimension(:,:), intent(in) :: xb_frag, vb_frag, rot_frag !! Fragment barycentric position, barycentric velocity, and rotation vectors + integer(I4B), intent(in) :: status !! Status flag to assign to adds + ! Internals + integer(I4B) :: i, ibiggest, nstart, nend, nfamily, nfrag + logical, dimension(system%pl%nbody) :: lmask class(symba_pl), allocatable :: plnew - + select type(pl => system%pl) class is (symba_pl) select type(pl_discards => system%pl_discards) class is (symba_merger) associate(pl_adds => system%pl_adds, cb => system%cb) - ! Collisional fragments will be uniformly distributed around the pre-impact barycenter - nfrag = NFRAG_SUPERCAT - allocate(m_frag(nfrag)) - allocate(rad_frag(nfrag)) - allocate(xb_frag(NDIM, nfrag)) - allocate(vb_frag(NDIM, nfrag)) - allocate(rot_frag(NDIM, nfrag)) - allocate(Ip_frag(NDIM, nfrag)) - - mtot = sum(mass(:)) - xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot - vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot - - ! Get mass weighted mean of Ip and average density - Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mtot - vol(:) = 4._DP / 3._DP * pi * radius(:)**3 - avg_dens = mtot / sum(vol(:)) - - ! If we are adding the first and largest fragment (lr), check to see if its mass is SMALLER than an equal distribution of - ! mass between all fragments. If so, we will just distribute the mass equally between the fragments - min_frag_mass = mtot / nfrag - if (mass_res(1) < min_frag_mass) then - m_frag(:) = min_frag_mass - else - m_frag(1) = mass_res(1) - m_frag(2:nfrag) = (mtot - mass_res(1)) / (nfrag - 1) - end if - ! Distribute any residual mass if there is any and set the radius - m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) - rad_frag(:) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) - + + ! Add the family bodies to the subtraction list + nfamily = size(family(:)) + nfrag = size(m_frag(:)) + lmask(:) = .false. + lmask(family(:)) = .true. + pl%status(family(:)) = MERGED + nstart = pl_discards%nbody + 1 + nend = pl_discards%nbody + nfamily + call pl_discards%append(pl, lmask) + pl%ldiscard(family(:)) = .true. + pl%lcollision(family(:)) = .true. + + ! Record how many bodies were subtracted in this event + pl_discards%ncomp(nstart:nend) = nfamily + + ! Setup new bodies + allocate(plnew, mold=pl) + call plnew%setup(nfrag, param) + ibiggest = maxloc(pl%Gmass(family(:)), dim=1) + + plnew%id(:) = id_frag(:) + system%maxid = system%maxid + nfrag + plnew%status(:) = ACTIVE + plnew%lcollision(:) = .false. + plnew%ldiscard(:) = .false. + plnew%xb(:,:) = xb_frag(:, :) + plnew%vb(:,:) = vb_frag(:, :) do i = 1, nfrag - Ip_frag(:, i) = Ip_new(:) + plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) + plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) end do + plnew%mass(:) = m_frag(:) + plnew%Gmass(:) = param%GU * m_frag(:) + plnew%radius(:) = rad_frag(:) + plnew%density(:) = m_frag(:) / rad_frag(:) - call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & - nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) - - if (lfailure) then - write(*,*) 'No fragment solution found, so treat as a pure hit-and-run' - status = ACTIVE - nfrag = 0 - else - ! Populate the list of new bodies - write(*,'("Generating ",I2.0," fragments")') nfrag - status = SUPERCATASTROPHIC - - ! Add the family bodies to the subtraction list - nfamily = size(family(:)) - lmask(:) = .false. - lmask(family(:)) = .true. - pl%status(family(:)) = MERGED - nstart = pl_discards%nbody + 1 - nend = pl_discards%nbody + nfamily - call pl_discards%append(pl, lmask) - ! Record how many bodies were subtracted in this event - pl_discards%ncomp(nstart:nend) = nfamily - - allocate(plnew, mold=pl) - call plnew%setup(nfrag, param) - - plnew%id(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)] - system%maxid = system%maxid + nfrag - plnew%status(:) = ACTIVE - plnew%lcollision(:) = .false. - plnew%ldiscard(:) = .false. - plnew%xb(:,:) = xb_frag(:, :) - plnew%vb(:,:) = vb_frag(:, :) - do i = 1, nfrag - plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) - plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) - end do - plnew%mass(:) = m_frag(:) - plnew%Gmass(:) = param%GU * m_frag(:) - plnew%density(:) = avg_dens - plnew%radius(:) = rad_frag(:) + select case(status) + case(DISRUPTION) + plnew%info(:)%origin_type = "Disruption" + case(SUPERCATASTROPHIC) plnew%info(:)%origin_type = "Supercatastrophic" + case(HIT_AND_RUN) + plnew%info(:)%origin_type = "Hit and run fragment" + case(MERGED) + plnew%info(1) = pl%info(ibiggest) + end select + + if (status /= MERGED) then plnew%info(:)%origin_time = param%t do i = 1, nfrag plnew%info(i)%origin_xh(:) = plnew%xh(:,i) plnew%info(i)%origin_vh(:) = plnew%vh(:,i) end do - if (param%lrotation) then - plnew%Ip(:,:) = Ip_frag(:,:) - plnew%rot(:,:) = rot_frag(:,:) - end if - if (param%ltides) then - ibiggest = maxloc(pl%Gmass(family(:)), dim=1) - plnew%Q = pl%Q(ibiggest) - plnew%k2 = pl%k2(ibiggest) - plnew%tlag = pl%tlag(ibiggest) - end if - - ! Append the new merged body to the list and record how many we made - nstart = pl_adds%nbody + 1 - nend = pl_adds%nbody + plnew%nbody - call pl_adds%append(plnew) - pl_adds%ncomp(nstart:nend) = plnew%nbody - - call plnew%setup(0, param) - deallocate(plnew) end if + + if (param%lrotation) then + plnew%Ip(:,:) = Ip_frag(:,:) + plnew%rot(:,:) = rot_frag(:,:) + end if + if (param%ltides) then + plnew%Q = pl%Q(ibiggest) + plnew%k2 = pl%k2(ibiggest) + plnew%tlag = pl%tlag(ibiggest) + end if + call plnew%set_mu(cb) + pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY + + ! Append the new merged body to the list and record how many we made + nstart = pl_adds%nbody + 1 + nend = pl_adds%nbody + plnew%nbody + call pl_adds%append(plnew, lsource_mask=[(.true., i=1, plnew%nbody)]) + pl_adds%ncomp(nstart:nend) = plnew%nbody + + call plnew%setup(0, param) + deallocate(plnew) end associate end select end select return - end function symba_fragmentation_casesupercatastrophic + end subroutine symba_fragmentation_mergeaddsub end submodule s_symba_fragmentation diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index ba972db8b..97e5798c2 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -2,26 +2,88 @@ use swiftest contains - module subroutine symba_io_dump_particle_info(self, param, msg) + module subroutine symba_io_dump_particle_info(system, param, lincludecb, tpidx, plidx) !! author: David A. Minton !! - !! Dumps the particle information data to a file + !! Dumps the particle information data to a file. + !! Pass a list of array indices for test particles (tpidx) and/or massive bodies (plidx) to append implicit none - class(symba_particle_info), intent(inout) :: self !! Swiftest base object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - character(*), optional, intent(in) :: msg !! Message to display with dump operation - end subroutine symba_io_dump_particle_info + ! Arguments + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA extensions + logical, optional, intent(in) :: lincludecb !! Set to true to include the central body (default is false) + integer(I4B), dimension(:), optional, intent(in) :: tpidx !! Array of test particle indices to append to the particle file + integer(I4B), dimension(:), optional, intent(in) :: plidx !! Array of massive body indices to append to the particle file + ! Internals + logical, save :: lfirst = .true. + integer(I4B), parameter :: LUN = 22 + integer(I4B) :: i, ierr + if (lfirst) then + select case(param%out_stat) + case('APPEND') + open(unit = LUN, file = param%particle_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr) + case('NEW', 'UNKNOWN', 'REPLACE') + open(unit = LUN, file = param%particle_out, status = param%out_stat, form = 'UNFORMATTED', iostat = ierr) + case default + write(*,*) 'Invalid status code',trim(adjustl(param%out_stat)) + call util_exit(FAILURE) + end select + if (ierr /= 0) then + write(*, *) "Swiftest error:" + write(*, *) " particle output file already exists or cannot be accessed" + call util_exit(FAILURE) + end if - module subroutine symba_io_initialize_particle_info(self, param) - !! author: David A. Minton - !! - !! Initializes a particle info data structure, either starting a new one or reading one in - !! from a file if it is a restarted run - implicit none - class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine symba_io_initialize_particle_info + lfirst = .false. + else + open(unit = LUN, file = param%particle_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr) + if (ierr /= 0) then + write(*, *) "Swiftest error:" + write(*, *) " unable to open binary output file for APPEND" + call util_exit(FAILURE) + end if + end if + + if (present(lincludecb)) then + if (lincludecb) then + select type(cb => system%cb) + class is (symba_cb) + write(LUN) cb%id + write(LUN) cb%info + end select + end if + end if + + if (present(plidx) .and. (system%pl%nbody > 0)) then + select type(pl => system%pl) + class is (symba_pl) + do i = 1, size(plidx) + write(LUN) pl%id(plidx(i)) + write(LUN) pl%info(plidx(i)) + end do + end select + end if + + if (present(tpidx) .and. (system%tp%nbody > 0)) then + select type(tp => system%tp) + class is (symba_tp) + do i = 1, size(tpidx) + write(LUN) tp%id(tpidx(i)) + write(LUN) tp%info(tpidx(i)) + end do + end select + end if + + close(unit = LUN, iostat = ierr) + if (ierr /= 0) then + write(*, *) "Swiftest error:" + write(*, *) " unable to close particle output file" + call util_exit(FAILURE) + end if + + return + end subroutine symba_io_dump_particle_info module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) @@ -68,6 +130,8 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms ifirst = ilast + 1 param_value = io_get_token(line_trim, ifirst, ilast, iostat) select case (param_name) + case ("PARTICLE_OUT") + param%particle_out = param_value case ("FRAGMENTATION") call io_toupper(param_value) if (param_value == "YES" .or. param_value == "T") self%lfragmentation = .true. @@ -166,7 +230,7 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms ! Special handling is required for writing the random number seed array as its size is not known until runtime ! For the "SEED" parameter line, the first value will be the size of the seed array and the rest will be the seed array elements - write(param_name, Afmt) "PARTICLE_FILE"; write(param_value, Afmt) trim(adjustl(param%particle_file)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "PARTICLE_OUT"; write(param_value, Afmt) trim(adjustl(param%particle_out)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "GMTINY"; write(param_value, Rfmt) param%Gmtiny; write(unit, Afmt) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "FRAGMENTATION"; write(param_value, Lfmt) param%lfragmentation; write(unit, Afmt) adjustl(param_name), adjustl(param_value) if (param%lfragmentation) then @@ -194,19 +258,74 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms end subroutine symba_io_param_writer - module subroutine symba_io_read_frame_info(self, iu, param, form, ierr) + module subroutine symba_io_read_particle(system, param) !! author: David A. Minton !! - !! Reads a single frame of a particle info data from a file. + !! Reads an old particle information file for a restartd run implicit none - class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - character(*), intent(in) :: form !! Input format code ("XV" or "EL") - integer(I4B), intent(out) :: ierr !! Error code + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system file + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA extensions + + ! Internals + integer(I4B), parameter :: LUN = 22 + integer(I4B) :: i, ierr, id, idx + logical :: lmatch + type(symba_particle_info) :: tmpinfo + + open(unit = LUN, file = param%particle_out, status = 'OLD', form = 'UNFORMATTED', iostat = ierr) + if (ierr /= 0) then + write(*, *) "Swiftest error:" + write(*, *) " unable to open binary particle file for reading" + call util_exit(FAILURE) + end if + + select type(cb => system%cb) + class is (symba_cb) + select type(pl => system%pl) + class is (symba_pl) + select type(tp => system%tp) + class is (symba_tp) + do + lmatch = .false. + read(LUN, iostat=ierr) id + if (ierr /=0) exit + + if (idx == cb%id) then + read(LUN) cb%info + lmatch = .true. + else + if (pl%nbody > 0) then + idx = findloc(pl%id(:), id, dim=1) + if (idx /= 0) then + read(LUN) pl%info(idx) + lmatch = .true. + end if + end if + if (.not.lmatch .and. tp%nbody > 0) then + idx = findloc(tp%id(:), id, dim=1) + if (idx /= 0) then + read(LUN) tp%info(idx) + lmatch = .true. + end if + end if + end if + if (.not.lmatch) then + write(*,*) 'Particle id ',id,' not found. Skipping' + read(LUN) tmpinfo + end if + end do + close(unit = LUN, iostat = ierr) + end select + end select + end select + if (ierr /= 0) then + write(*, *) "Swiftest error:" + write(*, *) " unable to close particle output file" + call util_exit(FAILURE) + end if - ierr = 0 - end subroutine symba_io_read_frame_info + return + end subroutine symba_io_read_particle module subroutine symba_io_write_discard(self, param) @@ -253,7 +372,7 @@ module subroutine symba_io_write_discard(self, param) nsub = pl_discards%ncomp(isub) do j = 1, nadd if (iadd <= pl_adds%nbody) then - write(LUN, NAMEFMT) ADD, pl_discards%id(iadd), pl_discards%status(iadd) + write(LUN, NAMEFMT) ADD, pl_adds%id(iadd), pl_adds%status(iadd) write(LUN, VECFMT) pl_adds%xh(1, iadd), pl_adds%xh(2, iadd), pl_adds%xh(3, iadd) write(LUN, VECFMT) pl_adds%vh(1, iadd), pl_adds%vh(2, iadd), pl_adds%vh(3, iadd) else @@ -280,13 +399,5 @@ module subroutine symba_io_write_discard(self, param) return end subroutine symba_io_write_discard - - module subroutine symba_io_write_frame_info(self, iu, param) - implicit none - class(symba_particle_info), intent(in) :: self !! SyMBA particle info object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine symba_io_write_frame_info - end submodule s_symba_io diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90 index ab8b5543e..e06fb20b5 100644 --- a/src/symba/symba_setup.f90 +++ b/src/symba/symba_setup.f90 @@ -2,6 +2,62 @@ use swiftest contains + module subroutine symba_setup_initialize_particle_info(system, param) + !! author: David A. Minton + !! + !! Initializes a new particle information data structure with initial conditions recorded + implicit none + ! Argumets + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA extensions + ! Internals + integer(I4B) :: i + integer(I4B), dimension(:), allocatable :: idx + + select type(cb => system%cb) + class is (symba_cb) + cb%info%origin_type = "Central body" + cb%info%origin_time = param%t0 + cb%info%origin_xh(:) = 0.0_DP + cb%info%origin_vh(:) = 0.0_DP + call symba_io_dump_particle_info(system, param, lincludecb=.true.) + end select + + select type(pl => system%pl) + class is (symba_pl) + do i = 1, pl%nbody + pl%info(i)%origin_type = "Initial conditions" + pl%info(i)%origin_time = param%t0 + pl%info(i)%origin_xh(:) = pl%xh(:,i) + pl%info(i)%origin_vh(:) = pl%vh(:,i) + end do + if (pl%nbody > 0) then + allocate(idx(pl%nbody)) + call symba_io_dump_particle_info(system, param, plidx=[(i, i=1, pl%nbody)]) + deallocate(idx) + end if + end select + + select type(tp => system%tp) + class is (symba_tp) + do i = 1, tp%nbody + tp%info(i)%origin_type = "Initial conditions" + tp%info(i)%origin_time = param%t0 + tp%info(i)%origin_xh(:) = tp%xh(:,i) + tp%info(i)%origin_vh(:) = tp%vh(:,i) + end do + if (tp%nbody > 0) then + allocate(idx(tp%nbody)) + call symba_io_dump_particle_info(system, param, tpidx=[(i, i=1, tp%nbody)]) + deallocate(idx) + end if + end select + + + return + end subroutine symba_setup_initialize_particle_info + + module subroutine symba_setup_initialize_system(self, param) !! author: David A. Minton !! @@ -27,6 +83,11 @@ module subroutine symba_setup_initialize_system(self, param) class is (symba_parameters) pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY pl%nplm = count(pl%lmtiny(:)) + if (param%lrestart) then + call symba_io_read_particle(system, param) + else + call symba_setup_initialize_particle_info(system, param) + end if end select end select end associate @@ -162,10 +223,12 @@ module subroutine symba_setup_tp(self, n, param) if (allocated(self%nplenc)) deallocate(self%nplenc) if (allocated(self%levelg)) deallocate(self%levelg) if (allocated(self%levelm)) deallocate(self%levelm) + if (allocated(self%info)) deallocate(self%info) allocate(self%nplenc(n)) allocate(self%levelg(n)) allocate(self%levelm(n)) + allocate(self%info(n)) self%nplenc(:) = 0 self%levelg(:) = -1 diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 028b0678c..4c0f256e3 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -2,7 +2,7 @@ use swiftest contains - module subroutine symba_util_append_arr_info(arr, source, lsource_mask) + module subroutine symba_util_append_arr_info(arr, source, nold, nsrc, lsource_mask) !! author: David A. Minton !! !! Append a single array of particle information type onto another. If the destination array is not allocated, or is not big enough, this will allocate space for it. @@ -10,38 +10,24 @@ module subroutine symba_util_append_arr_info(arr, source, lsource_mask) ! Arguments type(symba_particle_info), dimension(:), allocatable, intent(inout) :: arr !! Destination array type(symba_particle_info), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to - ! Internals - integer(I4B) :: narr, nsrc + integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to if (.not. allocated(source)) return - if (present(lsource_mask)) then - nsrc = count(lsource_mask) + if (.not.allocated(arr)) then + allocate(arr(nold+nsrc)) else - nsrc = size(source) + call util_resize(arr, nold + nsrc) end if - if (allocated(arr)) then - narr = size(arr) - else - allocate(arr(nsrc)) - narr = 0 - end if - - call util_resize(arr, narr + nsrc) - - if (present(lsource_mask)) then - arr(narr + 1:narr + nsrc) = pack(source(:), lsource_mask(:)) - else - arr(narr + 1:narr + nsrc) = source(:) - end if + arr(nold + 1:nold + nsrc) = pack(source(1:nsrc), lsource_mask(1:nsrc)) return end subroutine symba_util_append_arr_info - module subroutine symba_util_append_arr_kin(arr, source, lsource_mask) + module subroutine symba_util_append_arr_kin(arr, source, nold, nsrc, lsource_mask) !! author: David A. Minton !! !! Append a single array of kinship type onto another. If the destination array is not allocated, or is not big enough, this will allocate space for it. @@ -49,32 +35,18 @@ module subroutine symba_util_append_arr_kin(arr, source, lsource_mask) ! Arguments type(symba_kinship), dimension(:), allocatable, intent(inout) :: arr !! Destination array type(symba_kinship), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to - ! Internals - integer(I4B) :: narr, nsrc + integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to if (.not. allocated(source)) return - if (present(lsource_mask)) then - nsrc = count(lsource_mask) + if (.not.allocated(arr)) then + allocate(arr(nold+nsrc)) else - nsrc = size(source) + call util_resize(arr, nold + nsrc) end if - if (allocated(arr)) then - narr = size(arr) - else - allocate(arr(nsrc)) - narr = 0 - end if - - call util_resize(arr, narr + nsrc) - - if (present(lsource_mask)) then - arr(narr + 1:narr + nsrc) = pack(source(:), lsource_mask(:)) - else - arr(narr + 1:narr + nsrc) = source(:) - end if + arr(nold + 1:nold + nsrc) = pack(source(1:nsrc), lsource_mask(1:nsrc)) return end subroutine symba_util_append_arr_kin @@ -89,24 +61,26 @@ module subroutine symba_util_append_pl(self, source, lsource_mask) !! Arguments class(symba_pl), intent(inout) :: self !! SyMBA massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to select type(source) class is (symba_pl) - call util_append_pl(self, source, lsource_mask) ! Note: helio_pl does not have its own append method, so we skip back to the base class - - call util_append(self%lcollision, source%lcollision, lsource_mask) - call util_append(self%lencounter, source%lencounter, lsource_mask) - call util_append(self%lmtiny, source%lmtiny, lsource_mask) - call util_append(self%nplenc, source%nplenc, lsource_mask) - call util_append(self%ntpenc, source%ntpenc, lsource_mask) - call util_append(self%levelg, source%levelg, lsource_mask) - call util_append(self%levelm, source%levelm, lsource_mask) - call util_append(self%isperi, source%isperi, lsource_mask) - call util_append(self%peri, source%peri, lsource_mask) - call util_append(self%atp, source%atp, lsource_mask) - call util_append(self%kin, source%kin, lsource_mask) - call util_append(self%info, source%info, lsource_mask) + associate(nold => self%nbody, nsrc => source%nbody) + call util_append(self%lcollision, source%lcollision, nold, nsrc, lsource_mask) + call util_append(self%lencounter, source%lencounter, nold, nsrc, lsource_mask) + call util_append(self%lmtiny, source%lmtiny, nold, nsrc, lsource_mask) + call util_append(self%nplenc, source%nplenc, nold, nsrc, lsource_mask) + call util_append(self%ntpenc, source%ntpenc, nold, nsrc, lsource_mask) + call util_append(self%levelg, source%levelg, nold, nsrc, lsource_mask) + call util_append(self%levelm, source%levelm, nold, nsrc, lsource_mask) + call util_append(self%isperi, source%isperi, nold, nsrc, lsource_mask) + call util_append(self%peri, source%peri, nold, nsrc, lsource_mask) + call util_append(self%atp, source%atp, nold, nsrc, lsource_mask) + call util_append(self%kin, source%kin, nold, nsrc, lsource_mask) + call util_append(self%info, source%info, nold, nsrc, lsource_mask) + + call util_append_pl(self, source, lsource_mask) ! Note: helio_pl does not have its own append method, so we skip back to the base class + end associate class default write(*,*) "Invalid object passed to the append method. Source must be of class symba_pl or its descendents!" call util_exit(FAILURE) @@ -125,24 +99,30 @@ module subroutine symba_util_append_merger(self, source, lsource_mask) ! Arguments class(symba_merger), intent(inout) :: self !! SyMBA massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to ! Internals integer(I4B), dimension(:), allocatable :: ncomp_tmp !! Temporary placeholder for ncomp incase we are appending a symba_pl object to a symba_merger + integer(I4B) :: nold, nsrc + nold = self%nbody + nsrc = source%nbody select type(source) class is (symba_merger) + call util_append(self%ncomp, source%ncomp, nold, nsrc, lsource_mask) call symba_util_append_pl(self, source, lsource_mask) - call util_append(self%ncomp, source%ncomp, lsource_mask) class is (symba_pl) - call symba_util_append_pl(self, source, lsource_mask) allocate(ncomp_tmp, mold=source%id) ncomp_tmp(:) = 0 - call util_append(self%ncomp, ncomp_tmp, lsource_mask) + call util_append(self%ncomp, ncomp_tmp, nold, nsrc, lsource_mask) + call symba_util_append_pl(self, source, lsource_mask) class default write(*,*) "Invalid object passed to the append method. Source must be of class symba_pl or its descendents!" call util_exit(FAILURE) end select + ! Save the number of appended bodies + self%ncomp(nold+1:nold+nsrc) = nsrc + return end subroutine symba_util_append_merger @@ -156,15 +136,18 @@ module subroutine symba_util_append_tp(self, source, lsource_mask) !! Arguments class(symba_tp), intent(inout) :: self !! SyMBA test particle object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to select type(source) class is (symba_tp) - call util_append_tp(self, source, lsource_mask) ! Note: helio_tp does not have its own append method, so we skip back to the base class - - call util_append(self%nplenc, source%nplenc, lsource_mask) - call util_append(self%levelg, source%levelg, lsource_mask) - call util_append(self%levelm, source%levelm, lsource_mask) + associate(nold => self%nbody, nsrc => source%nbody) + call util_append(self%nplenc, source%nplenc, nold, nsrc, lsource_mask) + call util_append(self%levelg, source%levelg, nold, nsrc, lsource_mask) + call util_append(self%levelm, source%levelm, nold, nsrc, lsource_mask) + call util_append(self%info, source%info, nold, nsrc, lsource_mask) + + call util_append_tp(self, source, lsource_mask) ! Note: helio_tp does not have its own append method, so we skip back to the base class + end associate class default write(*,*) "Invalid object passed to the append method. Source must be of class symba_tp or its descendents!" call util_exit(FAILURE) @@ -271,6 +254,7 @@ module subroutine symba_util_fill_tp(self, inserts, lfill_list) call util_fill(keeps%nplenc, inserts%nplenc, lfill_list) call util_fill(keeps%levelg, inserts%levelg, lfill_list) call util_fill(keeps%levelm, inserts%levelm, lfill_list) + call util_fill(keeps%info, inserts%info, lfill_list) call util_fill_tp(keeps, inserts, lfill_list) ! Note: helio_tp does not have its own fill method, so we skip back to the base class class default @@ -380,11 +364,15 @@ module subroutine symba_util_rearray_pl(self, system, param) class(symba_parameters), intent(in) :: param !! Current run configuration parameters ! Internals class(symba_pl), allocatable :: tmp !! The discarded body list. + integer(I4B) :: i + logical, dimension(:), allocatable :: lmask associate(pl => self, pl_adds => system%pl_adds) allocate(tmp, mold=pl) ! Remove the discards and destroy the list, as the system already tracks pl_discards elsewhere - call pl%spill(tmp, lspill_list=(pl%ldiscard(:) .or. pl%status(:) == INACTIVE), ldestructive=.true.) + allocate(lmask, source=pl%ldiscard(:)) + lmask(:) = lmask(:) .or. pl%status(:) == INACTIVE + call pl%spill(tmp, lspill_list=lmask, ldestructive=.true.) call tmp%setup(0,param) deallocate(tmp) @@ -393,7 +381,10 @@ module subroutine symba_util_rearray_pl(self, system, param) if (allocated(pl%xend)) deallocate(pl%xend) ! Add in any new bodies - call pl%append(pl_adds) + if (pl_adds%nbody > 0) then + call pl%append(pl_adds, lsource_mask=[(.true., i=1, pl_adds%nbody)]) + call symba_io_dump_particle_info(system, param, plidx=[(i, i = 1, pl%nbody)]) + end if ! If there are still bodies in the system, sort by mass in descending order and re-index if (pl%nbody > 0) then @@ -486,10 +477,10 @@ module subroutine symba_util_resize_merger(self, nnew) class(symba_merger), intent(inout) :: self !! SyMBA massive body object integer(I4B), intent(in) :: nnew !! New size neded - call symba_util_resize_pl(self, nnew) - call util_resize(self%ncomp, nnew) + call symba_util_resize_pl(self, nnew) + return end subroutine symba_util_resize_merger @@ -503,8 +494,6 @@ module subroutine symba_util_resize_pl(self, nnew) class(symba_pl), intent(inout) :: self !! SyMBA massive body object integer(I4B), intent(in) :: nnew !! New size neded - call util_resize_pl(self, nnew) - call util_resize(self%lcollision, nnew) call util_resize(self%lencounter, nnew) call util_resize(self%lmtiny, nnew) @@ -518,6 +507,8 @@ module subroutine symba_util_resize_pl(self, nnew) call util_resize(self%kin, nnew) call util_resize(self%info, nnew) + call util_resize_pl(self, nnew) + return end subroutine symba_util_resize_pl @@ -531,11 +522,12 @@ module subroutine symba_util_resize_tp(self, nnew) class(symba_tp), intent(inout) :: self !! SyMBA test particle object integer(I4B), intent(in) :: nnew !! New size neded - call util_resize_tp(self, nnew) - call util_resize(self%nplenc, nnew) call util_resize(self%levelg, nnew) call util_resize(self%levelm, nnew) + call util_resize(self%info, nnew) + + call util_resize_tp(self, nnew) return end subroutine symba_util_resize_tp @@ -693,6 +685,7 @@ module subroutine symba_util_sort_rearrange_tp(self, ind) if (allocated(tp%nplenc)) tp%nplenc(1:ntp) = tp_sorted%nplenc(ind(1:ntp)) if (allocated(tp%levelg)) tp%levelg(1:ntp) = tp_sorted%levelg(ind(1:ntp)) if (allocated(tp%levelm)) tp%levelm(1:ntp) = tp_sorted%levelm(ind(1:ntp)) + if (allocated(tp%info)) tp%info(1:ntp) = tp_sorted%info(ind(1:ntp)) deallocate(tp_sorted) end associate @@ -850,6 +843,7 @@ module subroutine symba_util_spill_tp(self, discards, lspill_list, ldestructive) call util_spill(keeps%nplenc, discards%nplenc, lspill_list, ldestructive) call util_spill(keeps%levelg, discards%levelg, lspill_list, ldestructive) call util_spill(keeps%levelm, discards%levelm, lspill_list, ldestructive) + call util_spill(keeps%info, discards%info, lspill_list, ldestructive) call util_spill_tp(keeps, discards, lspill_list, ldestructive) class default diff --git a/src/util/util_append.f90 b/src/util/util_append.f90 index 0f7ac0bde..dc48f9861 100644 --- a/src/util/util_append.f90 +++ b/src/util/util_append.f90 @@ -2,7 +2,7 @@ use swiftest contains - module subroutine util_append_arr_char_string(arr, source, lsource_mask) + module subroutine util_append_arr_char_string(arr, source, nold, nsrc, lsource_mask) !! author: David A. Minton !! !! Append a single array of character string type onto another. If the destination array is not allocated, or is not big enough, this will allocate space for it. @@ -10,38 +10,24 @@ module subroutine util_append_arr_char_string(arr, source, lsource_mask) ! Arguments character(len=STRMAX), dimension(:), allocatable, intent(inout) :: arr !! Destination array character(len=STRMAX), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to - ! Internals - integer(I4B) :: narr, nsrc + integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to if (.not. allocated(source)) return - if (present(lsource_mask)) then - nsrc = count(lsource_mask) + if (.not.allocated(arr)) then + allocate(arr(nold+nsrc)) else - nsrc = size(source) + call util_resize(arr, nold + nsrc) end if - if (allocated(arr)) then - narr = size(arr) - else - allocate(arr(nsrc)) - narr = 0 - end if - - call util_resize(arr, narr + nsrc) - - if (present(lsource_mask)) then - arr(narr + 1:narr + nsrc) = pack(source(:), lsource_mask(:)) - else - arr(narr + 1:narr + nsrc) = source(:) - end if + arr(nold + 1:nold + nsrc) = pack(source(1:nsrc), lsource_mask(1:nsrc)) return end subroutine util_append_arr_char_string - module subroutine util_append_arr_DP(arr, source, lsource_mask) + module subroutine util_append_arr_DP(arr, source, nold, nsrc, lsource_mask) !! author: David A. Minton !! !! Append a single array of double precision type onto another. If the destination array is not allocated, or is not big enough, this will allocate space for it. @@ -49,38 +35,24 @@ module subroutine util_append_arr_DP(arr, source, lsource_mask) ! Arguments real(DP), dimension(:), allocatable, intent(inout) :: arr !! Destination array real(DP), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to - ! Internals - integer(I4B) :: narr, nsrc + integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to if (.not. allocated(source)) return - if (present(lsource_mask)) then - nsrc = count(lsource_mask) + if (.not.allocated(arr)) then + allocate(arr(nold+nsrc)) else - nsrc = size(source) + call util_resize(arr, nold + nsrc) end if - if (allocated(arr)) then - narr = size(arr) - else - allocate(arr(nsrc)) - narr = 0 - end if - - call util_resize(arr, narr + nsrc) - - if (present(lsource_mask)) then - arr(narr + 1:narr + nsrc) = pack(source(:), lsource_mask(:)) - else - arr(narr + 1:narr + nsrc) = source(:) - end if + arr(nold + 1:nold + nsrc) = pack(source(1:nsrc), lsource_mask(1:nsrc)) return end subroutine util_append_arr_DP - module subroutine util_append_arr_DPvec(arr, source, lsource_mask) + module subroutine util_append_arr_DPvec(arr, source, nold, nsrc, lsource_mask) !! author: David A. Minton !! !! Append a single array of double precision vector type of size (NDIM, n) onto another. If the destination array is not allocated, or is not big enough, this will allocate space for it. @@ -88,40 +60,26 @@ module subroutine util_append_arr_DPvec(arr, source, lsource_mask) ! Arguments real(DP), dimension(:,:), allocatable, intent(inout) :: arr !! Destination array real(DP), dimension(:,:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to - ! Internals - integer(I4B) :: narr, nsrc + integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to if (.not. allocated(source)) return - if (present(lsource_mask)) then - nsrc = count(lsource_mask) - else - nsrc = size(source, dim=2) - end if - - if (allocated(arr)) then - narr = size(arr, dim=2) + if (.not.allocated(arr)) then + allocate(arr(NDIM, nold+nsrc)) else - allocate(arr(NDIM, nsrc)) - narr = 0 + call util_resize(arr, nold + nsrc) end if - call util_resize(arr, narr + nsrc) - - if (present(lsource_mask)) then - arr(1, narr + 1:narr + nsrc) = pack(source(1,:), lsource_mask(:)) - arr(2, narr + 1:narr + nsrc) = pack(source(2,:), lsource_mask(:)) - arr(3, narr + 1:narr + nsrc) = pack(source(3,:), lsource_mask(:)) - else - arr(:, narr + 1:narr + nsrc) = source(:,:) - end if + arr(1, nold + 1:nold + nsrc) = pack(source(1,1:nsrc), lsource_mask(1:nsrc)) + arr(2, nold + 1:nold + nsrc) = pack(source(2,1:nsrc), lsource_mask(1:nsrc)) + arr(3, nold + 1:nold + nsrc) = pack(source(3,1:nsrc), lsource_mask(1:nsrc)) return end subroutine util_append_arr_DPvec - module subroutine util_append_arr_I4B(arr, source, lsource_mask) + module subroutine util_append_arr_I4B(arr, source, nold, nsrc, lsource_mask) !! author: David A. Minton !! !! Append a single array of integer(I4B) onto another. If the destination array is not allocated, or is not big enough, this will allocate space for it. @@ -129,38 +87,24 @@ module subroutine util_append_arr_I4B(arr, source, lsource_mask) ! Arguments integer(I4B), dimension(:), allocatable, intent(inout) :: arr !! Destination array integer(I4B), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to - ! Internals - integer(I4B) :: narr, nsrc + integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to if (.not. allocated(source)) return - if (present(lsource_mask)) then - nsrc = count(lsource_mask) - else - nsrc = size(source) - end if - - if (allocated(arr)) then - narr = size(arr) + if (.not.allocated(arr)) then + allocate(arr(nold+nsrc)) else - allocate(arr(nsrc)) - narr = 0 + call util_resize(arr, nold + nsrc) end if - call util_resize(arr, narr + nsrc) - - if (present(lsource_mask)) then - arr(narr + 1:narr + nsrc) = pack(source(:), lsource_mask(:)) - else - arr(narr + 1:narr + nsrc) = source(:) - end if + arr(nold + 1:nold + nsrc) = pack(source(1:nsrc), lsource_mask(1:nsrc)) return end subroutine util_append_arr_I4B - module subroutine util_append_arr_logical(arr, source, lsource_mask) + module subroutine util_append_arr_logical(arr, source, nold, nsrc, lsource_mask) !! author: David A. Minton !! !! Append a single array of logical type onto another. If the destination array is not allocated, or is not big enough, this will allocate space for it. @@ -168,32 +112,18 @@ module subroutine util_append_arr_logical(arr, source, lsource_mask) ! Arguments logical, dimension(:), allocatable, intent(inout) :: arr !! Destination array logical, dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to - ! Internals - integer(I4B) :: narr, nsrc + integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to if (.not. allocated(source)) return - if (allocated(arr)) then - narr = size(arr) - else - allocate(arr(nsrc)) - narr = 0 - end if - - if (present(lsource_mask)) then - nsrc = count(lsource_mask) + if (.not.allocated(arr)) then + allocate(arr(nold+nsrc)) else - nsrc = size(source) + call util_resize(arr, nold + nsrc) end if - call util_resize(arr, narr + nsrc) - - if (present(lsource_mask)) then - arr(narr + 1:narr + nsrc) = pack(source(:), lsource_mask(:)) - else - arr(narr + 1:narr + nsrc) = source(:) - end if + arr(nold + 1:nold + nsrc) = pack(source(1:nsrc), lsource_mask(1:nsrc)) return end subroutine util_append_arr_logical @@ -208,29 +138,31 @@ module subroutine util_append_body(self, source, lsource_mask) ! Arguments class(swiftest_body), intent(inout) :: self !! Swiftest body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to - - call util_append(self%name, source%name, lsource_mask) - call util_append(self%id, source%id, lsource_mask) - call util_append(self%status, source%status, lsource_mask) - call util_append(self%ldiscard, source%ldiscard, lsource_mask) - call util_append(self%lmask, source%lmask, lsource_mask) - call util_append(self%mu, source%mu, lsource_mask) - call util_append(self%xh, source%xh, lsource_mask) - call util_append(self%vh, source%vh, lsource_mask) - call util_append(self%xb, source%xb, lsource_mask) - call util_append(self%vb, source%vb, lsource_mask) - call util_append(self%ah, source%ah, lsource_mask) - call util_append(self%aobl, source%aobl, lsource_mask) - call util_append(self%atide, source%atide, lsource_mask) - call util_append(self%agr, source%agr, lsource_mask) - call util_append(self%ir3h, source%ir3h, lsource_mask) - call util_append(self%a, source%a, lsource_mask) - call util_append(self%e, source%e, lsource_mask) - call util_append(self%inc, source%inc, lsource_mask) - call util_append(self%capom, source%capom, lsource_mask) - call util_append(self%omega, source%omega, lsource_mask) - call util_append(self%capm, source%capm, lsource_mask) + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + + associate(nold => self%nbody, nsrc => source%nbody) + call util_append(self%name, source%name, nold, nsrc, lsource_mask) + call util_append(self%id, source%id, nold, nsrc, lsource_mask) + call util_append(self%status, source%status, nold, nsrc, lsource_mask) + call util_append(self%ldiscard, source%ldiscard, nold, nsrc, lsource_mask) + call util_append(self%lmask, source%lmask, nold, nsrc, lsource_mask) + call util_append(self%mu, source%mu, nold, nsrc, lsource_mask) + call util_append(self%xh, source%xh, nold, nsrc, lsource_mask) + call util_append(self%vh, source%vh, nold, nsrc, lsource_mask) + call util_append(self%xb, source%xb, nold, nsrc, lsource_mask) + call util_append(self%vb, source%vb, nold, nsrc, lsource_mask) + call util_append(self%ah, source%ah, nold, nsrc, lsource_mask) + call util_append(self%aobl, source%aobl, nold, nsrc, lsource_mask) + call util_append(self%atide, source%atide, nold, nsrc, lsource_mask) + call util_append(self%agr, source%agr, nold, nsrc, lsource_mask) + call util_append(self%ir3h, source%ir3h, nold, nsrc, lsource_mask) + call util_append(self%a, source%a, nold, nsrc, lsource_mask) + call util_append(self%e, source%e, nold, nsrc, lsource_mask) + call util_append(self%inc, source%inc, nold, nsrc, lsource_mask) + call util_append(self%capom, source%capom, nold, nsrc, lsource_mask) + call util_append(self%omega, source%omega, nold, nsrc, lsource_mask) + call util_append(self%capm, source%capm, nold, nsrc, lsource_mask) + end associate self%nbody = count(self%status(:) /= INACTIVE) @@ -247,26 +179,28 @@ module subroutine util_append_pl(self, source, lsource_mask) ! Arguments class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to select type(source) class is (swiftest_pl) - call util_append_body(self, source, lsource_mask) - - call util_append(self%mass, source%mass, lsource_mask) - call util_append(self%Gmass, source%Gmass, lsource_mask) - call util_append(self%rhill, source%rhill, lsource_mask) - call util_append(self%radius, source%radius, lsource_mask) - call util_append(self%xbeg, source%xbeg, lsource_mask) - call util_append(self%xend, source%xend, lsource_mask) - call util_append(self%vbeg, source%vbeg, lsource_mask) - call util_append(self%density, source%density, lsource_mask) - call util_append(self%Ip, source%Ip, lsource_mask) - call util_append(self%rot, source%rot, lsource_mask) - call util_append(self%k2, source%k2, lsource_mask) - call util_append(self%Q, source%Q, lsource_mask) - call util_append(self%tlag, source%tlag, lsource_mask) + associate(nold => self%nbody, nsrc => source%nbody) + call util_append(self%mass, source%mass, nold, nsrc, lsource_mask) + call util_append(self%Gmass, source%Gmass, nold, nsrc, lsource_mask) + call util_append(self%rhill, source%rhill, nold, nsrc, lsource_mask) + call util_append(self%radius, source%radius, nold, nsrc, lsource_mask) + call util_append(self%xbeg, source%xbeg, nold, nsrc, lsource_mask) + call util_append(self%xend, source%xend, nold, nsrc, lsource_mask) + call util_append(self%vbeg, source%vbeg, nold, nsrc, lsource_mask) + call util_append(self%density, source%density, nold, nsrc, lsource_mask) + call util_append(self%Ip, source%Ip, nold, nsrc, lsource_mask) + call util_append(self%rot, source%rot, nold, nsrc, lsource_mask) + call util_append(self%k2, source%k2, nold, nsrc, lsource_mask) + call util_append(self%Q, source%Q, nold, nsrc, lsource_mask) + call util_append(self%tlag, source%tlag, nold, nsrc, lsource_mask) + + call util_append_body(self, source, lsource_mask) + end associate call self%eucl_index() class default @@ -287,15 +221,17 @@ module subroutine util_append_tp(self, source, lsource_mask) ! Arguments class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to select type(source) class is (swiftest_tp) - call util_append_body(self, source, lsource_mask) + associate(nold => self%nbody, nsrc => source%nbody) + call util_append(self%isperi, source%isperi, nold, nsrc, lsource_mask) + call util_append(self%peri, source%peri, nold, nsrc, lsource_mask) + call util_append(self%atp, source%atp, nold, nsrc, lsource_mask) - call util_append(self%isperi, source%isperi, lsource_mask) - call util_append(self%peri, source%peri, lsource_mask) - call util_append(self%atp, source%atp, lsource_mask) + call util_append_body(self, source, lsource_mask) + end associate class default write(*,*) "Invalid object passed to the append method. Source must be of class swiftest_tp or its descendents" call util_exit(FAILURE) diff --git a/src/util/util_coord.f90 b/src/util/util_coord.f90 index c10dbace7..2a970d0dc 100644 --- a/src/util/util_coord.f90 +++ b/src/util/util_coord.f90 @@ -24,6 +24,7 @@ module subroutine util_coord_h2b_pl(self, cb) xtmp(:) = 0.0_DP vtmp(:) = 0.0_DP do i = 1, npl + if (pl%status(i) == INACTIVE) cycle Gmtot = Gmtot + pl%Gmass(i) xtmp(:) = xtmp(:) + pl%Gmass(i) * pl%xh(:,i) vtmp(:) = vtmp(:) + pl%Gmass(i) * pl%vh(:,i) @@ -31,6 +32,7 @@ module subroutine util_coord_h2b_pl(self, cb) cb%xb(:) = -xtmp(:) / Gmtot cb%vb(:) = -vtmp(:) / Gmtot do i = 1, npl + if (pl%status(i) == INACTIVE) cycle pl%xb(:,i) = pl%xh(:,i) + cb%xb(:) pl%vb(:,i) = pl%vh(:,i) + cb%vb(:) end do @@ -51,20 +53,15 @@ module subroutine util_coord_h2b_tp(self, cb) ! Arguments class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object class(swiftest_cb), intent(in) :: cb !! Swiftest central body object + ! Internals + integer(I4B) :: i if (self%nbody == 0) return - associate(ntp => self%nbody, xbcb => cb%xb, vbcb => cb%vb, status => self%status, & - xb => self%xb, xh => self%xh, vb => self%vb, vh => self%vh) - - where(status(1:ntp) /= INACTIVE) - xb(1, 1:ntp) = xh(1, 1:ntp) + xbcb(1) - xb(2, 1:ntp) = xh(2, 1:ntp) + xbcb(2) - xb(3, 1:ntp) = xh(3, 1:ntp) + xbcb(3) - - vb(1, 1:ntp) = vh(1, 1:ntp) + vbcb(1) - vb(2, 1:ntp) = vh(2, 1:ntp) + vbcb(2) - vb(3, 1:ntp) = vh(3, 1:ntp) + vbcb(3) - end where + associate(tp => self, ntp => self%nbody) + do concurrent (i = 1:ntp, tp%status(i) /= INACTIVE) + tp%xb(:, i) = tp%xh(:, i) + cb%xb(:) + tp%vb(:, i) = tp%vh(:, i) + cb%vb(:) + end do end associate return @@ -87,11 +84,10 @@ module subroutine util_coord_b2h_pl(self, cb) if (self%nbody == 0) return - associate(npl => self%nbody, xbcb => cb%xb, vbcb => cb%vb, xb => self%xb, xh => self%xh, & - vb => self%vb, vh => self%vh) - do i = 1, NDIM - xh(i, 1:npl) = xb(i, 1:npl) - xbcb(i) - vh(i, 1:npl) = vb(i, 1:npl) - vbcb(i) + associate(pl => self, npl => self%nbody) + do concurrent (i = 1:npl, pl%status(i) /= INACTIVE) + pl%xh(:, i) = pl%xb(:, i) - cb%xb(:) + pl%vh(:, i) = pl%vb(:, i) - cb%vb(:) end do end associate @@ -110,20 +106,16 @@ module subroutine util_coord_b2h_tp(self, cb) ! Arguments class(swiftest_tp), intent(inout) :: self !! Swiftest massive body object class(swiftest_cb), intent(in) :: cb !! Swiftest central body object + ! Internals + integer(I4B) :: i if (self%nbody == 0) return - associate(ntp => self%nbody, xbcb => cb%xb, vbcb => cb%vb, xb => self%xb, xh => self%xh, & - vb => self%vb, vh => self%vh, status => self%status) - where(status(1:ntp) /= INACTIVE) - xh(1, 1:ntp) = xb(1, 1:ntp) - xbcb(1) - xh(2, 1:ntp) = xb(2, 1:ntp) - xbcb(2) - xh(3, 1:ntp) = xb(3, 1:ntp) - xbcb(3) - - vh(1, 1:ntp) = vb(1, 1:ntp) - vbcb(1) - vh(2, 1:ntp) = vb(2, 1:ntp) - vbcb(2) - vh(3, 1:ntp) = vb(3, 1:ntp) - vbcb(3) - end where + associate(tp => self, ntp => self%nbody) + do concurrent(i = 1:ntp, tp%status(i) /= INACTIVE) + tp%xh(:, i) = tp%xb(:, i) - cb%xb(:) + tp%vh(:, i) = tp%vb(:, i) - cb%vb(:) + end do end associate return diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index 90f0d2242..fa7cda43d 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -15,12 +15,12 @@ module subroutine util_get_energy_momentum_system(self, param) ! Internals integer(I4B) :: i, j integer(I8B) :: k - real(DP) :: rmag, v2, rot2, oblpot, hx, hy, hz, hsx, hsy, hsz - real(DP) :: kecb, kespincb, Lcborbitx, Lcborbity, Lcborbitz, Lcbspinx, Lcbspiny, Lcbspinz + real(DP) :: oblpot, kecb, kespincb real(DP), dimension(self%pl%nbody) :: irh, kepl, kespinpl, pecb real(DP), dimension(self%pl%nbody) :: Lplorbitx, Lplorbity, Lplorbitz real(DP), dimension(self%pl%nbody) :: Lplspinx, Lplspiny, Lplspinz real(DP), dimension(self%pl%nplpl) :: pepl + real(DP), dimension(NDIM) :: Lcborbit, Lcbspin logical, dimension(self%pl%nplpl) :: lstatpl logical, dimension(self%pl%nbody) :: lstatus @@ -38,54 +38,50 @@ module subroutine util_get_energy_momentum_system(self, param) Lplspiny(:) = 0.0_DP Lplspinz(:) = 0.0_DP lstatus(1:npl) = pl%status(1:npl) /= INACTIVE - call pl%h2b(cb) + kecb = cb%mass * dot_product(cb%vb(:), cb%vb(:)) - hx = cb%xb(2) * cb%vb(3) - cb%xb(3) * cb%vb(2) - hy = cb%xb(3) * cb%vb(1) - cb%xb(1) * cb%vb(3) - hz = cb%xb(1) * cb%vb(2) - cb%xb(2) * cb%vb(1) - Lcborbitx = cb%mass * hx - Lcborbity = cb%mass * hy - Lcborbitz = cb%mass * hz - !!$omp simd private(v2, rot2, hx, hy, hz) - do i = 1, npl - v2 = dot_product(pl%vb(:,i), pl%vb(:,i)) - hx = pl%xb(2,i) * pl%vb(3,i) - pl%xb(3,i) * pl%vb(2,i) - hy = pl%xb(3,i) * pl%vb(1,i) - pl%xb(1,i) * pl%vb(3,i) - hz = pl%xb(1,i) * pl%vb(2,i) - pl%xb(2,i) * pl%vb(1,i) + Lcborbit(:) = cb%mass * cb%xb(:) .cross. cb%vb(:) + + do concurrent (i = 1:npl, lstatus(i)) + block ! We use a block construct to prevent generating temporary arrays for local variables + real(DP) :: v2, hx, hy, hz + v2 = dot_product(pl%vb(:,i), pl%vb(:,i)) + hx = pl%xb(2,i) * pl%vb(3,i) - pl%xb(3,i) * pl%vb(2,i) + hy = pl%xb(3,i) * pl%vb(1,i) - pl%xb(1,i) * pl%vb(3,i) + hz = pl%xb(1,i) * pl%vb(2,i) - pl%xb(2,i) * pl%vb(1,i) - ! Angular momentum from orbit - Lplorbitx(i) = pl%mass(i) * hx - Lplorbity(i) = pl%mass(i) * hy - Lplorbitz(i) = pl%mass(i) * hz + ! Angular momentum from orbit + Lplorbitx(i) = pl%mass(i) * hx + Lplorbity(i) = pl%mass(i) * hy + Lplorbitz(i) = pl%mass(i) * hz - ! Kinetic energy from orbit and spin - kepl(i) = pl%mass(i) * v2 + ! Kinetic energy from orbit and spin + kepl(i) = pl%mass(i) * v2 + end block end do if (param%lrotation) then kespincb = cb%mass * cb%Ip(3) * cb%radius**2 * dot_product(cb%rot(:), cb%rot(:)) - ! For simplicity, we always assume that the rotation pole is the 3rd principal axis - hsx = cb%Ip(3) * cb%radius**2 * cb%rot(1) - hsy = cb%Ip(3) * cb%radius**2 * cb%rot(2) - hsz = cb%Ip(3) * cb%radius**2 * cb%rot(3) - ! Angular momentum from spin - Lcbspinx = cb%mass * hsx - Lcbspiny = cb%mass * hsy - Lcbspinz = cb%mass * hsz + ! For simplicity, we always assume that the rotation pole is the 3rd principal axis + Lcbspin(:) = cb%Ip(3) * cb%mass * cb%radius**2 * cb%rot(:) - do i = 1, npl - rot2 = dot_product(pl%rot(:,i), pl%rot(:,i)) - ! For simplicity, we always assume that the rotation pole is the 3rd principal axis - hsx = pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(1,i) - hsy = pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(2,i) - hsz = pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(3,i) + do concurrent (i = 1:npl, lstatus(i)) + block + real(DP) :: rot2, hsx, hsy, hsz - ! Angular momentum from spin - Lplspinx(i) = pl%mass(i) * hsx - Lplspiny(i) = pl%mass(i) * hsy - Lplspinz(i) = pl%mass(i) * hsz - kespinpl(i) = pl%mass(i) * pl%Ip(3, i) * pl%radius(i)**2 * rot2 + rot2 = dot_product(pl%rot(:,i), pl%rot(:,i)) + ! For simplicity, we always assume that the rotation pole is the 3rd principal axis + hsx = pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(1,i) + hsy = pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(2,i) + hsz = pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(3,i) + + ! Angular momentum from spin + Lplspinx(i) = pl%mass(i) * hsx + Lplspiny(i) = pl%mass(i) * hsy + Lplspinz(i) = pl%mass(i) * hsz + kespinpl(i) = pl%mass(i) * pl%Ip(3, i) * pl%radius(i)**2 * rot2 + end block end do else kespincb = 0.0_DP @@ -93,45 +89,48 @@ module subroutine util_get_energy_momentum_system(self, param) end if ! Do the central body potential energy component first - !$omp simd - do i = 1, npl - associate(px => pl%xb(1,i), py => pl%xb(2,i), pz => pl%xb(3,i)) - pecb(i) = -cb%Gmass * pl%mass(i) / sqrt(px**2 + py**2 + pz**2) - end associate - end do + associate(px => pl%xb(1,:), py => pl%xb(2,:), pz => pl%xb(3,:)) + do concurrent(i = 1:npl, lstatus(i)) + pecb(i) = -cb%Gmass * pl%mass(i) / sqrt(px(i)**2 + py(i)**2 + pz(i)**2) + end do + end associate ! Do the potential energy between pairs of massive bodies - do k = 1, pl%nplpl - associate(ik => pl%k_plpl(1, k), jk => pl%k_plpl(2, k)) - pepl(k) = -pl%Gmass(ik) * pl%mass(jk) / norm2(pl%xb(:, jk) - pl%xb(:, ik)) - lstatpl(k) = (lstatus(ik) .and. lstatus(jk)) - end associate - end do + associate(indi => pl%k_plpl(1, :), indj => pl%k_plpl(2, :)) + do concurrent (k = 1:pl%nplpl) + lstatpl(k) = (lstatus(indi(k)) .and. lstatus(indj(k))) + end do + + do concurrent (k = 1:pl%nplpl, lstatpl(k)) + pepl(k) = -pl%Gmass(indi(k)) * pl%mass(indj(k)) / norm2(pl%xb(:, indi(k)) - pl%xb(:, indj(k))) + end do + end associate + system%pe = sum(pepl(:), lstatpl(:)) + sum(pecb(1:npl), lstatus(1:npl)) + system%ke_orbit = 0.5_DP * (kecb + sum(kepl(1:npl), lstatus(:))) if (param%lrotation) system%ke_spin = 0.5_DP * (kespincb + sum(kespinpl(1:npl), lstatus(:))) - system%pe = sum(pepl(:), lstatpl(:)) + sum(pecb(1:npl), lstatus(1:npl)) - ! Potential energy from the oblateness term if (param%loblatecb) then - !$omp simd - do i = 1, npl + do concurrent(i = 1:npl, lstatus(i)) irh(i) = 1.0_DP / norm2(pl%xh(:,i)) end do call obl_pot(npl, cb%Gmass, pl%mass, cb%j2rp2, cb%j4rp4, pl%xh, irh, oblpot) system%pe = system%pe + oblpot end if - system%Lorbit(1) = Lcborbitx + sum(Lplorbitx(1:npl), lstatus(1:npl)) - system%Lorbit(2) = Lcborbity + sum(Lplorbity(1:npl), lstatus(1:npl)) - system%Lorbit(3) = Lcborbitz + sum(Lplorbitz(1:npl), lstatus(1:npl)) + system%Lorbit(1) = Lcborbit(1) + sum(Lplorbitx(1:npl), lstatus(1:npl)) + system%Lorbit(2) = Lcborbit(2) + sum(Lplorbity(1:npl), lstatus(1:npl)) + system%Lorbit(3) = Lcborbit(3) + sum(Lplorbitz(1:npl), lstatus(1:npl)) if (param%lrotation) then - system%Lspin(1) = Lcbspinx + sum(Lplspinx(1:npl), lstatus(1:npl)) - system%Lspin(2) = Lcbspiny + sum(Lplspiny(1:npl), lstatus(1:npl)) - system%Lspin(3) = Lcbspinz + sum(Lplspinz(1:npl), lstatus(1:npl)) + system%Lspin(1) = Lcbspin(1) + sum(Lplspinx(1:npl), lstatus(1:npl)) + system%Lspin(2) = Lcbspin(2) + sum(Lplspiny(1:npl), lstatus(1:npl)) + system%Lspin(3) = Lcbspin(3) + sum(Lplspinz(1:npl), lstatus(1:npl)) end if + + system%te = system%ke_orbit + system%ke_spin + system%pe end associate return diff --git a/src/util/util_rescale.f90 b/src/util/util_rescale.f90 index 061ecf9a5..62a9409ec 100644 --- a/src/util/util_rescale.f90 +++ b/src/util/util_rescale.f90 @@ -33,6 +33,7 @@ module subroutine util_rescale_system(self, param, mscale, dscale, tscale) cb%radius = cb%radius / dscale cb%xb(:) = cb%xb(:) / dscale cb%vb(:) = cb%vb(:) / vscale + cb%rot(:) = cb%rot(:) * tscale pl%mass(1:npl) = pl%mass(1:npl) / mscale pl%Gmass(1:npl) = param%GU * pl%mass(1:npl) pl%radius(1:npl) = pl%radius(1:npl) / dscale diff --git a/src/whm/whm_util.f90 b/src/whm/whm_util.f90 index f3dc15d3e..cc84ba3d5 100644 --- a/src/whm/whm_util.f90 +++ b/src/whm/whm_util.f90 @@ -11,17 +11,19 @@ module subroutine whm_util_append_pl(self, source, lsource_mask) !! Arguments class(whm_pl), intent(inout) :: self !! WHM massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to select type(source) class is (whm_pl) - call util_append_pl(self, source, lsource_mask) - - call util_append(self%eta, source%eta, lsource_mask) - call util_append(self%muj, source%muj, lsource_mask) - call util_append(self%ir3j, source%ir3j, lsource_mask) - call util_append(self%xj, source%xj, lsource_mask) - call util_append(self%vj, source%vj, lsource_mask) + associate(nold => self%nbody, nsrc => source%nbody) + call util_append(self%eta, source%eta, nold, nsrc, lsource_mask) + call util_append(self%muj, source%muj, nold, nsrc, lsource_mask) + call util_append(self%ir3j, source%ir3j, nold, nsrc, lsource_mask) + call util_append(self%xj, source%xj, nold, nsrc, lsource_mask) + call util_append(self%vj, source%vj, nold, nsrc, lsource_mask) + + call util_append_pl(self, source, lsource_mask) + end associate class default write(*,*) "Invalid object passed to the append method. Source must be of class whm_pl or its descendents" call util_exit(FAILURE) @@ -74,14 +76,14 @@ module subroutine whm_util_resize_pl(self, nnew) class(whm_pl), intent(inout) :: self !! WHM massive body object integer(I4B), intent(in) :: nnew !! New size neded - call util_resize_pl(self, nnew) - call util_resize(self%eta, nnew) call util_resize(self%xj, nnew) call util_resize(self%vj, nnew) call util_resize(self%muj, nnew) call util_resize(self%ir3j, nnew) + call util_resize_pl(self, nnew) + return end subroutine whm_util_resize_pl