diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 6e80909c4..3eaa65ee4 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -342,7 +342,7 @@ def __init__(self,read_param: bool = False, read_data: bool = False, simdir: os. # If the file doesn't exist, flag it for now so we know to create it param_file_found = False if read_param or read_data: - if self.read_param(read_init_cond = read_data): + if self.read_param(read_init_cond = True): # We will add the parameter file to the kwarg list. This will keep the set_parameter method from # overriding everything with defaults when there are no arguments passed to Simulation() kwargs['param_file'] = self.param_file @@ -394,9 +394,9 @@ def _type_scrub(output_data): iloop = int((self.param['TSTART'] - self.param['T0']) / self.param['DT']) twidth = int(np.ceil(np.log10(self.param['TSTOP']/(self.param['DT'] * self.param['ISTEP_OUT'])))) pre_message = f"Time: {self.param['TSTART']:.{twidth}e} / {self.param['TSTOP']:.{twidth}e} {self.TU_name} " - post_message = f"npl: {self.data['npl'].values[0]} ntp: {self.data['ntp'].values[0]}" - if "nplm" in self.data: - post_message += f" nplm: {self.data['nplm'].values[0]}" + post_message = f"npl: {self.init_cond['npl'].values[0]} ntp: {self.init_cond['ntp'].values[0]}" + if "nplm" in self.init_cond: + post_message += f" nplm: {self.init_cond['nplm'].values[0]}" if self.param['ENERGY']: post_message += f" dL/L0: {0.0:.5e} dE/|E0|: {0.0:+.5e}" post_message += f" Wall time / step: {0.0:.5e} s" diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 4ef6ec215..f0fc6d9e7 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -172,7 +172,7 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) real(DP), intent(in) :: t !! The time of the collision ! Internals integer(I4B) :: i, j, k, ibiggest - real(DP), dimension(NDIM) :: L_spin_new + real(DP), dimension(NDIM) :: L_spin_new, L_residual real(DP) :: volume character(len=STRMAX) :: message @@ -228,6 +228,9 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) ! Get the energy of the system after the collision call self%get_energy_and_momentum(nbody_system, param, phase="after") + L_residual(:) = (self%L_total(:,2) - self%L_total(:,1)) + call collision_util_velocity_torque(-L_residual(:), fragments%mass(1), fragments%rb(:,1), fragments%vb(:,1)) + ! Update any encounter lists that have the removed bodies in them so that they instead point to the new body do k = 1, nbody_system%plpl_encounter%nenc do j = 1, impactors%ncoll diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index d6323f522..2957dffaf 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -521,6 +521,13 @@ module subroutine collision_util_setup_fragments(self, n) integer(I4B), intent(in) :: n !! Number of particles to allocate space for end subroutine collision_util_setup_fragments + module subroutine collision_util_velocity_torque(dL, mass, r, v) + implicit none + real(DP), dimension(:), intent(in) :: dL !! Change in angular momentum to apply + real(DP), intent(in) :: mass !! Mass of body + real(DP), dimension(:), intent(in) :: r !! Position of body wrt system center of mass + real(DP), dimension(:), intent(inout) :: v !! Velocity of body wrt system center of mass + end subroutine collision_util_velocity_torque end interface contains diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 63d623de4..7adbe2ec5 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -1028,4 +1028,32 @@ module subroutine collision_util_set_original_scale_factors(self) end subroutine collision_util_set_original_scale_factors + module subroutine collision_util_velocity_torque(dL, mass, r, v) + !! author: David A. Minton + !! + !! Applies a torque to a body's center of mass velocity given a change in angular momentum + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: dL !! Change in angular momentum to apply + real(DP), intent(in) :: mass !! Mass of body + real(DP), dimension(:), intent(in) :: r !! Position of body wrt system center of mass + real(DP), dimension(:), intent(inout) :: v !! Velocity of body wrt system center of mass + ! Internals + real(DP), dimension(NDIM) :: dL_unit, r_unit, r_lever, vapply + real(DP) :: rmag, r_lever_mag + + dL_unit(:) = .unit. dL + r_unit(:) = .unit.r(:) + rmag = .mag.r(:) + ! Project the position vector onto the plane defined by the angular momentum vector and the origin to get the "lever arm" distance + r_lever(:) = dL_unit(:) .cross. (r(:) .cross. dL_unit(:)) + r_lever_mag = .mag.r_lever(:) + if ((r_lever_mag > epsilon(1.0_DP)) .and. (rmag > epsilon(1.0_DP))) then + vapply(:) = (dL(:) .cross. r(:)) / (mass * rmag**2) + v(:) = v(:) + vapply(:) + end if + + return + end subroutine collision_util_velocity_torque + end submodule s_collision_util \ No newline at end of file diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index a5151d2d3..c91676621 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -26,6 +26,7 @@ module subroutine fraggle_generate(self, nbody_system, param, t) real(DP), intent(in) :: t !! Time of collision ! Internals integer(I4B) :: i, ibiggest, nfrag + real(DP), dimension(NDIM) :: L_residual, vbcom_orig, dvb character(len=STRMAX) :: message logical :: lfailure @@ -37,6 +38,9 @@ module subroutine fraggle_generate(self, nbody_system, param, t) ! Set the coordinate system of the impactors call impactors%set_coordinate_system() + + vbcom_orig(:) = impactors%vbcom(:) + select case (impactors%regime) case (COLLRESOLVE_REGIME_HIT_AND_RUN) call self%hitandrun(nbody_system, param, t) @@ -54,17 +58,14 @@ module subroutine fraggle_generate(self, nbody_system, param, t) end select call self%set_mass_dist(param) call self%disrupt(nbody_system, param, t, lfailure) - if (lfailure) then call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find an energy-losing solution. Simplifying the collisional model.") - impactors%mass_dist(1) = impactors%mass(1) impactors%mass_dist(2) = max(0.5_DP * impactors%mass(2), self%min_mfrag) impactors%mass_dist(3) = impactors%mass(2) - impactors%mass_dist(2) impactors%regime = COLLRESOLVE_REGIME_DISRUPTION call self%set_mass_dist(param) call self%disrupt(nbody_system, param, t, lfailure) - if (lfailure) then call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find an energy-losing solution. Treating this as a bounce.") call collision_util_bounce_one(impactors%rb(:,1),impactors%vb(:,1),impactors%rbcom(:),impactors%vbcom(:),impactors%radius(1)) @@ -77,16 +78,32 @@ module subroutine fraggle_generate(self, nbody_system, param, t) fragments%radius(1:2) = impactors%radius(1:2) fragments%rb(:,1:2) = impactors%rb(:,1:2) fragments%vb(:,1:2) = impactors%vb(:,1:2) + do concurrent(i = 1:2) + fragments%rc(:,i) = fragments%rb(:,i) - impactors%rbcom(:) + fragments%vc(:,i) = fragments%vb(:,i) - impactors%vbcom(:) + end do fragments%Ip(:,1:2) = impactors%Ip(:,1:2) fragments%rot(:,1:2) = impactors%rot(:,1:2) + fragments%mtot = sum(fragments%mass(1:2)) end associate - end if end if + ! Get the energy and momentum of the system before and after the collision + call self%get_energy_and_momentum(nbody_system, param, phase="before") + call self%get_energy_and_momentum(nbody_system, param, phase="after") + L_residual(:) = (self%L_total(:,2) - self%L_total(:,1)) + associate (fragments => self%fragments) - ! Populate the list of new bodies nfrag = fragments%nbody + + ! Put any residual angular momentum into orbital velocity + call collision_util_velocity_torque(-L_residual(:), fragments%mtot, impactors%rbcom(:), impactors%vbcom(:)) + dvb(:) = impactors%vbcom(:) - vbcom_orig(:) + do concurrent(i = 1:nfrag) + fragments%vb(:,i) = fragments%vb(:,i) + dvb(:) + end do + select case(impactors%regime) case(COLLRESOLVE_REGIME_DISRUPTION) status = DISRUPTED @@ -472,16 +489,16 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu real(DP) :: MOMENTUM_SUCCESS_METRIC = 10*epsilon(1.0_DP) !! Relative angular momentum error to accept as a success (should be *much* stricter than energy) integer(I4B) :: i, j, loop, try, istart, nfrag, nsteps, nsteps_best logical :: lhitandrun, lsupercat - real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual, L_residual_unit, dL, drot, rot_new - real(DP) :: vimp, vmag, vesc, dE, E_residual, E_residual_best, E_residual_last, ke_min, ke_avail, ke_remove, dE_best, fscale, dE_metric, mfrag, dL_metric, dL_best, rn + real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual, L_residual_unit, L_residual_best, dL, drot, rot_new, dL_metric + real(DP) :: vimp, vmag, vesc, dE, E_residual, E_residual_best, E_residual_last, ke_min, ke_avail, ke_remove, dE_best, fscale, dE_metric, mfrag, rn integer(I4B), dimension(collider%fragments%nbody) :: vsign real(DP), dimension(collider%fragments%nbody) :: vscale real(DP), parameter :: L_ROT_VEL_RATIO = 0.2_DP ! Ratio of angular momentum to put into rotation relative to velocity shear of fragments ! For the initial "guess" of fragment velocities, this is the minimum and maximum velocity relative to escape velocity that the fragments will have - real(DP) :: vmin_guess = 1.01_DP + real(DP) :: vmin_guess real(DP) :: vmax_guess real(DP) :: delta_v, GC - integer(I4B), parameter :: MAXINNER = 50 + integer(I4B), parameter :: MAXINNER = 100 integer(I4B), parameter :: MAXOUTER = 10 integer(I4B), parameter :: MAXANGMTM = 10000 class(collision_fraggle), allocatable :: collider_local @@ -521,12 +538,12 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu vimp = .mag. (impactors%vc(:,2) - impactors%vc(:,1)) vmax_guess = 1.1_DP * vimp + vmin_guess = 1.001_DP * vesc E_residual_best = huge(1.0_DP) lfailure = .false. dE_metric = huge(1.0_DP) dE_best = huge(1.0_DP) - dL_best = huge(1.0_DP) nsteps_best = 0 nsteps = 0 outer: do try = 1, MAXOUTER @@ -575,11 +592,11 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu angmtm: do j = 1, MAXANGMTM call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) - dL_metric = .mag.L_residual(:) / .mag.(collider_local%L_total(:,1)) + dL_metric(:) = abs(L_residual(:)) / .mag.(collider_local%L_total(:,1)) / MOMENTUM_SUCCESS_METRIC - if (dL_metric / MOMENTUM_SUCCESS_METRIC <= 1.0_DP) exit angmtm + if (all(dL_metric(:) <= 1.0_DP)) exit angmtm L_residual_unit(:) = .unit. L_residual(:) - do i = 1, fragments%nbody + do i = fragments%nbody,1,-1 mfrag = sum(fragments%mass(i:fragments%nbody)) drot(:) = -L_ROT_VEL_RATIO * L_residual(:) / (fragments%mtot * fragments%Ip(3,i) * fragments%radius(i)**2) rot_new(:) = fragments%rot(:,i) + drot(:) @@ -603,7 +620,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu end if dL(:) = -L_residual(:) * fragments%mass(i) / fragments%mtot - drot(:) * fragments%Ip(3,i) * fragments%mass(i) * fragments%radius(i)**2 - call fraggle_generate_velocity_torque(dL, fragments%mass(i), fragments%rc(:,i), fragments%vc(:,i)) + call collision_util_velocity_torque(dL, fragments%mass(i), fragments%rc(:,i), fragments%vc(:,i)) call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) fragments%vmag(i) = .mag.fragments%vc(:,i) @@ -613,8 +630,9 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) call fragments%set_coordinate_system() call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") + ke_avail = 0.0_DP - do i = 1, fragments%nbody + do i = fragments%nbody, 1, -1 ke_avail = ke_avail + 0.5_DP * fragments%mass(i) * max(fragments%vmag(i) - vesc,0.0_DP)**2 end do @@ -623,20 +641,16 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu E_residual = dE + impactors%Qloss L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) - dL_metric = .mag.L_residual(:) / .mag.collider_local%L_total(:,1) + dL_metric(:) = abs(L_residual(:)) / .mag.collider_local%L_total(:,1) / MOMENTUM_SUCCESS_METRIC ! Check if we've converged on our constraints - if (dL_metric < MOMENTUM_SUCCESS_METRIC) then + if (all(dL_metric(:) <= 1.0_DP)) then if ((abs(E_residual) < abs(E_residual_best)) .or. ((dE < 0.0_DP) .and. (dE_best >= 0.0_DP))) then ! This is our best case so far. Save it for posterity E_residual_best = E_residual + L_residual_best(:) = L_residual(:) dE_best = dE - dL_best = dL_metric nsteps_best = nsteps - do concurrent(i = 1:fragments%nbody) - fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:) - end do - if (allocated(collider%fragments)) deallocate(collider%fragments) allocate(collider%fragments, source=fragments) dE_metric = abs(E_residual) / impactors%Qloss @@ -668,22 +682,28 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu collider_local%fail_scale = collider_local%fail_scale * 1.001_DP ! Bring the minimum and maximum velocities closer together - delta_v = 0.125_DP * (vmax_guess - vmin_guess) + delta_v = (vmax_guess - vmin_guess) / 16.0_DP vmin_guess = vmin_guess + delta_v vmax_guess = vmax_guess - delta_v end do outer - lfailure = (dE_best > 0.0_DP) .or. (dL_best > MOMENTUM_SUCCESS_METRIC) + lfailure = (dE_best > 0.0_DP) + + call collision_util_velocity_torque(-L_residual_best(:), fragments%mtot, impactors%rbcom, impactors%vbcom) + + do concurrent(i = 1:collider%fragments%nbody) + collider%fragments%vb(:,i) = collider%fragments%vc(:,i) + impactors%vbcom(:) + end do write(message, *) nsteps if (lfailure) then call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle velocity calculation failed to converge after " // trim(adjustl(message)) // " steps. The best solution found had:") - write(message,*) "|dL|/|L0| = ",dL_best - call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) - write(message,*) " dE = ",dE_best - call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) else call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Fraggle velocity calculation converged after " // trim(adjustl(message)) // " steps.") end if + write(message,*) "dL/|L0| = ",(L_residual_best(:))/.mag.collider_local%L_total(:,1) + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + write(message,*) "dE/Qloss = ",-dE_best / impactors%Qloss + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) write(message,*) nsteps_best call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Best solution came after " // trim(adjustl(message)) // " steps.") @@ -693,33 +713,4 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu end subroutine fraggle_generate_vel_vec - subroutine fraggle_generate_velocity_torque(dL, mass, r, v) - !! author: David A. Minton - !! - !! Applies a torque to a body's center of mass velocity given a change in angular momentum - implicit none - ! Arguments - real(DP), dimension(:), intent(in) :: dL !! Change in angular momentum to apply - real(DP), intent(in) :: mass !! Mass of body - real(DP), dimension(:), intent(in) :: r !! Position of body wrt system center of mass - real(DP), dimension(:), intent(inout) :: v !! Velocity of body wrt system center of mass - ! Internals - real(DP), dimension(NDIM) :: dL_unit, r_unit, r_lever, vapply - real(DP) :: rmag, r_lever_mag - - dL_unit(:) = .unit. dL - r_unit(:) = .unit.r(:) - rmag = .mag.r(:) - ! Project the position vector onto the plane defined by the angular momentum vector and the origin to get the "lever arm" distance - r_lever(:) = dL_unit(:) .cross. (r(:) .cross. dL_unit(:)) - r_lever_mag = .mag.r_lever(:) - if ((r_lever_mag > epsilon(1.0_DP)) .and. (rmag > epsilon(1.0_DP))) then - vapply(:) = (dL(:) .cross. r(:)) / (mass * rmag**2) - v(:) = v(:) + vapply(:) - end if - - return - end subroutine fraggle_generate_velocity_torque - - end submodule s_fraggle_generate diff --git a/src/fraggle/fraggle_module.f90 b/src/fraggle/fraggle_module.f90 index 01412154d..fda3ebe1b 100644 --- a/src/fraggle/fraggle_module.f90 +++ b/src/fraggle/fraggle_module.f90 @@ -44,7 +44,7 @@ end subroutine fraggle_generate_disrupt module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t) implicit none - class(collision_fraggle), intent(inout) :: self !! Collision system object + class(collision_fraggle), intent(inout) :: self !! Collision system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions real(DP), intent(in) :: t !! Time of collision diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index af434b2c5..79daa0220 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -579,7 +579,7 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, param) integer(I4B) :: itmax, idmax, tslot real(DP), dimension(:), allocatable :: vals real(DP), dimension(1) :: rtemp - real(DP), dimension(NDIM) :: rot0, Ip0 + real(DP), dimension(NDIM) :: rot0, Ip0, L real(DP) :: mass0 associate (nc => self%system_history%nc, cb => self%cb) @@ -628,6 +628,8 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, param) rot0(:) = rot0(:) * DEG2RAD call netcdf_io_check( nf90_get_var(nc%id, nc%Ip_varid, Ip0, start=[1,1,tslot], count=[NDIM,1,1]), "netcdf_io_get_t0_values_system Ip_varid" ) cb%L0(:) = Ip0(3) * mass0 * cb%R0**2 * rot0(:) + L(:) = cb%Ip(3) * cb%mass * cb%radius**2 * cb%rot(:) + cb%dL(:) = L(:) - cb%L0 end if ! Retrieve the current bookkeeping variables @@ -2400,8 +2402,8 @@ module subroutine swiftest_io_param_writer(self, unit, iotype, v_list, iostat, i end if call io_param_writer_one("FIRSTKICK",param%lfirstkick, unit) - if (param%GMTINY > 0.0_DP) call io_param_writer_one("GMTINY",param%GMTINY, unit) - if (param%min_GMfrag > 0.0_DP) call io_param_writer_one("MIN_GMFRAG",param%min_GMfrag, unit) + if (param%GMTINY >= 0.0_DP) call io_param_writer_one("GMTINY",param%GMTINY, unit) + if (param%min_GMfrag >= 0.0_DP) call io_param_writer_one("MIN_GMFRAG",param%min_GMfrag, unit) call io_param_writer_one("COLLISION_MODEL",param%collision_model, unit) if (param%collision_model == "FRAGGLE" ) then nseeds = size(param%seed) diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 05f9383ee..74a795856 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -31,6 +31,8 @@ module subroutine swiftest_util_append_arr_char_string(arr, source, nold, nsrc, if (.not. allocated(source)) return nnew = count(lsource_mask(1:nsrc)) + if (nnew == 0) return + if (.not.allocated(arr)) then allocate(arr(nold+nnew)) else @@ -59,6 +61,8 @@ module subroutine swiftest_util_append_arr_DP(arr, source, nold, nsrc, lsource_m if (.not. allocated(source)) return nnew = count(lsource_mask(1:nsrc)) + if (nnew == 0) return + if (.not.allocated(arr)) then allocate(arr(nold+nnew)) else @@ -87,6 +91,8 @@ module subroutine swiftest_util_append_arr_DPvec(arr, source, nold, nsrc, lsourc if (.not. allocated(source)) return nnew = count(lsource_mask(1:nsrc)) + if (nnew == 0) return + if (.not.allocated(arr)) then allocate(arr(NDIM,nold+nnew)) else @@ -117,6 +123,8 @@ module subroutine swiftest_util_append_arr_I4B(arr, source, nold, nsrc, lsource_ if (.not. allocated(source)) return nnew = count(lsource_mask(1:nsrc)) + if (nnew == 0) return + if (.not.allocated(arr)) then allocate(arr(nold+nnew)) else @@ -146,6 +154,8 @@ module subroutine swiftest_util_append_arr_info(arr, source, nold, nsrc, lsource if (.not. allocated(source)) return nnew = count(lsource_mask(1:nsrc)) + if (nnew == 0) return + if (.not.allocated(arr)) then allocate(arr(nold+nnew)) else @@ -178,6 +188,8 @@ module subroutine swiftest_util_append_arr_kin(arr, source, nold, nsrc, lsource_ if (.not. allocated(source)) return nnew = count(lsource_mask(1:nsrc)) + if (nnew == 0) return + if (.not.allocated(arr)) then allocate(arr(nold+nnew)) else @@ -206,6 +218,8 @@ module subroutine swiftest_util_append_arr_logical(arr, source, nold, nsrc, lsou if (.not. allocated(source)) return nnew = count(lsource_mask(1:nsrc)) + if (nnew == 0) return + if (.not.allocated(arr)) then allocate(arr(nold+nnew)) else