diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 88ced4615..ebad58734 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -232,7 +232,7 @@ def data_stream(self, frame=0): # Set fragmentation parameters minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades - sim.set_parameter(collision_model="simple", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) + sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) sim.run(dt=1e-3, tstop=1.0e-3, istep_out=1, dump_cadence=0) print("Generating animation") diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 3683346fb..978dfdcff 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -141,7 +141,9 @@ module subroutine collision_generate_hitandrun(self, nbody_system, param, t) call self%set_mass_dist(param) ! Generate the position and velocity distributions of the fragments + call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) call self%disrupt(nbody_system, param, t, lpure) + call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) dpe = self%pe(2) - self%pe(1) nbody_system%Ecollisions = nbody_system%Ecollisions - dpe @@ -227,7 +229,9 @@ module subroutine collision_generate_simple(self, nbody_system, param, t) call util_exit(FAILURE) end select call self%set_mass_dist(param) + call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) call self%disrupt(nbody_system, param, t) + call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) dpe = self%pe(2) - self%pe(1) nbody_system%Ecollisions = nbody_system%Ecollisions - dpe @@ -376,12 +380,11 @@ module subroutine collision_generate_disrupt(self, nbody_system, param, t, lfail logical, optional, intent(out) :: lfailure ! Internals - call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) call collision_generate_simple_pos_vec(self) call self%set_coordinate_system() call collision_generate_simple_vel_vec(self) call collision_generate_simple_rot_vec(self) - call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) + return end subroutine collision_generate_disrupt @@ -481,28 +484,34 @@ module subroutine collision_generate_simple_rot_vec(collider) ! Arguments class(collision_basic), intent(inout) :: collider !! Fraggle collision system object ! Internals - real(DP), dimension(NDIM) :: Lresidual + real(DP), dimension(NDIM) :: Lresidual, Lbefore, Lafter, Lspin, rot integer(I4B) :: i associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) fragments%rot(:,:) = 0.0_DP - ! Keep the first two bodies spinning as before to start with - Lresidual(:) = 0.0_DP - do i = 1,2 - fragments%rot(:,i) = impactors%Lspin(:,i) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) - Lresidual(:) = Lresidual(:) + impactors%Lorbit(:,i) - end do + ! Torque the first body based on the change in angular momentum betwen the pre- and post-impact system + Lbefore(:) = impactors%mass(2) * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (impactors%vb(:,2) - impactors%vb(:,1)) - ! Compute the current orbital angular momentum - do i = 1, nfrag - Lresidual(:) = Lresidual(:) - fragments%mass(i) * (fragments%rc(:,i) .cross. fragments%vc(:,i)) + Lafter(:) = 0.0_DP + do i = 2, nfrag + Lafter(:) = Lafter(:) + fragments%mass(i) * (fragments%rb(:,i) - fragments%rb(:,1)) .cross. (fragments%vb(:,i) - fragments%vb(:,1)) end do + Lspin(:) = impactors%Lspin(:,1) + (Lbefore(:) - Lafter(:)) + + fragments%rot(:,1) = Lspin(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1)) + + Lresidual(:) = sum(impactors%Lspin(:,:) + impactors%Lorbit(:,:), dim=2) - Lspin(:) + + ! Randomize the rotational vector direction of the n>1 fragments and distribute the residual momentum amongst them + call random_number(fragments%rot(:,2:nfrag)) + rot(:) = Lresidual(:) / sum(fragments%mass(2:nfrag) * fragments%radius(2:nfrag)**2 * fragments%Ip(3,2:nfrag)) + + fragments%rot(:,2:nfrag) = .unit. fragments%rot(:,2:nfrag) * .mag. rot(:) - ! Distributed most of the remaining angular momentum amongst all the particles if (.mag.(Lresidual(:)) > tiny(1.0_DP)) then - do i = 1,nfrag - fragments%rot(:,i) = fragments%rot(:,i) + Lresidual(:) / (nfrag * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) + do i = 2,nfrag + fragments%rot(:,i) = fragments%rot(:,i) + Lresidual(:) / ((nfrag - 1) * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) end do end if diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index c0df619c5..cfac7b49f 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -498,8 +498,7 @@ subroutine collision_final_fragments(self) ! Arguments type(collision_fragments(*)), intent(inout) :: self - call self%reset() - + if (allocated(self%info)) deallocate(self%info) return end subroutine collision_final_fragments diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 76b48ea9e..4a8ee6b74 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -191,8 +191,7 @@ module subroutine collision_util_get_energy_momentum(self, nbody_system, param, ! Internals class(base_nbody_system), allocatable, save :: tmpsys class(base_parameters), allocatable, save :: tmpparam - integer(I4B) :: npl_before, npl_after, stage - + integer(I4B) :: npl_before, npl_after, stage,i select type(nbody_system) class is (swiftest_nbody_system) select type(param) @@ -226,7 +225,11 @@ module subroutine collision_util_get_energy_momentum(self, nbody_system, param, ! Build the exluded body logical mask for the *after* case: Only the new bodies are used to compute energy and momentum call self%add_fragments(tmpsys, tmpparam) tmpsys%pl%status(impactors%id(1:impactors%ncoll)) = INACTIVE - tmpsys%pl%status(npl_before+1:npl_after) = ACTIVE + do concurrent(i = npl_before+1:npl_after) + tmpsys%pl%status(i) = ACTIVE + tmpsys%pl%rot(:,i) = 0.0_DP + tmpsys%pl%vb(:,i) = impactors%vbcom(:) + end do end select end if select type(tmpsys) @@ -234,7 +237,8 @@ module subroutine collision_util_get_energy_momentum(self, nbody_system, param, if (param%lflatten_interactions) call tmpsys%pl%flatten(param) - call tmpsys%get_energy_and_momentum(param) + call tmpsys%get_energy_and_momentum(param) + ! Calculate the current fragment energy and momentum balances if (lbefore) then @@ -248,8 +252,7 @@ module subroutine collision_util_get_energy_momentum(self, nbody_system, param, self%ke_orbit(stage) = tmpsys%ke_orbit self%ke_spin(stage) = tmpsys%ke_spin self%pe(stage) = tmpsys%pe - self%Etot(stage) = tmpsys%te - if (stage == 2) self%Etot(stage) = self%Etot(stage) - (self%pe(2) - self%pe(1)) ! Gotta be careful with PE when number of bodies changes. + self%Etot(stage) = tmpsys%ke_orbit + tmpsys%ke_spin end select end associate end select diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 21a297c08..7ca7cac06 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -12,9 +12,8 @@ use symba 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), parameter :: F_SPIN_FIRST = 0.05_DP !! The initial try value of the fraction of energy or momenum in spin (whichever has the lowest kinetic energy) - real(DP), parameter :: FRAGGLE_LTOL = 10 * epsilon(1.0_DP) - real(DP), parameter :: FRAGGLE_ETOL = 1e-8_DP + real(DP), parameter :: FRAGGLE_LTOL = 1e-4_DP !10 * epsilon(1.0_DP) + real(DP), parameter :: FRAGGLE_ETOL = 1e-1_DP contains @@ -31,10 +30,9 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur real(DP), intent(in) :: t !! Time of collision logical, optional, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? ! Internals - integer(I4B) :: i integer(I4B) :: try - real(DP) :: r_max_start, f_spin, dEtot, dLmag - integer(I4B), parameter :: MAXTRY = 100 + real(DP) :: dEtot, dLmag + integer(I4B), parameter :: MAXTRY = 10 logical :: lk_plpl, lfailure_local logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes logical, dimension(size(IEEE_USUAL)) :: fpe_flag @@ -69,66 +67,38 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur else lk_plpl = .false. end if + call ieee_set_flag(ieee_all, .false.) ! Set all fpe flags to quiet call self%set_natural_scale() call fragments%reset() - ! Calculate the initial energy of the nbody_system without the collisional family - call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) - - ! Start out the fragments close to the initial separation distance. This will be increased if there is any overlap or we fail to find a solution - r_max_start = 1.1_DP * .mag.(impactors%rb(:,2) - impactors%rb(:,1)) lfailure_local = .false. try = 1 - f_spin = F_SPIN_FIRST + do while (try < MAXTRY) write(message,*) try call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle try " // trim(adjustl(message))) if (lfailure_local) then - call fragments%restructure(impactors, try, f_spin, r_max_start) + !call fragments%restructure(impactors, try, f_spin, r_max_start) call fragments%reset() try = try + 1 end if - lfailure_local = .false. - call ieee_set_flag(ieee_all, .false.) ! Set all fpe flags to quiet - - call collision_generate_simple_pos_vec(self) - call self%set_coordinate_system() - - ! Initial velocity guess will be the barycentric velocity of the colliding nbody_system so that the budgets are based on the much smaller collisional-frame velocities - do concurrent (i = 1:nfrag) - fragments%vb(:, i) = impactors%vbcom(:) - end do + ! Use the simple collision model to generate initial conditions + ! Compute the "before" energy/momentum and the budgets + call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) + call self%collision_simple_disruption%disrupt(nbody_system, param, t) call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) call self%set_budgets() - call collision_generate_simple_vel_vec(self) - - call fraggle_generate_tan_vel(self, lfailure_local) - if (lfailure_local) then - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find tangential velocities") - cycle - end if - - call fraggle_generate_rad_vel(self, lfailure_local) - if (lfailure_local) then - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find radial velocities") - cycle - end if - - call fraggle_generate_spins(self, f_spin, lfailure_local) - if (lfailure_local) then - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find spins") - cycle - end if + ! Minimize difference between energy/momentum and budgets + call fraggle_generate_minimize(self, lfailure_local) call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) dEtot = self%Etot(2) - self%Etot(1) dLmag = .mag. (self%Ltot(:,2) - self%Ltot(:,1)) - exit lfailure_local = ((abs(dEtot + impactors%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) if (lfailure_local) then @@ -136,6 +106,8 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to high energy error: " // & trim(adjustl(message))) cycle + lfailure_local = .false. + exit end if lfailure_local = ((abs(dLmag) / (.mag.self%Ltot(:,1))) > FRAGGLE_LTOL) @@ -178,408 +150,143 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur end subroutine fraggle_generate_disrupt - subroutine fraggle_generate_spins(collider, f_spin, lfailure) - !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Calculates the spins of a collection of fragments such that they conserve angular momentum without blowing the fragment kinetic energy budget. - !! - !! A failure will trigger a restructuring of the fragments so we will try new values of the radial position distribution. - implicit none - ! Arguments - class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object - real(DP), intent(in) :: f_spin !! Fraction of energy or momentum that goes into spin (whichever gives the lowest kinetic energy) - logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds! - ! Internals - real(DP), dimension(NDIM) :: L_remainder, rot_L, rot_ke, L - real(DP), dimension(NDIM,collider%fragments%nbody) :: frot_rand ! The random rotation factor applied to fragments - real(DP), parameter :: frot_rand_mag = 1.50_DP ! The magnitude of the rotation variation to apply to the fragments - integer(I4B) :: i - character(len=STRMAX) :: message - real(DP) :: ke_remainder, ke - - associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) - select type(fragments => collider%fragments) - class is (fraggle_fragments(*)) - lfailure = .false. - L_remainder(:) = fragments%L_budget(:) - ke_remainder = fragments%ke_budget - - ! Add a fraction of the orbit angular momentum of the second body to the spin of the first body to kick things off - L(:) = impactors%Lspin(:,1) + f_spin * (impactors%Lorbit(:,2) + impactors%Lspin(:,2)) - fragments%rot(:,1) = L(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1)) - L_remainder(:) = L_remainder(:) - L(:) - - ! Partition the spin momentum of the second body into the spin of the fragments, with some random variation - L(:) = impactors%Lspin(:,2) / (nfrag - 1) - - call random_number(frot_rand(:,2:nfrag)) - frot_rand(:,2:nfrag) = 2 * (frot_rand(:,2:nfrag) - 0.5_DP) * frot_rand_mag - - do i = 2, nfrag - rot_L(:) = L(:) + frot_rand(:,i) * .mag.L(:) - fragments%rot(:,i) = rot_L(:) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,1)) - L_remainder(:) = L_remainder(:) - rot_L(:) - end do - - ! Make sure we didn't blow our kinetic energy budget - do i = 1, nfrag - ke_remainder = ke_remainder - 0.5_DP * fragments%mass(i) * fragments%Ip(3, i) * fragments%radius(i)**2 * .mag.(fragments%rot(:, i)) - end do - - ! Distributed most of the remaining angular momentum amongst all the particles - fragments%ke_spin = 0.0_DP - if (.mag.(L_remainder(:)) > FRAGGLE_LTOL) then - do i = nfrag, 1, -1 - ! 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_remainder / (i * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3, i))) * L_remainder(:) / .mag.(L_remainder(:)) - rot_L(:) = f_spin * L_remainder(:) / (i * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3, i)) - if (.mag.(rot_ke) < .mag.(rot_L)) then - fragments%rot(:,i) = fragments%rot(:,i) + rot_ke(:) - else - fragments%rot(:, i) = fragments%rot(:,i) + rot_L(:) - end if - ke = 0.5_DP * fragments%mass(i) * fragments%Ip(3, i) * fragments%radius(i)**2 * norm2(fragments%rot(:, i)) - L(:) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3, i) * fragments%rot(:, i) - ke_remainder = ke_remainder - ke - L_remainder(:) = L_remainder(:) - L(:) - fragments%ke_spin = fragments%ke_spin + ke - end do - end if - - lfailure = ((fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit) < 0.0_DP) - - if (lfailure) then - call swiftest_io_log_one_message(COLLISION_LOG_OUT, " ") - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Spin failure diagnostics") - write(message, *) fragments%ke_budget - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_budget : " // trim(adjustl(message))) - write(message, *) fragments%ke_spin - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_spin : " // trim(adjustl(message))) - write(message, *) fragments%ke_orbit - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_orbit : " // trim(adjustl(message))) - write(message, *) fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_remainder : " // trim(adjustl(message))) - end if - - end select - end associate - - return - end subroutine fraggle_generate_spins - - - subroutine fraggle_generate_tan_vel(collider, lfailure) + subroutine fraggle_generate_minimize(collider, lfailure) !! 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 (fragments%ke_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. + !! Minimizes the differences between the energy and momentum and the budget !! - !! A failure will trigger a restructuring of the fragments so we will try new values of the radial position distribution. + !! A failure will trigger a restructuring of the fragments so we will try a new random distribution implicit none ! Arguments class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds ! Internals - integer(I4B) :: i, try - real(DP), parameter :: TOL_MIN = 1e0_DP ! This doesn't have to be very accurate, as we really just want a tangential velocity distribution with less kinetic energy than our initial guess. - real(DP), parameter :: TOL_INIT = 1e-12_DP - real(DP), parameter :: VNOISE_MAG = 1e-1_DP !! Magnitude of the noise to apply to initial conditions to help minimizer find a solution in case of failure - integer(I4B), parameter :: MAXLOOP = 10 - integer(I4B), parameter :: MAXTRY = 100 - real(DP) :: tol - real(DP), dimension(:), allocatable :: v_t_initial, v_t_output - real(DP), dimension(collider%fragments%nbody) :: kefrag, vnoise - type(lambda_obj_err) :: objective_function - real(DP), dimension(NDIM) :: L_frag_tot - character(len=STRMAX) :: message - real(DP) :: ke_diff + real(DP), parameter :: TOL_MIN = 1.0e-6_DP + real(DP), parameter :: TOL_INIT = 1e-8_DP + integer(I4B), parameter :: MAXLOOP = 30 + real(DP), dimension(collider%fragments%nbody) :: input_v + real(DP), dimension(:), allocatable :: output_v + type(lambda_obj) :: Efunc + real(DP) :: tol, fval + integer(I4B) :: i, nelem + + associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) + select type(fragments => collider%fragments) + class is (fraggle_fragments(*)) + + nelem = 2 * (nfrag - 1) + lfailure = .false. + ! Find the local kinetic energy minimum for the nbody_system that conserves linear and angular momentum + Efunc = lambda_obj(E_objective_function) + tol = TOL_INIT - associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) - select type(fragments => collider%fragments) - class is (fraggle_fragments(*)) - lfailure = .false. + do while(tol < TOL_MIN) - allocate(v_t_initial, mold=fragments%v_t_mag) - do try = 1, MAXTRY - v_t_initial(1) = dot_product(impactors%vb(:,1),fragments%v_t_unit(:,1)) - do i = 2, nfrag - v_t_initial(i) = dot_product(impactors%vb(:,2), fragments%v_t_unit(:,i)) - end do - fragments%v_t_mag(:) = v_t_initial + fragments%v_r_unit(:,:) = .unit. fragments%vc(:,:) + fragments%vmag(:) = .mag. fragments%vc(:,:) + + input_v(:) = fragments%vmag(1:nfrag) + fval = E_objective_function(input_v) + call minimize_bfgs(Efunc, nelem, input_v, tol, MAXLOOP, lfailure, output_v) + fval = E_objective_function(output_v) + input_v(:) = output_v(:) - ! Find the local kinetic energy minimum for the nbody_system that conserves linear and angular momentum - objective_function = lambda_obj(tangential_objective_function, lfailure) + fragments%vmag(1:nfrag) = output_v(1:nfrag) - tol = TOL_INIT - do while(tol < TOL_MIN) - call minimize_bfgs(objective_function, nfrag-6, v_t_initial(7:nfrag), tol, MAXLOOP, lfailure, v_t_output) - fragments%v_t_mag(7:nfrag) = v_t_output(:) - ! 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) = fragments%v_t_mag(7:nfrag) - if (.not.lfailure) exit - tol = tol * 2_DP ! Keep increasing the tolerance until we converge on a solution - call random_number(vnoise(1:nfrag)) ! Adding a bit of noise to the initial conditions helps it find a solution more often - vnoise(:) = 1.0_DP + VNOISE_MAG * (2 * vnoise(:) - 1._DP) - v_t_initial(:) = v_t_initial(:) * vnoise(:) - end do - fragments%v_t_mag(1:nfrag) = solve_fragment_tan_vel(v_t_mag_input=v_t_initial(7:nfrag), lfailure=lfailure) + do concurrent(i=2:nfrag) + fragments%vc(:,i) = abs(fragments%vmag(i)) * fragments%v_r_unit(:,i) + end do - ! Perform one final shift of the radial velocity vectors to align with the center of mass of the collisional system (the origin) - fragments%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%v_r_unit(:,1:nfrag), fragments%v_t_mag(1:nfrag), & - fragments%v_t_unit(:,1:nfrag), fragments%mass(1:nfrag), impactors%vbcom(:)) - do concurrent (i = 1:nfrag) - fragments%vc(:,i) = fragments%vb(:,i) - impactors%vbcom(:) - end do + call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) + ! Set spins in order to conserve angular momentum + call fraggle_generate_set_spins(fragments) - ! 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) - kefrag(i) = fragments%mass(i) * dot_product(fragments%vb(:, i), fragments%vb(:, i)) + if (.not.lfailure) exit + tol = tol * 2 ! Keep increasing the tolerance until we converge on a solution end do - fragments%ke_orbit = 0.5_DP * sum(kefrag(:)) - - ! If we are over the energy budget, flag this as a failure so we can try again - ke_diff = fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit - lfailure = ke_diff < 0.0_DP - if (.not.lfailure) exit - fragments%rc(:,:) = fragments%rc(:,:) * 1.1_DP - end do - if (lfailure) then - call swiftest_io_log_one_message(COLLISION_LOG_OUT, " ") - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Tangential velocity failure diagnostics") - call fragments%get_angular_momentum() - L_frag_tot = fragments%Lspin(:) + fragments%Lorbit(:) - write(message, *) .mag.(fragments%L_budget(:) - L_frag_tot(:)) / (.mag.collider%Ltot(:,1)) - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "|L_remainder| : " // trim(adjustl(message))) - write(message, *) fragments%ke_budget - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_budget : " // trim(adjustl(message))) - write(message, *) fragments%ke_spin - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_spin : " // trim(adjustl(message))) - write(message, *) fragments%ke_orbit - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_tangential : " // trim(adjustl(message))) - write(message, *) fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_radial : " // trim(adjustl(message))) - end if - end select + end select end associate return contains - function solve_fragment_tan_vel(lfailure, 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) :: lfailure !! 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 - - select type(fragments => collider%fragments) - class is (fraggle_fragments(*)) - associate(nfrag => fragments%nbody) - lfailure = .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) = fragments%mass(i) * fragments%v_t_unit(:, i) - A(4:6, i) = fragments%mass(i) * fragments%rmag(i) * (fragments%v_r_unit(:, i) .cross. fragments%v_t_unit(:, i)) - else if (present(v_t_mag_input)) then - vtmp(:) = v_t_mag_input(i - 6) * fragments%v_t_unit(:, i) - L_lin_others(:) = L_lin_others(:) + fragments%mass(i) * vtmp(:) - L(:) = fragments%mass(i) * (fragments%rc(:, i) .cross. vtmp(:)) - L_orb_others(:) = L_orb_others(:) + L(:) - end if - end do - b(1:3) = -L_lin_others(:) - b(4:6) = fragments%L_budget(:) - fragments%Lspin(:) - L_orb_others(:) - allocate(v_t_mag_output(nfrag)) - v_t_mag_output(1:6) = solve_linear_system(A, b, 6, lfailure) - if (present(v_t_mag_input)) v_t_mag_output(7:nfrag) = v_t_mag_input(:) - end associate - end select - return - end function solve_fragment_tan_vel - - function tangential_objective_function(v_t_mag_input, lfailure) result(fval) + function E_objective_function(val_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 + !! Objective function for evaluating how close our fragment trajectories are to the energy budget 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) :: lfailure !! Error flag + real(DP), dimension(:), intent(in) :: val_input !! Flattened velocity and rotation arrays ! Result - real(DP) :: fval + real(DP) :: fval !! The objective function result, which is the sum of the squares of the difference between the calculated fragment kinetic energy and the components of angular and linear momentum + !! Minimizing this brings us closer to our objective ! Internals integer(I4B) :: i - real(DP), dimension(NDIM,collider%fragments%nbody) :: v_shift - real(DP), dimension(collider%fragments%nbody) :: v_t_new, kearr - real(DP) :: keo - - select type(fragments => collider%fragments) - class is (fraggle_fragments(*)) - associate(impactors => collider%impactors, nfrag => fragments%nbody) - lfailure = .false. - - v_t_new(:) = solve_fragment_tan_vel(v_t_mag_input=v_t_mag_input(:), lfailure=lfailure) - v_shift(:,:) = fraggle_util_vmag_to_vb(fragments%v_r_mag, fragments%v_r_unit, v_t_new, fragments%v_t_unit, fragments%mass, impactors%vbcom) - - kearr = 0.0_DP - do concurrent(i = 1:nfrag) - kearr(i) = fragments%mass(i) * dot_product(v_shift(:, i), v_shift(:, i)) - end do - keo = 0.5_DP * sum(kearr(:)) - fval = keo - lfailure = .false. + type(fraggle_fragments(:)), allocatable :: tmp_frag + real(DP) :: deltaE + + associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) + select type(fragments => collider%fragments) + class is (fraggle_fragments(*)) + allocate(tmp_frag, source=fragments) + tmp_frag%vmag(1:nfrag) = val_input(1:nfrag) + do concurrent(i=1:nfrag) + tmp_frag%vc(:,i) = abs(tmp_frag%vmag(i)) * tmp_frag%v_r_unit(:,i) + end do + + call collision_util_shift_vector_to_origin(tmp_frag%mass, tmp_frag%vc) + ! Set spins in order to conserve angular momentum + call fraggle_generate_set_spins(tmp_frag) + + ! Get the current kinetic energy of the system + call fragments%get_kinetic_energy() + deltaE = fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin) + + ! Use the deltaE as the basis of our objective function, with a higher penalty for having excess kinetic energy compared with having a deficit + if (deltaE < 0.0_DP) then + fval = deltaE**4 + else + fval = deltaE**2 + end if + + end select end associate - end select - + return - end function tangential_objective_function + end function E_objective_function - end subroutine fraggle_generate_tan_vel + end subroutine fraggle_generate_minimize - subroutine fraggle_generate_rad_vel(collider, lfailure) - !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + + subroutine fraggle_generate_set_spins(fragments) + !! Author: 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 + !! Distributes any residual angular momentum into the spins of the n>1 fragments implicit none ! Arguments - class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object - logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds! + class(fraggle_fragments(*)), intent(inout) :: fragments !! Fraggle fragment system object ! Internals - real(DP), parameter :: TOL_MIN = FRAGGLE_ETOL ! This needs to be more accurate than the tangential step, as we are trying to minimize the total residual energy - real(DP), parameter :: TOL_INIT = 1e-14_DP - real(DP), parameter :: VNOISE_MAG = 1e-10_DP !! Magnitude of the noise to apply to initial conditions to help minimizer find a solution in case of failure - integer(I4B), parameter :: MAXLOOP = 100 - real(DP) :: ke_radial, tol - integer(I4B) :: i - real(DP), dimension(:), allocatable :: v_r_initial, v_r_output - real(DP), dimension(collider%fragments%nbody) :: vnoise - type(lambda_obj) :: objective_function - character(len=STRMAX) :: message - - associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) - select type(fragments => collider%fragments) - class is (fraggle_fragments(*)) - ! Set the "target" ke for the radial component - - allocate(v_r_initial, source=fragments%v_r_mag) - ! Initialize radial velocity magnitudes with a random value that related to equipartition of kinetic energy with some noise and scaled with respect to the initial distance - v_r_initial(1) = dot_product(impactors%vb(:,1),fragments%v_r_unit(:,1)) - fragments%ke_orbit = 0.5_DP * fragments%mass(1) * v_r_initial(1)**2 - - ke_radial = fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit - call random_number(vnoise(1:nfrag)) - vnoise(:) = 1.0_DP + VNOISE_MAG * (2 * vnoise(:) - 1.0_DP) - v_r_initial(2:nfrag) = -sqrt(abs(2 * ke_radial) / (fragments%mass(1:nfrag) * nfrag)) * vnoise(1: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) - tol = TOL_INIT - do while(tol < TOL_MIN) - call minimize_bfgs(objective_function, nfrag, v_r_initial, tol, MAXLOOP, lfailure, v_r_output) - fragments%v_r_mag(1:nfrag) = v_r_output(1:nfrag) - if (.not.lfailure) exit - tol = tol * 2 ! Keep increasing the tolerance until we converge on a solution - v_r_initial(:) = fragments%v_r_mag(:) - call random_number(vnoise(1:nfrag)) ! Adding a bit of noise to the initial conditions helps it find a solution more often - vnoise(:) = 1.0_DP + VNOISE_MAG * (2 * vnoise(:) - 1._DP) - v_r_initial(:) = v_r_initial(:) * vnoise(:) - end do - - ! Shift the radial velocity vectors to align with the center of mass of the collisional system (the origin) - fragments%ke_orbit = 0.0_DP - fragments%ke_spin = 0.0_DP - fragments%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%v_r_unit(:,1:nfrag), & - fragments%v_t_mag(1:nfrag), fragments%v_t_unit(:,1:nfrag), fragments%mass(1:nfrag), impactors%vbcom(:)) - do i = 1, nfrag - fragments%vc(:, i) = fragments%vb(:, i) - impactors%vbcom(:) - fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * norm2(fragments%vb(:, i)) - fragments%ke_spin = fragments%ke_spin + fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * norm2(fragments%rot(:,i)) + integer(I4B) :: i + real(DP), dimension(NDIM) :: Lresidual + + call fragments%get_angular_momentum() + Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit(:) + fragments%Lspin(:)) + + ! Distribute residual angular momentum amongst the fragments + if (.mag.(Lresidual(:)) > tiny(1.0_DP)) then + do i = 2,fragments%nbody + fragments%rot(:,i) = fragments%rot(:,i) + Lresidual(:) / ((fragments%nbody - 1) * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) end do - fragments%ke_orbit = 0.5_DP * fragments%ke_orbit - fragments%ke_spin = 0.5_DP * fragments%ke_spin - - lfailure = abs((fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin)) / fragments%ke_budget) > FRAGGLE_ETOL - if (lfailure) then - call swiftest_io_log_one_message(COLLISION_LOG_OUT, " ") - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Radial velocity failure diagnostics") - write(message, *) fragments%ke_budget - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_budget : " // trim(adjustl(message))) - write(message, *) fragments%ke_spin - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_spin : " // trim(adjustl(message))) - write(message, *) fragments%ke_orbit - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_orbit : " // trim(adjustl(message))) - write(message, *) fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin) - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_remainder : " // trim(adjustl(message))) - end if + end if - end select - end associate return + end subroutine fraggle_generate_set_spins + + - contains - 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(collider%fragments%nbody) :: kearr - real(DP) :: keo, ke_radial, rotmag2, vmag2 - - associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) - select type(fragments => collider%fragments) - class is (fraggle_fragments(*)) - allocate(v_shift, mold=fragments%vb) - v_shift(:,:) = fraggle_util_vmag_to_vb(v_r_mag_input, fragments%v_r_unit, fragments%v_t_mag, fragments%v_t_unit, fragments%mass, impactors%vbcom) - do i = 1,fragments%nbody - rotmag2 = fragments%rot(1,i)**2 + fragments%rot(2,i)**2 + fragments%rot(3,i)**2 - vmag2 = v_shift(1,i)**2 + v_shift(2,i)**2 + v_shift(3,i)**2 - kearr(i) = fragments%mass(i) * (fragments%Ip(3, i) * fragments%radius(i)**2 * rotmag2 + vmag2) - end do - keo = 2 * fragments%ke_budget - sum(kearr(:)) - ke_radial = fragments%ke_budget - fragments%ke_orbit - fragments%ke_spin - ! 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 - end select - end associate - - return - end function radial_objective_function - end subroutine fraggle_generate_rad_vel end submodule s_fraggle_generate diff --git a/src/fraggle/fraggle_module.f90 b/src/fraggle/fraggle_module.f90 index fe0eaf8a3..f2ad578c4 100644 --- a/src/fraggle/fraggle_module.f90 +++ b/src/fraggle/fraggle_module.f90 @@ -27,10 +27,11 @@ module fraggle real(DP) :: ke_budget !! Kinetic energy budget for computing fragment trajectories real(DP), dimension(NDIM) :: L_budget !! Angular momentum budget for computing fragment trajectories contains - procedure :: get_angular_momentum => fraggle_util_get_angular_momentum !! Calcualtes the current angular momentum of the fragments - procedure :: reset => fraggle_util_reset_fragments !! Resets all position and velocity-dependent fragment quantities in order to do a fresh calculation (does not reset mass, radius, or other values that get set prior to the call to fraggle_generate) - procedure :: restructure => fraggle_util_restructure !! Restructure the inputs after a failed attempt failed to find a set of positions and velocities that satisfy the energy and momentum constraints - final :: fraggle_final_fragments !! Finalizer will deallocate all allocatables + procedure :: get_angular_momentum => fraggle_util_get_angular_momentum !! Calcualtes the current angular momentum of the fragments + procedure :: get_kinetic_energy => fraggle_util_get_kinetic_energy !! Calcualtes the current kinetic energy of the fragments + procedure :: reset => fraggle_util_reset_fragments !! Resets all position and velocity-dependent fragment quantities in order to do a fresh calculation (does not reset mass, radius, or other values that get set prior to the call to fraggle_generate) + procedure :: restructure => fraggle_util_restructure !! Restructure the inputs after a failed attempt failed to find a set of positions and velocities that satisfy the energy and momentum constraints + final :: fraggle_final_fragments !! Finalizer will deallocate all allocatables end type fraggle_fragments @@ -84,6 +85,11 @@ module subroutine fraggle_util_get_angular_momentum(self) class(fraggle_fragments(*)), intent(inout) :: self !! Fraggle fragment system object end subroutine fraggle_util_get_angular_momentum + module subroutine fraggle_util_get_kinetic_energy(self) + implicit none + class(fraggle_fragments(*)), intent(inout) :: self !! Fraggle fragment system object + end subroutine fraggle_util_get_kinetic_energy + module subroutine fraggle_util_reset_fragments(self) implicit none class(fraggle_fragments(*)), intent(inout) :: self @@ -139,7 +145,7 @@ subroutine fraggle_final_fragments(self) ! Arguments type(fraggle_fragments(*)), intent(inout) :: self !! Fraggle encountar storage object - call self%collision_fragments%reset() + if (allocated(self%info)) deallocate(self%info) return end subroutine fraggle_final_fragments diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 43823642c..07d845053 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -40,7 +40,7 @@ end subroutine fraggle_util_construct_temporary_system module subroutine fraggle_util_get_angular_momentum(self) !! Author: David A. Minton !! - !! Calcualtes the current angular momentum of the fragments + !! Calculates the current angular momentum of the fragments implicit none ! Arguments class(fraggle_fragments(*)), intent(inout) :: self !! Fraggle fragment system object @@ -61,6 +61,33 @@ module subroutine fraggle_util_get_angular_momentum(self) end subroutine fraggle_util_get_angular_momentum + module subroutine fraggle_util_get_kinetic_energy(self) + !! Author: David A. Minton + !! + !! Calculates the current kinetic energy of the fragments + implicit none + ! Argument + class(fraggle_fragments(*)), intent(inout) :: self !! Fraggle fragment system object + ! Internals + integer(I4B) :: i + + associate(fragments => self, nfrag => self%nbody) + fragments%ke_orbit = 0.0_DP + fragments%ke_spin = 0.0_DP + + do i = 1, nfrag + fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i)) + fragments%ke_spin = fragments%ke_spin + fragments%mass(i) * fragments%Ip(3,i) * dot_product(fragments%rot(:,i),fragments%rot(:,i) ) + end do + + fragments%ke_orbit = fragments%ke_orbit / 2 + fragments%ke_spin = fragments%ke_spin / 2 + + end associate + + return + end subroutine fraggle_util_get_kinetic_energy + module subroutine fraggle_util_reset_fragments(self) !! author: David A. Minton !! @@ -177,7 +204,7 @@ module subroutine fraggle_util_set_budgets(self) dL(:) = self%Ltot(:,2) - self%Ltot(:,1) fragments%L_budget(:) = -dL(:) - fragments%ke_budget = -(dEtot - 0.5_DP * fragments%mtot * dot_product(impactors%vbcom(:), impactors%vbcom(:))) - impactors%Qloss + fragments%ke_budget = -(dEtot - impactors%Qloss) end select end associate @@ -214,7 +241,10 @@ module subroutine fraggle_util_set_natural_scale_factors(self) impactors%bounce_unit(:) = impactors%bounce_unit(:) / collider%vscale impactors%rb(:,:) = impactors%rb(:,:) / collider%dscale impactors%vb(:,:) = impactors%vb(:,:) / collider%vscale + impactors%rc(:,:) = impactors%rc(:,:) / collider%dscale + impactors%vc(:,:) = impactors%vc(:,:) / collider%vscale impactors%mass(:) = impactors%mass(:) / collider%mscale + impactors%Mcb = impactors%Mcb / collider%mscale impactors%radius(:) = impactors%radius(:) / collider%dscale impactors%Lspin(:,:) = impactors%Lspin(:,:) / collider%Lscale impactors%Lorbit(:,:) = impactors%Lorbit(:,:) / collider%Lscale @@ -256,11 +286,16 @@ module subroutine fraggle_util_set_original_scale_factors(self) impactors%rbimp(:) = impactors%rbimp(:) * collider%dscale impactors%bounce_unit(:) = impactors%bounce_unit(:) * collider%vscale - impactors%mass = impactors%mass * collider%mscale - impactors%radius = impactors%radius * collider%dscale - impactors%rb = impactors%rb * collider%dscale - impactors%vb = impactors%vb * collider%vscale - impactors%Lspin = impactors%Lspin * collider%Lscale + impactors%mass = impactors%mass * collider%mscale + impactors%Mcb = impactors%Mcb * collider%mscale + impactors%mass_dist = impactors%mass_dist * collider%mscale + impactors%radius = impactors%radius * collider%dscale + impactors%rb = impactors%rb * collider%dscale + impactors%vb = impactors%vb * collider%vscale + impactors%rc = impactors%rc * collider%dscale + impactors%vc = impactors%vc * collider%vscale + impactors%Lspin = impactors%Lspin * collider%Lscale + impactors%Lorbit = impactors%Lorbit * collider%Lscale do i = 1, 2 impactors%rot(:,i) = impactors%Lspin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i)) end do diff --git a/src/misc/minimizer_module.f90 b/src/misc/minimizer_module.f90 index 772db418a..667f061db 100644 --- a/src/misc/minimizer_module.f90 +++ b/src/misc/minimizer_module.f90 @@ -128,7 +128,7 @@ module subroutine minimize_bfgs(f, N, x0, eps, maxloop, lerr, x1) real(DP), dimension(:), intent(out), allocatable :: x1 ! Internals integer(I4B) :: i, j, k, l, conv - real(DP), parameter :: graddelta = 1e-4_DP !! Delta x for gradient calculations + real(DP), parameter :: graddelta = 1e-6_DP !! Delta x for gradient calculations real(DP), dimension(N) :: S !! Direction vectors real(DP), dimension(N,N) :: H !! Approximated inverse Hessian matrix real(DP), dimension(N) :: grad1 !! gradient of f @@ -137,6 +137,7 @@ module subroutine minimize_bfgs(f, N, x0, eps, maxloop, lerr, x1) real(DP), dimension(N) :: y, P real(DP), dimension(N,N) :: PP, PyH, HyP real(DP), save :: yHy, Py + real(DP) :: Hnorm call ieee_get_status(original_fpe_status) ! Save the original floating point exception status call ieee_set_flag(ieee_all, .false.) ! Set all flags to quiet @@ -181,7 +182,7 @@ module subroutine minimize_bfgs(f, N, x0, eps, maxloop, lerr, x1) end do !$omp end do simd ! prevent divide by zero (convergence) - if (abs(Py) < tiny(Py)) exit + if (abs(Py) < N**2 * tiny(Py)) exit ! set up update PyH(:,:) = 0._DP HyP(:,:) = 0._DP @@ -201,13 +202,21 @@ module subroutine minimize_bfgs(f, N, x0, eps, maxloop, lerr, x1) ! update H matrix H(:,:) = H(:,:) + ((1._DP - yHy / Py) * PP(:,:) - PyH(:,:) - HyP(:,:)) / Py ! Normalize to prevent it from blowing up if it takes many iterations to find a solution - H(:,:) = H(:,:) / norm2(H(:,:)) + Hnorm = 0.0_DP + do concurrent (j = 1:N,k=1:N,abs(H(j,k))>sqrt(10*tiny(1.0_DP))) + Hnorm = Hnorm + H(j,k)**2 + end do + Hnorm = sqrt(Hnorm) ! Stop everything if there are any exceptions to allow the routine to fail gracefully call ieee_get_flag(ieee_usual, fpe_flag) if (any(fpe_flag)) exit if (i == maxloop) then lerr = .true. end if + where(abs(H(:,:)) < sqrt(10*tiny(1.0_DP)) ) + H(:,:) = 0.0_DP + endwhere + H(:,:) = H(:,:) / Hnorm end do call ieee_get_flag(ieee_usual, fpe_flag) lerr = lerr .or. any(fpe_flag)