From 4767b5ee5afbfe1151e37d4aba26f9dfdcc12bca Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 16 Aug 2021 17:54:54 -0400 Subject: [PATCH] Fixed restructuring code to do a better job interpolating bad tries to converge on good tries --- src/fragmentation/fragmentation.f90 | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index 6f9f2bfee..169c42e2c 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -35,7 +35,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, real(DP) :: Etot_before, ke_orbit_before, ke_spin_before, pe_before, Lmag_before real(DP) :: Etot_after, ke_orbit_after, ke_spin_after, pe_after, Lmag_after, dEtot, dLmag real(DP), dimension(NDIM) :: L_frag_tot, L_frag_orb, L_frag_spin, L_frag_budget, Lorbit_before, Lorbit_after, Lspin_before, Lspin_after, dL - real(DP) :: ke_frag_budget, ke_frag_orbit, ke_frag_spin, ke_avg_deficit, ke_avg_deficit_old + real(DP) :: ke_frag_budget, ke_frag_orbit, ke_frag_spin, ke_tot_deficit, ke_avg_deficit, ke_avg_deficit_old real(DP), dimension(NDIM) :: x_col_unit, y_col_unit, z_col_unit character(len=*), parameter :: fmtlabel = "(A14,10(ES11.4,1X,:))" integer(I4B) :: try @@ -58,11 +58,8 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, fpe_quiet_modes(:) = .false. call ieee_set_halting_mode(IEEE_ALL,fpe_quiet_modes) - f_spin = 0.05_DP - allocate(x_frag, source=xb_frag) allocate(v_frag, source=vb_frag) - allocate(rmag(nfrag)) allocate(rotmag(nfrag)) allocate(v_r_mag(nfrag)) @@ -99,6 +96,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, L_frag_orb(:) = L_orb(:, 1) + L_orb(:, 2) L_frag_spin(:) = L_spin(:, 1) + L_spin(:, 2) L_frag_budget(:) = L_frag_orb(:) + L_frag_spin(:) + f_spin = norm2(L_frag_spin(:)) / norm2(L_frag_budget(:)) call define_coordinate_system() call construct_temporary_system() @@ -109,7 +107,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, r_max_start = 1 * norm2(x(:,2) - x(:,1)) try = 1 lfailure = .false. - ke_avg_deficit = 0.0_DP + ke_tot_deficit = 0.0_DP do while (try < MAXTRY) lfailure = .false. call reset_fragments() @@ -130,7 +128,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, if (lfailure) then write(*,*) 'Failed to find tangential velocities' - ke_avg_deficit = ke_avg_deficit - (ke_frag_orbit + ke_frag_spin) + else call set_fragment_radial_velocities(lfailure) if (lfailure) write(*,*) 'Failed to find radial velocities' @@ -969,10 +967,12 @@ subroutine restructure_failed_fragments() real(DP) :: delta_r, delta_r_max real(DP), parameter :: ke_avg_deficit_target = 0.0_DP + ke_tot_deficit = ke_tot_deficit - (ke_frag_budget - ke_frag_orbit - ke_frag_spin) + ke_avg_deficit = ke_tot_deficit / try ! Introduce a bit of noise in the radius determination so we don't just flip flop between similar failed positions call random_number(delta_r_max) delta_r_max = sum(radius(:)) * (1.0_DP + 2e-1_DP * (delta_r_max - 0.5_DP)) - if (try > 2) then + if (try > 1) then ! Linearly interpolate the last two failed solution ke deficits to find a new distance value to try delta_r = (r_max_start - r_max_start_old) * (ke_avg_deficit_target - ke_avg_deficit_old) / (ke_avg_deficit - ke_avg_deficit_old) if (abs(delta_r) > delta_r_max) delta_r = sign(delta_r_max, delta_r) @@ -981,6 +981,7 @@ subroutine restructure_failed_fragments() end if r_max_start_old = r_max_start r_max_start = r_max_start + delta_r ! The larger lever arm can help if the problem is in the angular momentum step + ke_avg_deficit_old = ke_avg_deficit if (f_spin > epsilon(1.0_DP)) then ! Try reducing the fraction in spin f_spin = f_spin / 2