Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Made substantial improvements to the convergence of Fraggle by isolat…
Browse files Browse the repository at this point in the history 
…ing the kinetic energy and momentum of the collision system. Prior to this, spin kinetic energy of the Sun dominated the calculation of the test cases.
  • Loading branch information
@daminton
daminton committed Jan 10, 2023
1 parent 6947240 commit 0362420
Show file tree
Hide file tree
Showing 5 changed files with 206 additions and 155 deletions.
77 changes: 50 additions & 27 deletions examples/Fragmentation/Fragmentation_Movie.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,13 +31,12 @@
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from pathlib import Path

# ----------------------------------------------------------------------------------------------------------------------
# Define the names and initial conditions of the various fragmentation simulation types
# ----------------------------------------------------------------------------------------------------------------------
available_movie_styles = ["disruption_headon", "disruption_off_axis", "supercatastrophic_headon", "supercatastrophic_off_axis","hitandrun_disrupt", "hitandrun_pure"]
movie_title_list = ["Head-on Disruption", "Off-axis Disruption", "Head-on Supercatastrophic", "Off-axis Supercatastrophic", "Hit and Run w/ Runner Disruption", "Pure Hit and Run"]
available_movie_styles = ["disruption_headon", "disruption_off_axis", "supercatastrophic_headon", "supercatastrophic_off_axis","hitandrun_disrupt", "hitandrun_pure", "merge"]
movie_title_list = ["Head-on Disruption", "Off-axis Disruption", "Head-on Supercatastrophic", "Off-axis Supercatastrophic", "Hit and Run w/ Runner Disruption", "Pure Hit and Run", "Merge"]
movie_titles = dict(zip(available_movie_styles, movie_title_list))
num_movie_frames = 1200

Expand All @@ -54,12 +53,14 @@
"hitandrun_disrupt" : [np.array([1.0, -4.2e-05, 0.0]),
np.array([1.0, 4.2e-05, 0.0])],
"hitandrun_pure" : [np.array([1.0, -4.2e-05, 0.0]),
np.array([1.0, 4.2e-05, 0.0])]
np.array([1.0, 4.2e-05, 0.0])],
"merge" : [np.array([1.0, -5.0e-05, 0.0]),
np.array([1.0, 5.0e-05 ,0.0])]
}

vel_vectors = {"disruption_headon" : [np.array([ 0.00, 6.280005, 0.0]),
np.array([ 0.00, -6.280005, 0.0])],
"disruption_off_axis" : [np.array([ 0.00, 6.280005, 0.0]),
"disruption_off_axis" : [np.array([ 0.00, 6.280005, 0.0]),
np.array([ 0.50, -6.280005, 0.0])],
"supercatastrophic_headon": [np.array([ 0.00, 6.28, 0.0]),
np.array([ 0.00, -6.28, 0.0])],
Expand All @@ -68,7 +69,9 @@
"hitandrun_disrupt" : [np.array([ 0.00, 6.28, 0.0]),
np.array([-1.45, -6.28, 0.0])],
"hitandrun_pure" : [np.array([ 0.00, 6.28, 0.0]),
np.array([-1.51, -6.28, 0.0])]
np.array([-1.51, -6.28, 0.0])],
"merge" : [np.array([ 0.00, 0.0, 0.0]),
np.array([ 0.01, -0.100005, 0.0])]
}

rot_vectors = {"disruption_headon" : [np.array([0.0, 0.0, 0.0]),
Expand All @@ -82,17 +85,29 @@
"hitandrun_disrupt" : [np.array([0.0, 0.0, 6.0e3]),
np.array([0.0, 0.0, 1.0e4])],
"hitandrun_pure" : [np.array([0.0, 0.0, 6.0e3]),
np.array([0.0, 0.0, 1.0e4])]
np.array([0.0, 0.0, 1.0e4])],
"merge" : [np.array([0.0, 0.0, -6.0e3]),
np.array([0.0, 0.0, 1.0e4])]
}

body_Gmass = {"disruption_headon" : [1e-7, 1e-10],
"disruption_off_axis" : [1e-7, 1e-10],
"supercatastrophic_headon" : [1e-7, 1e-8],
"supercatastrophic_off_axis": [1e-7, 1e-8],
"hitandrun_disrupt" : [1e-7, 7e-10],
"hitandrun_pure" : [1e-7, 7e-10]
"hitandrun_pure" : [1e-7, 7e-10],
"merge" : [1e-7, 1e-8]
}

tstop = {"disruption_headon" : 5.0e-4,
"disruption_off_axis" : 5.0e-4,
"supercatastrophic_headon" : 5.0e-4,
"supercatastrophic_off_axis": 5.0e-4,
"hitandrun_disrupt" : 2.0e-4,
"hitandrun_pure" : 2.0e-4,
"merge" : 2.0e-3,
}

