From 0362420ebefc52e4ca80adb5714cf72e077fe1d1 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 10 Jan 2023 13:40:08 -0500 Subject: [PATCH] Made substantial improvements to the convergence of Fraggle by isolating 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. --- examples/Fragmentation/Fragmentation_Movie.py | 77 ++++++---- src/collision/collision_module.f90 | 19 ++- src/collision/collision_util.f90 | 132 ++++++++++-------- src/fraggle/fraggle_generate.f90 | 81 +++++++---- src/swiftest/swiftest_util.f90 | 52 +++---- 5 files changed, 206 insertions(+), 155 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 8d8e11698..a480064d1 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -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 @@ -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])], @@ -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]), @@ -82,7 +85,9 @@ "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], @@ -90,9 +95,19 @@ "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() @@ -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() @@ -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) @@ -182,18 +200,16 @@ 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, @@ -201,11 +217,17 @@ 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__": @@ -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") @@ -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) diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 3ce48a518..361a27bd7 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -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 diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 3f90f39c6..6e2627126 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -61,63 +61,77 @@ module subroutine collision_util_add_fragments_to_collider(self, nbody_system, p end subroutine collision_util_add_fragments_to_collider - module subroutine collision_util_construct_after_system(collider, nbody_system, param, after_system) + module subroutine collision_util_construct_constraint_system(collider, nbody_system, param, constraint_system, phase) !! Author: David A. Minton !! - !! Constructs a temporary internal system consisting of active bodies with impacting bodies replaced with additional fragments. This is used to compute the energy and momentum - !! of the system after the collision has been resolved + !! Constructs a temporary system that is used to evaluate the convergence on energy and angular momentum constraints. + !! The rotations of all bodies other than those involved in the collision are set to 0 so that the changes in spin kinetic + !! energy and momentum are isolated to the collision system. 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 + 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 ! Internals integer(I4B) :: i, status class(swiftest_nbody_system), allocatable :: tmpsys - !class(swiftest_parameters), allocatable :: tmpparam select type(nbody_system) class is (swiftest_nbody_system) select type(param) class is (swiftest_parameters) - associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody, pl => nbody_system%pl, npl => nbody_system%pl%nbody, cb => nbody_system%cb) - - ! Update barycentric vector values - do concurrent(i = 1:nfrag) - fragments%rb(:,i) = fragments%rc(:,i) + impactors%rbcom(:) - fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:) - end do + associate(pl => nbody_system%pl, npl => nbody_system%pl%nbody, cb => nbody_system%cb, nfrag => collider%fragments%nbody) ! Set up a new temporary system based on the original allocate(tmpsys, mold=nbody_system) allocate(tmpsys%cb, source=cb) allocate(tmpsys%pl, source=pl) - allocate(tmpsys%pl_adds, mold=pl) - allocate(tmpsys%pl_discards, mold=pl) allocate(tmpsys%collider, source=collider) call tmpsys%collider%set_original_scale() - select case(impactors%regime) - case(COLLRESOLVE_REGIME_DISRUPTION) - status = DISRUPTED - case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - status = SUPERCATASTROPHIC - case(COLLRESOLVE_REGIME_HIT_AND_RUN) - status = HIT_AND_RUN_DISRUPT - end select + ! Remove spins and velocities from all bodies other than the new fragments so that we can isolate the kinetic energy and momentum of the collision system, but still be able to compute + ! the potential energy correctly + tmpsys%cb%rot(:) = 0.0_DP + tmpsys%pl%rot(:,:) = 0.0_DP + tmpsys%pl%vb(:,:) = 0.0_DP - call collision_resolve_mergeaddsub(tmpsys, param, nbody_system%t, status) - call tmpsys%pl%rearray(tmpsys, param) + if (phase == "before") then + ! Put back the spins and velocities of the colliding bodies to compute pre-impact KE and L + tmpsys%pl%rot(:,collider%impactors%id(:)) = pl%rot(:,collider%impactors%id(:)) + tmpsys%pl%vb(:,collider%impactors%id(:)) = pl%vb(:,collider%impactors%id(:)) + else if (phase == "after") then + allocate(tmpsys%pl_adds, mold=pl) + allocate(tmpsys%pl_discards, mold=pl) + associate(impactors => tmpsys%collider%impactors, fragments => tmpsys%collider%fragments) ! Be sure to select the temporary version because its unit system has been updated + ! Update barycentric vector values + do concurrent(i = 1:nfrag) + fragments%rb(:,i) = fragments%rc(:,i) + impactors%rbcom(:) + fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:) + end do - call move_alloc(tmpsys, after_system) + select case(impactors%regime) + case(COLLRESOLVE_REGIME_DISRUPTION) + status = DISRUPTED + case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + status = SUPERCATASTROPHIC + case(COLLRESOLVE_REGIME_HIT_AND_RUN) + status = HIT_AND_RUN_DISRUPT + end select + end associate + + call collision_resolve_mergeaddsub(tmpsys, param, nbody_system%t, status) + call tmpsys%pl%rearray(tmpsys, param) + end if + call move_alloc(tmpsys, constraint_system) end associate end select end select return - end subroutine collision_util_construct_after_system + end subroutine collision_util_construct_constraint_system module subroutine collision_util_get_idvalues_snapshot(self, idvals) @@ -179,7 +193,7 @@ module subroutine collision_util_get_energy_and_momentum(self, nbody_system, par !! Calculates total system energy in either the pre-collision outcome state (phase = "before") or the post-collision outcome state (lbefore = .false.) !! This subrourtine works by building a temporary internal massive body object 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 massive body object in a temporary system object called after_system to feed it into symba_energy + !! This will temporarily expand the massive body object in a temporary system object called constraint_system to feed it into symba_energy implicit none ! Arguments class(collision_basic), intent(inout) :: self !! Encounter collision system object @@ -187,52 +201,52 @@ module subroutine collision_util_get_energy_and_momentum(self, nbody_system, par class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameters character(len=*), intent(in) :: phase !! One of "before" or "after", indicating which phase of the calculation this needs to be done ! Internals - class(base_nbody_system), allocatable :: after_system - integer(I4B) :: i + class(base_nbody_system), allocatable :: constraint_system + integer(I4B) :: i, phase_val + + select case(phase) + case("before") + phase_val = 1 + case("after") + phase_val = 2 + case default + write(*,*) "Unknown value of phase argument passed to collision_util_get_energy_and_momentum: ",trim(adjustl(phase)) + return + end select select type(nbody_system) class is (swiftest_nbody_system) select type(param) class is (swiftest_parameters) associate(fragments => self%fragments, impactors => self%impactors, nfrag => self%fragments%nbody, pl => nbody_system%pl, cb => nbody_system%cb) - if (phase == "before") then - call nbody_system%get_energy_and_momentum(param) - self%L_orbit(:,1) = nbody_system%L_orbit(:) / self%Lscale - self%L_spin(:,1) = nbody_system%L_spin(:) / self%Lscale - self%L_total(:,1) = nbody_system%L_total(:) / self%Lscale - self%ke_orbit(1) = nbody_system%ke_orbit / self%Escale - self%ke_spin(1) = nbody_system%ke_spin / self%Escale - self%pe(1) = nbody_system%pe / self%Escale - self%be(1) = nbody_system%be / self%Escale - self%te(1) = nbody_system%te / self%Escale - else if (phase == "after") then - call collision_util_construct_after_system(self, nbody_system, param, after_system) - select type(after_system) - class is (swiftest_nbody_system) - call after_system%get_energy_and_momentum(param) - self%L_orbit(:,2) = after_system%L_orbit(:) / self%Lscale - self%L_spin(:,2) = after_system%L_spin(:) / self%Lscale - self%L_total(:,2) = after_system%L_total(:) / self%Lscale - self%ke_orbit(2) = after_system%ke_orbit / self%Escale - self%ke_spin(2) = after_system%ke_spin / self%Escale - self%pe(2) = after_system%pe / self%Escale - self%be(2) = after_system%be / self%Escale - self%te(2) = after_system%te / self%Escale - + call collision_util_construct_constraint_system(self, nbody_system, param, constraint_system, phase) + select type(constraint_system) + class is (swiftest_nbody_system) + call constraint_system%get_energy_and_momentum(param) + self%L_orbit(:,phase_val) = constraint_system%L_orbit(:) / self%Lscale + self%L_spin(:,phase_val) = constraint_system%L_spin(:) / self%Lscale + self%L_total(:,phase_val) = constraint_system%L_total(:) / self%Lscale + self%ke_orbit(phase_val) = constraint_system%ke_orbit / self%Escale + self%ke_spin(phase_val) = constraint_system%ke_spin / self%Escale + self%pe(phase_val) = constraint_system%pe / self%Escale + self%be(phase_val) = constraint_system%be / self%Escale + self%te(phase_val) = constraint_system%te / self%Escale + + if (phase_val == 2) then do concurrent(i = 1:nfrag) fragments%ke_orbit(i) = 0.5_DP * fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i)) fragments%ke_spin(i) = 0.5_DP * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * dot_product(fragments%rot(:,i), fragments%rot(:,i)) fragments%L_orbit(:,i) = fragments%mass(i) * fragments%rc(:,i) .cross. fragments%vc(:,i) fragments%L_spin(:,i) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * fragments%rot(:,i) end do - call swiftest_util_get_potential_energy(nfrag, [(.true., i = 1, nfrag)], after_system%cb%Gmass, fragments%Gmass, fragments%mass, fragments%rb, fragments%pe) + call swiftest_util_get_potential_energy(nfrag, [(.true., i = 1, nfrag)], constraint_system%cb%Gmass, fragments%Gmass, fragments%mass, fragments%rb, fragments%pe) fragments%be = sum(-3*fragments%Gmass(:)*fragments%mass(:)/(5*fragments%radius(:))) fragments%L_orbit_tot(:) = sum(fragments%L_orbit(:,:),dim=2) fragments%L_spin_tot(:) = sum(fragments%L_spin(:,:),dim=2) fragments%ke_orbit_tot = sum(fragments%ke_orbit(:)) fragments%ke_spin_tot = sum(fragments%ke_spin(:)) - end select - end if + end if + end select end associate end select diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 13457740f..e7b0aabba 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -93,7 +93,8 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur ! Internals logical :: lk_plpl, lfailure_local logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes - real(DP) :: dE, dL + real(DP) :: dE + real(DP), dimension(NDIM) :: dL character(len=STRMAX) :: message real(DP), parameter :: fail_scale_initial = 1.001_DP @@ -128,16 +129,34 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur call fraggle_generate_rot_vec(self, nbody_system, param) call fraggle_generate_vel_vec(self, nbody_system, param, lfailure_local) call self%get_energy_and_momentum(nbody_system, param, phase="after") - call self%set_original_scale() - - dE = self%te(2) - self%te(1) - dL = .mag.(self%L_total(:,2) - self%L_total(:,1)) - write(message,*) dE - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Estimated energy change: " // trim(adjustl(message))) - write(message,*) dL - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Estimated angular momentum change: " // trim(adjustl(message))) + dL = self%L_total(:,2)- self%L_total(:,1) + dE = self%te(2) - self%te(1) + + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "All quantities in collision system natural units") + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "* Conversion factors (collision system units / nbody system units):") + write(message,*) "* Mass: ", self%mscale + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + write(message,*) "* Distance: ", self%dscale + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + write(message,*) "* Time: ", self%tscale + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + write(message,*) "* Velocity: ", self%vscale + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + write(message,*) "* Energy: ",self%Escale + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + write(message,*) "* Momentum: ", self%Lscale + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Energy constraint") + write(message,*) "Expected: Qloss = ", -impactors%Qloss + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + write(message,*) "Actual : dE = ",dE + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + write(message,*) "Actual : dL = ",dL + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + call self%set_original_scale() ! Restore the big array if (lk_plpl) call pl%flatten(param) @@ -394,14 +413,14 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu integer(I4B) :: i, j, loop, try, istart, nfrag logical :: lhitandrun, lsupercat real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual - real(DP) :: vmag, vesc, E_residual, ke_avail, ke_remove, E_residual_min, fscale, f_spin, f_orbit + real(DP) :: vmag, vesc, dE, E_residual, ke_avail, ke_remove, dE_best, E_residual_best, fscale, f_spin, f_orbit, dE_metric integer(I4B), dimension(collider%fragments%nbody) :: vsign real(DP), dimension(collider%fragments%nbody) :: vscale, ke_rot_remove real(DP), parameter :: vmin_initial_factor = 1.5_DP ! For the initial "guess" of fragment velocities, this is the maximum velocity relative to escape velocity that the fragments will have real(DP), parameter :: vmax_initial_factor = 5.0_DP ! For the initial "guess" of fragment velocities, this is the maximum velocity relative to escape velocity that the fragments will have - integer(I4B), parameter :: MAXLOOP = 100 + integer(I4B), parameter :: MAXLOOP = 200 integer(I4B), parameter :: MAXTRY = 20 - real(DP), parameter :: TOL = 1.0_DP + real(DP), parameter :: SUCCESS_METRIC = 1.0_DP class(collision_fraggle), allocatable :: collider_local character(len=STRMAX) :: message @@ -470,28 +489,35 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu ! Every time the collision-frame velocities are altered, we need to be sure to shift everything back to the center-of-mass frame call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) call fragments%set_coordinate_system() - E_residual_min = huge(1.0_DP) + E_residual_best = huge(1.0_DP) lfailure = .false. + dE_metric = huge(1.0_DP) outer: do try = 1, MAXTRY do loop = 1, MAXLOOP call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") ! Check for any residual angular momentum, and if there is any, put it into spin of the largest body - E_residual = collider_local%te(2) - collider_local%te(1) - E_residual = E_residual + impactors%Qloss L_residual(:) = collider_local%L_total(:,2) - collider_local%L_total(:,1) fragments%L_spin(:,1) = fragments%L_spin(:,1) - L_residual(:) fragments%rot(:,1) = fragments%L_spin(:,1) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(:,1)) call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") - E_residual = collider_local%te(2) + impactors%Qloss - collider_local%te(1) L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) - if ((abs(E_residual) < abs(E_residual_min)) .or. ((E_residual <= 0.0_DP) .and. (E_residual_min >= 0.0_DP))) then ! This is our best case so far. Save it for posterity + dE = collider_local%te(2) - collider_local%te(1) + E_residual = dE + impactors%Qloss + + if ((abs(E_residual) < abs(E_residual_best)) .or. ((dE < 0.0_DP) .and. (E_residual_best >= 0.0_DP))) then ! This is our best case so far. Save it for posterity + E_residual_best = E_residual + dE_best = dE + + do concurrent(i = 1:nfrag) + fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:) + end do + if (allocated(collider%fragments)) deallocate(collider%fragments) allocate(collider%fragments, source=fragments) - E_residual_min = E_residual - if (abs(E_residual) <= tiny(0.0_DP)) exit outer - if ((E_residual < 0.0_DP) .and. (impactors%Qloss / abs(E_residual)) > TOL) exit outer + dE_metric = abs(E_residual) / impactors%Qloss end if + if ((dE_best < 0.0_DP) .and. (dE_metric <= SUCCESS_METRIC * try)) exit outer ! As the tries increase, we relax the success metric. What was once a failure might become a success ke_avail = fragments%ke_orbit_tot - 0.5_DP * fragments%mtot * vesc**2 @@ -500,13 +526,15 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu f_orbit = 1.0_DP - f_spin ke_remove = min(f_orbit * E_residual, ke_avail) - fscale = sqrt((fragments%ke_orbit_tot - f_orbit * ke_remove)/fragments%ke_orbit_tot) + f_orbit = ke_remove / E_residual + fscale = sqrt((fragments%ke_orbit_tot - ke_remove)/fragments%ke_orbit_tot) fragments%vc(:,:) = fscale * fragments%vc(:,:) - ke_remove = min(E_residual - ke_remove, fragments%ke_spin_tot) + f_spin = 1.0_DP - f_orbit + ke_remove = min(f_spin * E_residual, fragments%ke_spin_tot) ke_rot_remove(:) = ke_remove * (fragments%ke_spin(:) / fragments%ke_spin_tot) where(ke_rot_remove(:) > fragments%ke_spin(:)) ke_rot_remove(:) = fragments%ke_spin(:) - do concurrent(i = 1:nfrag, fragments%ke_spin(i) > epsilon(1.0_DP)) + do concurrent(i = 1:nfrag, fragments%ke_spin(i) > 10*sqrt(tiny(1.0_DP))) fscale = sqrt((fragments%ke_spin(i) - ke_rot_remove(i))/fragments%ke_spin(i)) fragments%rot(:,i) = fscale * fragments%rot(:,i) end do @@ -521,11 +549,12 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu call fraggle_generate_rot_vec(collider_local, nbody_system, param) collider_local%fail_scale = collider_local%fail_scale*1.01_DP end do outer - lfailure = E_residual < 0.0_DP + lfailure = dE_best > 0.0_DP + if (lfailure) then - write(message,*) "Fraggle failed to converge after ",try*loop," steps. This collision may have added energy." + write(message,*) "Fraggle velocity calculation failed to converge after ",try*loop," steps. This collision may have added energy." else - write(message,*) "Fraggle succeeded after ",try*loop," steps." + write(message,*) "Fraggle velocity calculation converged after ",try*loop," steps." end if call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index fd37eaa3e..b83ffc8e1 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -1367,13 +1367,13 @@ module subroutine swiftest_util_get_energy_and_momentum_system(self, param) class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals - integer(I4B) :: i + integer(I4B) :: i,j real(DP) :: kecb, kespincb real(DP), dimension(self%pl%nbody) :: kepl, kespinpl - real(DP), dimension(self%pl%nbody) :: Lplorbitx, Lplorbity, Lplorbitz - real(DP), dimension(self%pl%nbody) :: Lplspinx, Lplspiny, Lplspinz + real(DP), dimension(NDIM,self%pl%nbody) :: Lplorbit + real(DP), dimension(NDIM,self%pl%nbody) :: Lplspin real(DP), dimension(NDIM) :: Lcborbit, Lcbspin - real(DP) :: hx, hy, hz + real(DP), dimension(NDIM) :: h associate(nbody_system => self, pl => self%pl, npl => self%pl%nbody, cb => self%cb) nbody_system%L_orbit(:) = 0.0_DP @@ -1383,12 +1383,8 @@ module subroutine swiftest_util_get_energy_and_momentum_system(self, param) nbody_system%ke_spin = 0.0_DP kepl(:) = 0.0_DP - Lplorbitx(:) = 0.0_DP - Lplorbity(:) = 0.0_DP - Lplorbitz(:) = 0.0_DP - Lplspinx(:) = 0.0_DP - Lplspiny(:) = 0.0_DP - Lplspinz(:) = 0.0_DP + Lplorbit(:,:) = 0.0_DP + Lplspin(:,:) = 0.0_DP pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE @@ -1397,14 +1393,10 @@ module subroutine swiftest_util_get_energy_and_momentum_system(self, param) Lcborbit(:) = cb%mass * (cb%rb(:) .cross. cb%vb(:)) do concurrent (i = 1:npl, pl%lmask(i)) - hx = pl%rb(2,i) * pl%vb(3,i) - pl%rb(3,i) * pl%vb(2,i) - hy = pl%rb(3,i) * pl%vb(1,i) - pl%rb(1,i) * pl%vb(3,i) - hz = pl%rb(1,i) * pl%vb(2,i) - pl%rb(2,i) * pl%vb(1,i) + h(:) = pl%rb(:,i) .cross. pl%vb(:,i) ! Angular momentum from orbit - Lplorbitx(i) = pl%mass(i) * hx - Lplorbity(i) = pl%mass(i) * hy - Lplorbitz(i) = pl%mass(i) * hz + Lplorbit(:,i) = pl%mass(i) * h(:) ! Kinetic energy from orbit kepl(i) = pl%mass(i) * dot_product(pl%vb(:,i), pl%vb(:,i)) @@ -1419,16 +1411,20 @@ module subroutine swiftest_util_get_energy_and_momentum_system(self, param) do concurrent (i = 1:npl, pl%lmask(i)) ! Currently we assume that the rotation pole is the 3rd principal axis ! Angular momentum from spin - Lplspinx(i) = pl%mass(i) * pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(1,i) - Lplspiny(i) = pl%mass(i) * pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(2,i) - Lplspinz(i) = pl%mass(i) * pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(3,i) + Lplspin(:,i) = pl%mass(i) * pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(:,i) ! Kinetic energy from spin kespinpl(i) = pl%mass(i) * pl%Ip(3,i) * pl%radius(i)**2 * dot_product(pl%rot(:,i), pl%rot(:,i)) end do + + nbody_system%ke_spin = 0.5_DP * (kespincb + sum(kespinpl(1:npl), pl%lmask(1:npl))) + + do concurrent (j = 1:NDIM) + nbody_system%L_spin(j) = Lcbspin(j) + sum(Lplspin(j,1:npl), pl%lmask(1:npl)) + end do else - kespincb = 0.0_DP - kespinpl(:) = 0.0_DP + nbody_system%ke_spin = 0.0_DP + nbody_system%L_spin(:) = 0.0_DP end if if (param%lflatten_interactions) then @@ -1444,20 +1440,12 @@ module subroutine swiftest_util_get_energy_and_momentum_system(self, param) end if nbody_system%ke_orbit = 0.5_DP * (kecb + sum(kepl(1:npl), pl%lmask(1:npl))) - if (param%lrotation) nbody_system%ke_spin = 0.5_DP * (kespincb + sum(kespinpl(1:npl), pl%lmask(1:npl))) - nbody_system%L_orbit(1) = Lcborbit(1) + sum(Lplorbitx(1:npl), pl%lmask(1:npl)) - nbody_system%L_orbit(2) = Lcborbit(2) + sum(Lplorbity(1:npl), pl%lmask(1:npl)) - nbody_system%L_orbit(3) = Lcborbit(3) + sum(Lplorbitz(1:npl), pl%lmask(1:npl)) - - if (param%lrotation) then - nbody_system%L_spin(1) = Lcbspin(1) + sum(Lplspinx(1:npl), pl%lmask(1:npl)) - nbody_system%L_spin(2) = Lcbspin(2) + sum(Lplspiny(1:npl), pl%lmask(1:npl)) - nbody_system%L_spin(3) = Lcbspin(3) + sum(Lplspinz(1:npl), pl%lmask(1:npl)) - end if + do concurrent (j = 1:NDIM) + nbody_system%L_orbit(j) = Lcborbit(j) + sum(Lplorbit(j,1:npl), pl%lmask(1:npl)) + end do nbody_system%be = sum(-3*pl%Gmass(1:npl)*pl%mass(1:npl)/(5*pl%radius(1:npl)), pl%lmask(1:npl)) - nbody_system%te = nbody_system%ke_orbit + nbody_system%ke_spin + nbody_system%pe + nbody_system%be nbody_system%L_total(:) = nbody_system%L_orbit(:) + nbody_system%L_spin(:) end associate