diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index a5151d2d3..21ea232e6 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -469,11 +469,11 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu logical, intent(out) :: lfailure !! Did the velocity computation fail? ! Internals real(DP), parameter :: ENERGY_SUCCESS_METRIC = 1.0e-4_DP !! Relative energy error to accept as a success (success also must be energy-losing in addition to being within the metric amount) - 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) + real(DP) :: MOMENTUM_SUCCESS_METRIC = 2*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, dL, drot, rot_new, dL_metric, dL_best + 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 @@ -526,7 +526,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu lfailure = .false. dE_metric = huge(1.0_DP) dE_best = huge(1.0_DP) - dL_best = huge(1.0_DP) + dL_best(:) = huge(1.0_DP) nsteps_best = 0 nsteps = 0 outer: do try = 1, MAXOUTER @@ -575,11 +575,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(:) @@ -614,7 +614,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu 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,14 +623,14 @@ 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 dE_best = dE - dL_best = dL_metric + dL_best(:) = dL_metric(:) nsteps_best = nsteps do concurrent(i = 1:fragments%nbody) @@ -672,12 +672,12 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu 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) .or. (any(dL_best(:) > 1.0_DP)) 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 + write(message,*) "|dL|/|L0| = ",.mag.(dL_best) * MOMENTUM_SUCCESS_METRIC 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)