density = 3000 * swiftest.AU2M**3 / swiftest.MSun
GU = swiftest.GMSun * swiftest.YR2S**2 / swiftest.AU2M**3
body_radius = body_Gmass.copy()
Expand All @@ -114,7 +129,7 @@ def encounter_combiner(sim):
"""

# Only keep a minimal subset of necessary data from the simulation and encounter datasets
keep_vars = ['rh','Gmass','radius']
keep_vars = ['rh','vh','Gmass','radius']
data = sim.data[keep_vars]
enc = sim.encounters[keep_vars].load()

Expand Down Expand Up @@ -167,6 +182,9 @@ def setup_plot(self):
rhy2 = self.ds['rh'].sel(name="Projectile",space='y').isel(time=0).values[()]

scale_frame = abs(rhy1) + abs(rhy2)
if "hitandrun" in style:
scale_frame *= 2

ax = plt.Axes(fig, [0.1, 0.1, 0.8, 0.8])
self.ax_pt_size = self.figsize[0] * 72 / scale_frame
ax.set_xlim(-scale_frame, scale_frame)
Expand All @@ -182,30 +200,34 @@ def setup_plot(self):
return fig, ax

def update_plot(self, frame):
# Define a function to calculate the center of mass of the system.
def center(Gmass, x, y):
x = x[~np.isnan(x)]
y = y[~np.isnan(y)]
Gmass = Gmass[~np.isnan(Gmass)]
x_com = np.sum(Gmass * x) / np.sum(Gmass)
y_com = np.sum(Gmass * y) / np.sum(Gmass)
return x_com, y_com

Gmass, rh, point_rad = next(self.data_stream(frame))
x_com, y_com = center(Gmass, rh[:,0], rh[:,1])
self.scatter_artist.set_offsets(np.c_[rh[:,0] - x_com, rh[:,1] - y_com])

# Define a function to calculate a reference frame for the animation
# This will be based on the initial velocity of the Target body
def reference_frame(r_ref, v_ref, t):
coord_pos = r_ref + v_ref * t
return coord_pos.values[0], coord_pos.values[1]

t, Gmass, rh, point_rad = next(self.data_stream(frame))
x_ref, y_ref = reference_frame(self.r_ref,self.v_ref, t)
self.scatter_artist.set_offsets(np.c_[rh[:,0] - x_ref, rh[:,1] - y_ref])
self.scatter_artist.set_sizes(point_rad**2)
return self.scatter_artist,

def data_stream(self, frame=0):
while True:
ds = self.ds.isel(time=frame)
ds = ds.where(ds['name'] != "Sun", drop=True)
t = ds['time'].values[()]
radius = ds['radius'].values
Gmass = ds['Gmass'].values
rh = ds['rh'].values
point_rad = radius * self.ax_pt_size
yield Gmass, rh, point_rad

# Save the initial velocity of body 1 to use as a reference
if frame == 0:
self.r_ref = ds.sel(name="Target")['rh']
self.v_ref = ds.sel(name="Target")['vh']
yield t, Gmass, rh, point_rad

if __name__ == "__main__":

Expand All @@ -216,10 +238,11 @@ def data_stream(self, frame=0):
print("4. Off-axis supercatastrophic")
print("5. Hit and run with disruption of the runner")
print("6. Pure hit and run")
print("7. All of the above")
print("7. Merge")
print("8. All of the above")
user_selection = int(input("? "))

if user_selection > 0 and user_selection < 7:
if user_selection > 0 and user_selection < 8:
movie_styles = [available_movie_styles[user_selection-1]]
else:
print("Generating all movie styles")
Expand All @@ -234,10 +257,10 @@ def data_stream(self, frame=0):
sim.add_body(name=names, Gmass=body_Gmass[style], radius=body_radius[style], rh=pos_vectors[style], vh=vel_vectors[style], rot=rot_vectors[style])

# Set fragmentation parameters
minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body
gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades
minimum_fragment_gmass = 0.01 * body_Gmass[style][1]
gmtiny = 0.10 * body_Gmass[style][1]
sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False)
sim.run(dt=5e-4, tstop=5.0e-4, istep_out=1, dump_cadence=0)
sim.run(dt=5e-4, tstop=tstop[style], istep_out=1, dump_cadence=0)

print("Generating animation")
anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1)
19 changes: 8 additions & 11 deletions src/collision/collision_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -385,17 +385,14 @@ module subroutine collision_util_add_fragments_to_collider(self, nbody_system, p
class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters
end subroutine collision_util_add_fragments_to_collider

module subroutine collision_util_construct_after_system(collider, nbody_system, param, after_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
implicit none
! Arguments
class(collision_basic), intent(inout) :: collider !! Collision system object
class(base_nbody_system), intent(in) :: nbody_system !! Original Swiftest nbody system object
class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters
class(base_nbody_system), allocatable, intent(out) :: after_system !! Output temporary Swiftest nbody system object
end subroutine collision_util_construct_after_system
module subroutine collision_util_construct_constraint_system(collider, nbody_system, param, constraint_system, phase)
implicit none
class(collision_basic), intent(inout) :: collider !! Collision system object
class(base_nbody_system), intent(in) :: nbody_system !! Original Swiftest nbody system object
class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters
class(base_nbody_system), allocatable, intent(out) :: constraint_system !! Output temporary Swiftest nbody system object
character(len=*), intent(in) :: phase !! One of "before" or "after", indicating which phase of the calculation this needs to be done
end subroutine collision_util_construct_constraint_system

module subroutine collision_util_dealloc_fragments(self)
implicit none
Expand Down
Loading

0 comments on commit 0362420

Please sign in to comment.