diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index ac0ef2b27..7562d7b89 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -46,7 +46,7 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi 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) :: r_max_start real(DP), parameter :: Ltol = 10 * epsilon(1.0_DP) - real(DP), parameter :: Etol = 1e-2_DP + real(DP), parameter :: Etol = 1e-10_DP integer(I4B), parameter :: MAXTRY = 200 logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes @@ -99,8 +99,11 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi if (lmerge) write(*,*) 'Failed to find radial velocities' if (.not.lmerge) then call calculate_system_energy(linclude_fragments=.true.) - if ((abs((dEtot - Qloss) / Qloss) > Etol) .or. (dEtot > 0.0_DP)) then - write(*,*) 'Failed due to high energy error: ', abs((dEtot - Qloss) / Qloss), dEtot / abs(Etot_before) + write(*,*) 'Qloss : ',Qloss + write(*,*) '-dEtot: ',-dEtot + write(*,*) 'delta : ',abs((dEtot + Qloss)) + if ((abs(dEtot + Qloss) > Etol) .or. (dEtot > 0.0_DP)) then + write(*,*) 'Failed due to high energy error: ' lmerge = .true. else if (abs(dLmag) / Lmag_before > Ltol) then write(*,*) 'Failed due to high angular momentum error: ', dLmag / Lmag_before @@ -156,7 +159,7 @@ subroutine set_scale_factors() ! Set scale factors mscale = mtot !! Because of the implied G, mass is actually G*mass with units of distance**3 / time**2 Escale = ke_family - vscale = sqrt(2 * Escale / mscale) + vscale = sqrt(Escale / mscale) rscale = mscale**2 / Escale tscale = rscale / vscale Lscale = mscale * rscale * vscale @@ -292,15 +295,15 @@ subroutine define_pre_collisional_family() !! Defines the pre-collisional "family" consisting of all of the bodies involved in the collision. These bodies need to be identified so that they can be excluded from energy calculations once !! the collisional products (the fragments) are created. integer(I4B) :: i - associate(vbpl => symba_plA%helio%swiftest%vb, Mpl => symba_plA%helio%swiftest%mass) + associate(vbpl => symba_plA%helio%swiftest%vb, Mpl => symba_plA%helio%swiftest%mass, & + Ipl => symba_plA%helio%swiftest%Ip, radpl => symba_plA%helio%swiftest%radius, rotpl => symba_plA%helio%swiftest%rot) fam_size = size(family) ! We need the original kinetic energy of just the pre-impact family members in order to balance the energy later ke_family = 0.0_DP do i = 1, fam_size - ke_family = ke_family + Mpl(family(i)) * dot_product(vbpl(:,family(i)), vbpl(:,family(i))) + ke_family = ke_family + 0.5_DP * Mpl(family(i)) * (dot_product(vbpl(:,family(i)), vbpl(:,family(i))) + Ipl(3,family(i)) * radpl(family(i))**2 * dot_product(rotpl(:,family(i)), rotpl(:, family(i)))) end do - ke_family = 0.5_DP * ke_family end associate return end subroutine define_pre_collisional_family @@ -370,8 +373,8 @@ subroutine calculate_system_energy(linclude_fragments) call move_alloc(ltmp, lexclude) end if - where (lexclude(1:npl)) - symba_plwksp%helio%swiftest%status(1:npl) = INACTIVE + where (lexclude(1:npl_new)) + symba_plwksp%helio%swiftest%status(1:npl_new) = INACTIVE end where nplm = count(symba_plwksp%helio%swiftest%mass > param%mtiny / mscale) @@ -527,6 +530,21 @@ subroutine set_fragment_tangential_velocities(lerr) call calculate_system_energy(linclude_fragments=.true.) ke_frag_budget = ke_family + (ke_spin_before - ke_spin_after) + (pe_before - pe_after) - Qloss + !ke_frag_budget = -dEtot - Qloss + write(*,*) '***************************************************' + write(*,*) 'Energy balance so far: ' + write(*,*) 'ke_family : ',ke_family + write(*,*) 'ke_frag_budget : ',ke_frag_budget + write(*,*) 'dEtot : ',dEtot + write(*,*) 'ke_frag_budget + dEtot ~ 0?', ke_frag_budget + dEtot + write(*,*) '***************************************************' + write(*,*) 'ke_orbit_before: ',ke_orbit_before + write(*,*) 'ke_orbit_after : ',ke_orbit_after + write(*,*) 'ke_spin_before : ',ke_spin_before + write(*,*) 'ke_spin_after : ',ke_spin_after + write(*,*) 'pe_before : ',pe_before + write(*,*) 'pe_after : ',pe_after + write(*,*) 'Qloss : ',Qloss if (ke_frag_budget < 0.0_DP) then write(*,*) 'Negative ke_frag_budget: ',ke_frag_budget lerr = .true. @@ -535,7 +553,7 @@ subroutine set_fragment_tangential_velocities(lerr) allocate(v_t_initial, mold=v_t_mag) - f_spin = 0.5_DP + f_spin = 0.0_DP L_frag_spin(:) = 0.0_DP do i = 1, nfrag rot_frag(:,i) = sqrt(2 * f_spin * ke_frag_budget / (nfrag * m_frag(i) * rad_frag(i)**2 * Ip_frag(3, i))) * v_h_unit(:, i) @@ -548,7 +566,7 @@ subroutine set_fragment_tangential_velocities(lerr) do i = 7, nfrag call util_crossproduct(v_r_unit(:, i), z_col_unit, r_cross_L(:)) rbar = rmag(i) * norm2(r_cross_L(:)) - v_t_initial(i) = L_orb_mag / (m_frag(i) * rbar * nfrag) + v_t_initial(i) = 0.5_DP * L_orb_mag / (m_frag(i) * rbar * nfrag) end do objective_function = lambda_obj(tangential_objective_function, lerr) v_t_mag(7:nfrag) = util_minimize_bfgs(objective_function, nfrag-6, v_t_initial(7:nfrag), TOL, lerr) @@ -676,7 +694,7 @@ subroutine set_fragment_radial_velocities(lerr) ! Arguments logical, intent(out) :: lerr ! Internals - real(DP), parameter :: TOL = 1e-12_DP + real(DP), parameter :: TOL = 2e-8_DP integer(I4B) :: i, j real(DP), dimension(:), allocatable :: v_r_initial, v_r_sigma real(DP), dimension(:,:), allocatable :: v_r @@ -716,6 +734,7 @@ subroutine set_fragment_radial_velocities(lerr) end do ke_frag = 0.5_DP * ke_frag write(*,*) 'Radial' + write(*,*) 'Failure? ',lerr write(*,*) 'ke_frag_budget: ',ke_frag_budget write(*,*) 'ke_frag : ',ke_frag write(*,*) 'ke_frag_spin : ',ke_frag_spin @@ -745,7 +764,7 @@ function radial_objective_function(v_r_mag_input) result(fval) fval = fval - m_frag(i) * (Ip_frag(3, i) * rad_frag(i)**2 * dot_product(rot_frag(:, i), rot_frag(:, i)) + dot_product(v_shift(:, i), v_shift(:, i))) end do ! The following ensures that fval = 0 is a local minimum, which is what the BFGS method is searching for - fval = (fval / (2 * ke_frag_budget))**2 + fval = (fval / (2 * ke_radial))**2 return diff --git a/src/util/util_minimize_bfgs.f90 b/src/util/util_minimize_bfgs.f90 index 585d7dc76..2e12640a9 100644 --- a/src/util/util_minimize_bfgs.f90 +++ b/src/util/util_minimize_bfgs.f90 @@ -60,13 +60,13 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) !check for convergence conv = count(abs(grad1(:)) > eps) if (conv == 0) then - !write(*,*) "BFGS converged on gradient after ",i," iterations" + write(*,*) "BFGS converged on gradient after ",i," iterations" exit end if S(:) = -matmul(H(:,:), grad1(:)) astar = minimize1D(f, x1, S, N, gradeps, lerr) if (lerr) then - !write(*,*) "Exiting BFGS with error in minimize1D step" + write(*,*) "Exiting BFGS with error in minimize1D step" exit end if ! Get new x values @@ -86,7 +86,7 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) end do ! prevent divide by zero (convergence) if (abs(Py) < tiny(Py)) then - !write(*,*) "BFGS Converged on tiny Py after ",i," iterations" + write(*,*) "BFGS Converged on tiny Py after ",i," iterations" exit end if ! set up update @@ -110,12 +110,12 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) if (any(fpe_flag)) exit if (i == MAXLOOP) then lerr = .true. - !write(*,*) "BFGS ran out of loops!" + write(*,*) "BFGS ran out of loops!" end if end do call ieee_get_flag(ieee_usual, fpe_flag) lerr = lerr .or. any(fpe_flag) - !if (any(fpe_flag)) write(*,*) 'BFGS did not converge due to fpe' + if (any(fpe_flag)) write(*,*) 'BFGS did not converge due to fpe' !if (lerr) write(*,*) "BFGS did not converge!" call ieee_set_status(original_fpe_status)