diff --git a/src/modules/module_interfaces.f90 b/src/modules/module_interfaces.f90 index c0905e2d8..cae1120a4 100644 --- a/src/modules/module_interfaces.f90 +++ b/src/modules/module_interfaces.f90 @@ -726,7 +726,7 @@ function symba_casedisruption (symba_plA, family, nmergeadd, mergeadd_list, x, v real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip real(DP), dimension(:), intent(inout) :: mass, radius, mass_res type(user_input_parameters),intent(inout) :: param - real(DP), intent(in) :: Qloss + real(DP), intent(inout) :: Qloss integer(I4B) :: status end function symba_casedisruption @@ -744,7 +744,7 @@ function symba_casehitandrun (symba_plA, family, nmergeadd, mergeadd_list, name, real(DP), dimension(:,:), intent(inout) :: x, v, Lspin, Ip real(DP), dimension(:), intent(inout) :: mass, radius, mass_res type(user_input_parameters),intent(inout) :: param - real(DP), intent(in) :: Qloss + real(DP), intent(inout) :: Qloss integer(I4B) :: status end function symba_casehitandrun @@ -776,7 +776,7 @@ function symba_casesupercatastrophic (symba_plA, family, nmergeadd, mergeadd_lis real(DP), dimension(:,:), intent(in) :: x, v, lspin, Ip real(DP), dimension(:), intent(in) :: mass, radius, mass_res type(user_input_parameters),intent(inout) :: param - real(DP), intent(in) :: Qloss + real(DP), intent(inout) :: Qloss integer(I4B) :: status end function symba_casesupercatastrophic end interface @@ -929,7 +929,7 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi type(user_input_parameters), intent(in) :: param type(symba_pl), intent(inout) :: symba_plA integer(I4B), dimension(:), intent(in) :: family - real(DP), intent(in) :: Qloss + real(DP), intent(inout) :: Qloss real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip real(DP), dimension(:), intent(inout) :: mass, radius, m_frag, rad_frag real(DP), dimension(:,:), intent(inout) :: Ip_frag diff --git a/src/modules/module_symba.f90 b/src/modules/module_symba.f90 index 30b29efc4..a93651dc3 100644 --- a/src/modules/module_symba.f90 +++ b/src/modules/module_symba.f90 @@ -103,78 +103,4 @@ MODULE module_symba end type symba_merger - type, public, extends(lambda_obj) :: symba_frag_lambda - procedure(abstract_objective_func), pointer, nopass :: ke_objective_func_ptr => null() - real(DP), dimension(:), allocatable :: m_frag - real(DP), dimension(:,:), allocatable :: x_frag - real(DP), dimension(:,:), allocatable :: v_frag - real(DP), dimension(NDIM) :: L_target - real(DP) :: ke_target - contains - generic :: init => ke_objective_func_init - procedure :: eval => ke_objective_func_eval - procedure, nopass :: ke_objective_func_init - final :: ke_objective_func_destroy - end type symba_frag_lambda - interface symba_frag_lambda - module procedure ke_objective_func_init - end interface - - abstract interface - function abstract_objective_func(vflat, v_frag, x_frag, m_frag, L_target, ke_target) result(fnorm) - ! Template for the kinetic energy constraint function used for minimizing - import DP - real(DP), dimension(:), intent(in) :: vflat !! Unrolled unknown velocity vectors - real(DP), dimension(:,:), intent(in) :: v_frag, x_frag !! Velocity and position vectors - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:), intent(in) :: L_target !! Target orbital momentum - real(DP), intent(in) :: ke_target !! Target kinetic energ - real(DP) :: fnorm !! The objective function result: norm of the vector composed of the tangential momentum and energy - end function - end interface - - contains - type(symba_frag_lambda) function ke_objective_func_init(lambda, v_frag, x_frag, m_frag, L_target, ke_target) - implicit none - ! Arguments - procedure(abstract_objective_func) :: lambda !! The lambda function - real(DP), dimension(:,:), intent(in) :: v_frag, x_frag !! Velocity and position vectors - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:), intent(in) :: L_target !! Target orbital momentum - real(DP), intent(in) :: ke_target !! Target kinetic energ - ! Internals - associate(self => ke_objective_func_init) - self%ke_objective_func_ptr => lambda - allocate(self%m_frag, source=m_frag) - allocate(self%v_frag, source=v_frag) - allocate(self%x_frag, source=x_frag) - self%L_target(:) = L_target(:) - self%ke_target = ke_target - end associate - return - end function ke_objective_func_init - - subroutine ke_objective_func_destroy(self) - implicit none - type(symba_frag_lambda) :: self - if (allocated(self%m_frag)) deallocate(self%m_frag) - if (allocated(self%v_frag)) deallocate(self%v_frag) - if (allocated(self%x_frag)) deallocate(self%x_frag) - if (associated(self%ke_objective_func_ptr)) nullify(self%ke_objective_func_ptr) - end subroutine ke_objective_func_destroy - - function ke_objective_func_eval(self, x) result(fnorm) - implicit none - ! Arguments - class(symba_frag_lambda), intent(in) :: self - real(DP), dimension(:), intent(in) :: x - ! Result - real(DP) :: fnorm - - if (associated(self%ke_objective_func_ptr)) then - fnorm = self%ke_objective_func_ptr(x, self%v_frag, self%x_frag, self%m_frag, self%L_target, self%ke_target) - else - error stop "KE Objective function was not initialized." - end if - end function ke_objective_func_eval END MODULE module_symba \ No newline at end of file diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index 21eddeca9..71c38b17f 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -1,3 +1,97 @@ +module symba_frag_pos_lambda_implementation + use swiftest_globals + use lambda_function + implicit none + + ! Define the class and interface used to implement the lambda function + type, public, extends(lambda_obj) :: ke_constraint + procedure(abstract_objective_func), pointer, nopass :: ke_objective_func_ptr => null() + real(DP), dimension(:), allocatable :: v_t_mag + real(DP), dimension(:,:), allocatable :: v_r_unit, v_t_unit + real(DP), dimension(:), allocatable :: m_frag + real(DP), dimension(:,:), allocatable :: x_frag + real(DP), dimension(NDIM) :: L_target + real(DP) :: ke_target + contains + generic :: init => ke_objective_func_init + procedure :: eval => ke_objective_func_eval + procedure, nopass :: ke_objective_func_init + final :: ke_objective_func_destroy + end type ke_constraint + interface ke_constraint + module procedure ke_objective_func_init + end interface + + abstract interface + function abstract_objective_func(v_r_mag, v_r_unit, v_t_mag, v_t_unit, x_frag, m_frag, L_target, ke_target) result(fnorm) + ! Template for the kinetic energy constraint function used for minimizing + import DP + real(DP), dimension(:), intent(in) :: v_r_mag !! Radial velocity mangitude + real(DP), dimension(:,:), intent(in) :: v_r_unit !! Radial velocity unit vector + real(DP), dimension(:), intent(in) :: v_t_mag !! Tangential velocity magnitude + real(DP), dimension(:,:), intent(in) :: v_t_unit !! Tangential velocity unit vector + real(DP), dimension(:,:), intent(in) :: x_frag !! Position vectors + real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses + real(DP), dimension(:), intent(in) :: L_target !! Target orbital momentum + real(DP), intent(in) :: ke_target !! Target kinetic energ + real(DP) :: fnorm !! The objective function result: norm of the vector composed of the tangential momentum and energy + end function + end interface + + contains + type(ke_constraint) function ke_objective_func_init(lambda, v_r_unit, v_t_mag, v_t_unit, x_frag, m_frag, L_target, ke_target) + implicit none + ! Arguments + procedure(abstract_objective_func) :: lambda !! The lambda function + real(DP), dimension(:,:), intent(in) :: v_r_unit !! Radial velocity unit vector + real(DP), dimension(:), intent(in) :: v_t_mag !! Tangential velocity magnitude + real(DP), dimension(:,:), intent(in) :: v_t_unit !! Tangential velocity unit vector + real(DP), dimension(:,:), intent(in) :: x_frag !! Position vectors + real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses + real(DP), dimension(:), intent(in) :: L_target !! Target orbital momentum + real(DP), intent(in) :: ke_target !! Target kinetic energ + ! Internals + associate(self => ke_objective_func_init) + self%ke_objective_func_ptr => lambda + allocate(self%v_r_unit, source=v_r_unit) + allocate(self%v_t_mag, source=v_t_mag) + allocate(self%v_t_unit, source=v_t_unit) + allocate(self%m_frag, source=m_frag) + allocate(self%x_frag, source=x_frag) + self%L_target(:) = L_target(:) + self%ke_target = ke_target + end associate + return + end function ke_objective_func_init + + subroutine ke_objective_func_destroy(self) + implicit none + type(ke_constraint) :: self + if (allocated(self%v_r_unit)) deallocate(self%v_r_unit) + if (allocated(self%v_t_mag)) deallocate(self%v_t_mag) + if (allocated(self%v_t_unit)) deallocate(self%v_t_unit) + if (allocated(self%x_frag)) deallocate(self%x_frag) + if (allocated(self%m_frag)) deallocate(self%m_frag) + if (associated(self%ke_objective_func_ptr)) nullify(self%ke_objective_func_ptr) + end subroutine ke_objective_func_destroy + + function ke_objective_func_eval(self, x) result(fval) + implicit none + ! Arguments + class(ke_constraint), intent(in) :: self + real(DP), dimension(:), intent(in) :: x + ! Result + real(DP) :: fval + + if (associated(self%ke_objective_func_ptr)) then + fval = self%ke_objective_func_ptr(x, self%v_r_unit, self%v_t_mag, self%v_t_unit, self%x_frag, self%m_frag, self%L_target, self%ke_target) + else + error stop "KE Objective function was not initialized." + end if + end function ke_objective_func_eval + +end module symba_frag_pos_lambda_implementation + subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, radius, & Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, lmerge, Qloss) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton @@ -9,37 +103,118 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad use module_helio use module_symba use module_swiftestalloc - use module_interfaces, EXCEPT_THIS_ONE => symba_frag_pos + use module_interfaces, except_this_one => symba_frag_pos implicit none ! Arguments type(user_input_parameters), intent(in) :: param type(symba_pl), intent(inout) :: symba_plA integer(I4B), dimension(:), intent(in) :: family - real(DP), intent(in) :: Qloss + real(DP), intent(inout) :: Qloss real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip real(DP), dimension(:), intent(inout) :: mass, radius, m_frag, rad_frag real(DP), dimension(:,:), intent(inout) :: Ip_frag real(DP), dimension(:,:), intent(inout) :: xb_frag, vb_frag, rot_frag logical, intent(out) :: lmerge ! Answers the question: Should this have been a merger instead? ! Internals - real(DP) :: mscale, rscale, vscale, Lscale, tscale , Qloss_scaled ! Scale factors that reduce quantities to O(~1) in the collisional system + real(DP), parameter :: f_spin = 0.20_DP !! Fraction of pre-impact orbital angular momentum that is converted to fragment spin + real(DP) :: mscale, rscale, vscale, Lscale, tscale ! Scale factors that reduce quantities to O(~1) in the collisional system + integer(I4B) :: i, j, nfrag, fam_size + logical, dimension(:), allocatable :: lexclude real(DP), dimension(NDIM, 2) :: rot, L_orb - integer(I4B) :: i, j, nfrag, fam_size, istart - real(DP), dimension(:,:), allocatable :: x_frag, v_frag - real(DP), dimension(NDIM) :: xcom, vcom, Ltot_before, Ltot_after, L_residual, L_spin_frag, L_target, L_frag + real(DP), dimension(:,:), allocatable :: x_frag, v_frag, v_r_unit, v_t_unit + real(DP), dimension(:), allocatable :: rmag, v_r_mag, v_t_mag + real(DP), dimension(NDIM) :: xcom, vcom, Ltot_before, Ltot_after, L_residual + real(DP), dimension(NDIM) :: L_frag_spin, L_frag_tot, L_frag_orb real(DP) :: mtot, Lmag_before, Lmag_after - real(DP) :: Etot_before, Etot_after, ke_before, pe_before - real(DP) :: pe_after, ke_spin_before, ke_spin_after, ke_after, ke_family, ke_target, ke_frag - real(DP), dimension(NDIM) :: h, dx, xc, vc, vcom_scaled - real(DP) :: rmag - logical, dimension(:), allocatable :: lexclude + real(DP) :: Etot_before, Etot_after, ke_orb_before, pe_before + real(DP) :: pe_after, ke_spin_before, ke_spin_after, ke_orb_after, ke_family, ke_target, ke_frag + real(DP), dimension(NDIM) :: x_col_unit, y_col_unit, z_col_unit character(len=*), parameter :: fmtlabel = "(A14,10(ES9.2,1X,:))" - real(DP), dimension(NDIM) :: x_cross_v, v_t_unit, h_unit, v_r_unit - real(DP), dimension(:,:), allocatable :: v_r, v_t + + allocate(x_frag, source=xb_frag) + allocate(v_frag, source=vb_frag) + nfrag = size(m_frag) + + associate(npl => symba_plA%helio%swiftest%nbody, status => symba_plA%helio%swiftest%status) + allocate(lexclude(npl)) + where (status(1:npl) == INACTIVE) ! Safety check in case one of the included bodies has been previously deactivated + lexclude(1:npl) = .true. + elsewhere + lexclude(1:npl) = .false. + end where + end associate + + call set_scale_factors() + call define_coordinate_system() + + call calculate_system_energy(ke_orb_before, ke_spin_before, pe_before, Ltot_before, linclude_fragments=.false.) + Lmag_before = norm2(Ltot_before(:)) + Etot_before = ke_orb_before + ke_spin_before + pe_before + + call define_pre_collisional_family() + call set_fragment_position_vectors() + + write(*, "(' ---------------------------------------------------------------------------')") + write(*, "(' Energy normalized by |Etot_before|')") + write(*, "(' | T_orb T_spin T pe Etot Ltot')") + write(*, "(' ---------------------------------------------------------------------------')") + write(*, "(' ---------------------------------------------------------------------------')") + write(*, "(' First pass to get angular momentum ')") + write(*, "(' ---------------------------------------------------------------------------')") + + call set_fragment_tangential_velocities() - associate(xhpl => symba_plA%helio%swiftest%xh, xbpl => symba_plA%helio%swiftest%xh, vbpl => symba_plA%helio%swiftest%vb, & - Mpl => symba_plA%helio%swiftest%mass, Ippl => symba_plA%helio%swiftest%Ip, radpl => symba_plA%helio%swiftest%radius, & - rotpl => symba_plA%helio%swiftest%rot, status => symba_plA%helio%swiftest%status, npl => symba_plA%helio%swiftest%nbody, name => symba_plA%helio%swiftest%id) + call calculate_system_energy(ke_orb_after, ke_spin_after, pe_after, Ltot_after, linclude_fragments=.true.) + Etot_after = ke_orb_after + ke_spin_after + pe_after + Lmag_after = norm2(Ltot_after(:)) + + write(*,fmtlabel) ' change |',(ke_orb_after - ke_orb_before) / abs(Etot_before), & + (ke_spin_after - ke_spin_before)/ abs(Etot_before), & + (ke_orb_after + ke_spin_after - ke_orb_before - ke_spin_before)/ abs(Etot_before), & + (pe_after - pe_before) / abs(Etot_before), & + (Etot_after - Etot_before) / abs(Etot_before), & + norm2(Ltot_after - Ltot_before) / Lmag_before + + call set_fragment_radial_velocities(lmerge) + write(*, "(' ---------------------------------------------------------------------------')") + write(*, "(' Second pass to get energy ')") + write(*, "(' ---------------------------------------------------------------------------')") + write(*,fmtlabel) ' Qloss |',-Qloss / abs(Etot_before) + write(*, "(' ---------------------------------------------------------------------------')") + write(*,fmtlabel) ' T_family |',ke_family / abs(Etot_before) + write(*,fmtlabel) ' T_frag targ |',ke_target / abs(Etot_before) + write(*, "(' ---------------------------------------------------------------------------')") + write(*,fmtlabel) ' T_frag calc |',ke_frag / abs(Etot_before) + write(*,fmtlabel) ' residual |',(ke_frag - ke_target) / abs(Etot_before) + + call calculate_system_energy(ke_orb_after, ke_spin_after, pe_after, Ltot_after, linclude_fragments=.true.) + Etot_after = ke_orb_after + ke_spin_after + pe_after + Lmag_after = norm2(Ltot_after(:)) + + write(*, "(' ---------------------------------------------------------------------------')") + write(*,fmtlabel) ' change |',(ke_orb_after - ke_orb_before) / abs(Etot_before), & + (ke_spin_after - ke_spin_before)/ abs(Etot_before), & + (ke_orb_after + ke_spin_after - ke_orb_before - ke_spin_before)/ abs(Etot_before), & + (pe_after - pe_before) / abs(Etot_before), & + (Etot_after - Etot_before) / abs(Etot_before), & + norm2(Ltot_after - Ltot_before) / Lmag_before + + lmerge = lmerge .or. ((Etot_after - Etot_before) / abs(Etot_before) > 0._DP) + + call restore_scale_factors() + + return + + contains + + ! Because of the complexity of this procedure, we have chosen to break it up into a series of nested subroutines. + + subroutine set_scale_factors() + !! author: David A. Minton + !! + !! Scales dimenional quantities to ~O(1) with respect to the collisional system. This scaling makes it easier for the non-linear minimization + !! to converge on a solution + implicit none mtot = 1.0_DP @@ -56,23 +231,43 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad v = v / vscale L_spin = L_spin / Lscale - xb_frag = xb_frag / rscale - vb_frag = vb_frag / rscale m_frag = m_frag / mscale - rot_frag = rot_frag * tscale rad_frag = rad_frag / rscale - Qloss_scaled = Qloss / (mscale * vscale**2) + Qloss = Qloss / (mscale * vscale**2) + return + end subroutine set_scale_factors - allocate(x_frag, source=xb_frag) - allocate(v_frag, source=vb_frag) - allocate(v_r, mold=v_frag) - allocate(v_t, mold=v_frag) + subroutine restore_scale_factors() + !! author: David A. Minton + !! + !! Restores dimenional quantities back to the system units + implicit none + mass = mass * mscale + radius = radius * rscale + x = x * rscale + v = v * vscale + L_spin = L_spin * Lscale + + xb_frag = xb_frag * rscale + vb_frag = vb_frag * vscale + m_frag = m_frag * mscale + rot_frag = rot_frag / tscale + rad_frag = rad_frag * rscale + return + end subroutine restore_scale_factors - fam_size = size(family) + subroutine define_coordinate_system() + !! author: David A. Minton + !! + !! Defines the collisional coordinate system, including the unit vectors of both the system and individual fragments. + implicit none + integer(I4B) :: i + real(DP), dimension(NDIM) :: x_cross_v, xc, vc, delta_r, delta_v + real(DP) :: r_col_norm, v_col_norm ! Find the center of mass of the collisional system - xcom(:) = mass(1) * x(:,1) + mass(2) * x(:,2) - vcom(:) = mass(1) * v(:,1) + mass(2) * v(:,2) + xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot + vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot L_orb(:, :) = 0.0_DP ! Compute orbital angular momentum of pre-impact system @@ -83,254 +278,257 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad L_orb(:, j) = mass(j) * x_cross_v(:) end do - allocate(lexclude(npl)) - where (status(1:npl) == INACTIVE) ! Safety check in case one of the included bodies has been previously deactivated - lexclude(1:npl) = .true. - elsewhere - lexclude(1:npl) = .false. - end where + ! Compute orbital angular momentum of pre-impact system. This will be the normal vector to the collision fragment plane + L_frag_tot(:) = L_spin(:,1) + L_spin(:,2) + L_orb(:, 1) + L_orb(:, 2) + L_frag_spin(:) = L_spin(:,1) + L_spin(:, 2) + f_spin * (L_orb(:, 1) + L_orb(:, 2)) + L_frag_orb(:) = L_frag_tot - L_frag_spin - call symba_frag_pos_energy_calc(npl, symba_plA, lexclude, ke_before, ke_spin_before, pe_before, Ltot_before) - Lmag_before = norm2(Ltot_before(:)) - Etot_before = ke_before + ke_spin_before + pe_before + delta_v(:) = v(:, 2) - v(:, 1) + v_col_norm = norm2(delta_v(:)) + delta_r(:) = x(:, 2) - x(:, 1) + r_col_norm = norm2(delta_r(:)) - ! 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)) / mscale * dot_product(vbpl(:,family(i)), vbpl(:,family(i))) / vscale**2 - lexclude(family(i)) = .true. ! For all subsequent energy calculations the pre-impact family members will be replaced by the fragments - end do - ke_family = 0.5_DP * ke_family + ! We will initialize fragments on a plane defined by the pre-impact system, with the z-axis aligned with the angular momentum vector + ! and the y-axis aligned with the pre-impact distance vector. + y_col_unit(:) = delta_r(:) / r_col_norm + z_col_unit(:) = L_frag_tot(:) / norm2(L_frag_tot) + ! The cross product of the z- by x-axis will give us the y-axis + call util_crossproduct(y_col_unit, z_col_unit, x_col_unit) - nfrag = size(m_frag) - ! Initialize positions and velocities of fragments that conserve angular momentum - call symba_frag_pos_initialize_fragments(x, v, L_orb, L_spin, mass, radius, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag, lmerge) - if (lmerge) return + return + end subroutine define_coordinate_system - ! Energy calculation requires the fragments to be in the system barcyentric frame, so we need to temporarily shift them - do i = 1, nfrag - xb_frag(:,i) = x_frag(:,i) + xcom(:) - vb_frag(:,i) = v_frag(:,i) + vcom(:) - end do - - call symba_frag_pos_energy_calc(npl, symba_plA, lexclude, ke_after, ke_spin_after, pe_after, Ltot_after, & - nfrag=nfrag, Ip_frag=Ip_frag, m_frag=m_frag, rad_frag=rad_frag, xb_frag=xb_frag, vb_frag=vb_frag, rot_frag=rot_frag) - Etot_after = ke_after + ke_spin_after + pe_after - Lmag_after = norm2(Ltot_after(:)) - - write(*, "(' ---------------------------------------------------------------------------')") - write(*, "(' Energy normalized by |Etot_before|')") - write(*, "(' | T_orb T_spin T pe Etot Ltot')") - write(*, "(' ---------------------------------------------------------------------------')") - write(*, "(' ---------------------------------------------------------------------------')") - write(*, "(' First pass to get angular momentum ')") - write(*, "(' ---------------------------------------------------------------------------')") - write(*,fmtlabel) ' change |',(ke_after - ke_before) / abs(Etot_before), & - (ke_spin_after - ke_spin_before)/ abs(Etot_before), & - (ke_after + ke_spin_after - ke_before - ke_spin_before)/ abs(Etot_before), & - (pe_after - pe_before) / abs(Etot_before), & - (Etot_after - Etot_before) / abs(Etot_before), & - norm2(Ltot_after - Ltot_before) / Lmag_before - write(*, "(' ---------------------------------------------------------------------------')") - write(*, "(' Second pass to get energy ')") - write(*, "(' ---------------------------------------------------------------------------')") - write(*,fmtlabel) ' Q_loss |',-Qloss_scaled / abs(Etot_before) - write(*, "(' ---------------------------------------------------------------------------')") + subroutine define_pre_collisional_family() + !! author: David A. Minton + !! + !! 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. + associate(vbpl => symba_plA%helio%swiftest%vb, Mpl => symba_plA%helio%swiftest%mass) + 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))) + lexclude(family(i)) = .true. ! For all subsequent energy calculations the pre-impact family members will be replaced by the fragments + end do + ke_family = 0.5_DP * ke_family / (mscale * vscale**2) + end associate + return + end subroutine define_pre_collisional_family - - ! Set the "target" ke_after (the value of the orbital kinetic energy that the fragments ought to have) - ke_target = ke_family + (ke_spin_before - ke_spin_after) + (pe_before - pe_after) - Qloss_scaled - L_target(:) = 0.0_DP - do i = 1, nfrag - call util_crossproduct(x_frag(:, i), v_frag(:, i), L_frag(:)) - L_target(:) = L_target(:) + L_frag(:) - end do - call symba_frag_pos_kinetic_energy(m_frag, ke_target, L_target, x_frag, v_frag, lmerge) - - write(*, "(' ---------------------------------------------------------------------------')") - write(*,fmtlabel) ' T_family |',ke_family / abs(Etot_before) - write(*,fmtlabel) ' T_frag targ |',ke_target / abs(Etot_before) + subroutine calculate_system_energy(ke_orb, ke_spin, pe, Ltot, linclude_fragments) + !! Author: David A. Minton + !! + !! Calculates total system energy, including all bodies in the symba_plA list that do not have a corresponding value of the lexclude array that is true + !! and optionally including fragments. + use module_swiftestalloc + implicit none + ! Arguments + real(DP), intent(out) :: ke_orb, ke_spin, pe + real(DP), dimension(:), intent(out) :: Ltot + logical, intent(in) :: linclude_fragments + ! Internals + integer(I4B) :: i, npl_new, nplm + logical, dimension(:), allocatable :: ltmp + logical :: lk_plpl + real(DP) :: te + type(symba_pl) :: symba_plwksp - ! Shift the fragments into the system barycenter frame - do i = 1, nfrag - xb_frag(:,i) = x_frag(:, i) + xcom(:) - vb_frag(:,i) = v_frag(:, i) + vcom(:) - end do + ! Because we're making a copy of symba_pl with the excludes/fragments appended, we need to deallocate the + ! big k_plpl array and recreate it when we're done, otherwise we run the risk of blowing up the memory by + ! allocating two of these ginormous arrays simulteouously. This is not particularly efficient, but as this + ! subroutine should be called relatively infrequently, it shouldn't matter too much. + !if (allocated(symba_plA%helio%swiftest%k_plpl)) deallocate(symba_plA%helio%swiftest%k_plpl) - ke_frag = 0._DP - do i = 1, nfrag - ke_frag = ke_frag + 0.5_DP * m_frag(i) * dot_product(vb_frag(:, i), vb_frag(:, i)) - end do - write(*, "(' ---------------------------------------------------------------------------')") - write(*,fmtlabel) ' T_frag calc |',ke_frag / abs(Etot_before) - write(*,fmtlabel) ' residual |',(ke_frag - ke_target) / abs(Etot_before) - - ! Calculate the new energy of the system of fragments - call symba_frag_pos_energy_calc(npl, symba_plA, lexclude, ke_after, ke_spin_after, pe_after, Ltot_after,& - nfrag=nfrag, Ip_frag=Ip_frag, m_frag=m_frag, rad_frag=rad_frag, xb_frag=xb_frag, vb_frag=vb_frag, rot_frag=rot_frag) - Etot_after = ke_after + ke_spin_after + pe_after - Lmag_after = norm2(Ltot_after(:)) - - write(*, "(' ---------------------------------------------------------------------------')") - write(*,fmtlabel) ' change |',(ke_after - ke_before) / abs(Etot_before), & - (ke_spin_after - ke_spin_before)/ abs(Etot_before), & - (ke_after + ke_spin_after - ke_before - ke_spin_before)/ abs(Etot_before), & - (pe_after - pe_before) / abs(Etot_before), & - (Etot_after - Etot_before) / abs(Etot_before), & - norm2(Ltot_after - Ltot_before) / Lmag_before + ! Build the internal planet list out of the non-excluded bodies and optionally with fragments appended. This + ! will get passed to the energy calculation subroutine so that energy is computed exactly the same way is it + ! is in the main program. + associate(npl => symba_plA%helio%swiftest%nbody) + lk_plpl = allocated(symba_plA%helio%swiftest%k_plpl) + if (lk_plpl) deallocate(symba_plA%helio%swiftest%k_plpl) + if (linclude_fragments) then ! Temporarily expand the planet list to feed it into symba_energy + npl_new = npl + nfrag + else + npl_new = npl + end if + call symba_pl_allocate(symba_plwksp, npl_new) + + ! Copy over old data + symba_plwksp%helio%swiftest%id(1:npl) = symba_plA%helio%swiftest%id(1:npl) + symba_plwksp%helio%swiftest%status(1:npl) = symba_plA%helio%swiftest%status(1:npl) + symba_plwksp%helio%swiftest%mass(1:npl) = symba_plA%helio%swiftest%mass(1:npl) / mscale + symba_plwksp%helio%swiftest%radius(1:npl) = symba_plA%helio%swiftest%radius(1:npl) / rscale + symba_plwksp%helio%swiftest%xh(:,1:npl) = symba_plA%helio%swiftest%xh(:,1:npl) / rscale + symba_plwksp%helio%swiftest%vh(:,1:npl) = symba_plA%helio%swiftest%vh(:,1:npl) / vscale + symba_plwksp%helio%swiftest%rhill(1:npl) = symba_plA%helio%swiftest%rhill(1:npl) / rscale + symba_plwksp%helio%swiftest%xb(:,1:npl) = symba_plA%helio%swiftest%xb(:,1:npl) / rscale + symba_plwksp%helio%swiftest%vb(:,1:npl) = symba_plA%helio%swiftest%vb(:,1:npl) / vscale + symba_plwksp%helio%swiftest%rot(:,1:npl) = symba_plA%helio%swiftest%rot(:,1:npl) * tscale + symba_plwksp%helio%swiftest%Ip(:,1:npl) = symba_plA%helio%swiftest%Ip(:,1:npl) + + if (linclude_fragments) then ! Append the fragments if they are included + ! Energy calculation requires the fragments to be in the system barcyentric frame, s + symba_plwksp%helio%swiftest%Ip(:,npl+1:npl_new) = Ip_frag(:,:) + symba_plwksp%helio%swiftest%mass(npl+1:npl_new) = m_frag(:) + symba_plwksp%helio%swiftest%radius(npl+1:npl_new) = rad_frag(:) + symba_plwksp%helio%swiftest%xb(:,npl+1:npl_new) = xb_frag(:,:) + symba_plwksp%helio%swiftest%vb(:,npl+1:npl_new) = vb_frag(:,:) + symba_plwksp%helio%swiftest%rot(:,npl+1:npl_new) = rot_frag(:,:) + symba_plwksp%helio%swiftest%status(npl+1:npl_new) = COLLISION + call coord_b2h(npl_new, symba_plwksp%helio%swiftest) + allocate(ltmp(npl_new)) + ltmp(1:npl) = lexclude(1:npl) + ltmp(npl+1:npl_new) = .false. + call move_alloc(ltmp, lexclude) + + ke_frag = 0._DP + do i = 1, nfrag + ke_frag = ke_frag + 0.5_DP * m_frag(i) * dot_product(vb_frag(:, i), vb_frag(:, i)) + end do + end if + + where (lexclude(1:npl)) + symba_plwksp%helio%swiftest%status(1:npl) = INACTIVE + end where - lmerge = lmerge .or. ((Etot_after - Etot_before) / abs(Etot_before) > 0._DP) + nplm = count(symba_plwksp%helio%swiftest%mass > param%mtiny) + call util_dist_index_plpl(npl_new, nplm, symba_plwksp) + call symba_energy_eucl(npl_new, symba_plwksp, param%j2rp2, param%j4rp4, ke_orb, ke_spin, pe, te, Ltot) + + ! Restore the big array + deallocate(symba_plwksp%helio%swiftest%k_plpl) + nplm = count(symba_plA%helio%swiftest%mass > param%mtiny) + if (lk_plpl) call util_dist_index_plpl(npl, nplm, symba_plA) + end associate + return + end subroutine calculate_system_energy + + subroutine shift_vector_to_origin(m_frag, vec_frag) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Adjusts the position or velocity of the fragments as needed to align them with the origin + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses + real(DP), dimension(:,:), intent(inout) :: vec_frag !! Fragment positions or velocities in the center of mass frame - mass = mass * mscale - radius = radius * rscale - x = x * rscale - v = v * vscale - L_spin = L_spin * Lscale + ! Internals + real(DP), dimension(NDIM) :: mvec_frag, COM_offset + integer(I4B) :: i, nfrag + real(DP) :: mtot - xb_frag = xb_frag * rscale - vb_frag = vb_frag * rscale - m_frag = m_frag * mscale - rot_frag = rot_frag / tscale - rad_frag = rad_frag * rscale + nfrag = size(m_frag) + mtot = sum(m_frag) + mvec_frag(:) = 0.0_DP - end associate - return + do i = 1, nfrag + mvec_frag = mvec_frag(:) + vec_frag(:,i) * m_frag(i) + end do + COM_offset(:) = -mvec_frag(:) / mtot + do i = 1, nfrag + vec_frag(:, i) = vec_frag(:, i) + COM_offset(:) + end do - contains + return + end subroutine shift_vector_to_origin - subroutine symba_frag_pos_initialize_fragments(x, v, L_orb, L_spin, mass, radius, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag, lmerge) + subroutine set_fragment_position_vectors() !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Initializes the orbits of the fragments around the center of mass. The fragments are initially placed on a plane defined by the !! pre-impact angular momentum. They are distributed on an ellipse surrounding the center of mass. !! The initial positions do not conserve energy or momentum, so these need to be adjusted later. - implicit none - ! Arguments - real(DP), dimension(:,:), intent(in) :: x, v, L_orb, L_spin !! Pre-impact angular momentum vectors - real(DP), dimension(:), intent(in) :: mass, radius !! Pre-impact masses and radii - real(DP), dimension(:), intent(in) :: m_frag, rad_frag !! Fragment masses and radii - real(DP), dimension(:,:), intent(in) :: Ip_frag !! Fragment prinicpal moments of inertia - real(DP), dimension(:,:), intent(out) :: x_frag, v_frag, rot_frag !! Fragment position, velocities, and spin states - logical, intent(out) :: lmerge - ! Internals - real(DP) :: theta, v_frag_norm, r_frag_norm, v_col_norm, r_col_norm - real(DP), dimension(NDIM) :: Ltot, delta_r, delta_v - real(DP), dimension(NDIM) :: x_col_unit, y_col_unit, z_col_unit - integer(I4B) :: i - - ! Now create the fragment distribution - ! Compute orbital angular momentum of pre-impact system. This will be the normal vector to the collision fragment plane - Ltot = L_spin(:,1) + L_spin(:,2) + L_orb(:, 1) + L_orb(:, 2) - - delta_v(:) = v(:, 2) - v(:, 1) - v_col_norm = norm2(delta_v(:)) - delta_r(:) = x(:, 2) - x(:, 1) - r_col_norm = norm2(delta_r(:)) - - ! We will initialize fragments on a plane defined by the pre-impact system, with the y-axis aligned with the angular momentum vector - ! and the z-axis aligned with the pre-impact distance vector. - y_col_unit = Ltot(:) / norm2(Ltot(:)) - z_col_unit(:) = delta_r(:) / r_col_norm - ! The cross product of the z- by x-axis will give us the y-axis - call util_crossproduct(y_col_unit, z_col_unit, x_col_unit) + implicit none + real(DP) :: r_max + real(DP), dimension(NDIM) :: L, L_sigma - ! The angular spacing of fragments on the ellipse - We will only use half the ellipse - theta = 2 * PI / nfrag + allocate(rmag(nfrag)) + allocate(v_r_mag(nfrag)) + allocate(v_t_mag(nfrag)) + allocate(v_r_unit(NDIM,nfrag)) + allocate(v_t_unit(NDIM,nfrag)) ! Re-normalize position and velocity vectors by the fragment number so that for our initial guess we weight each ! fragment position by the mass and assume equipartition of energy for the velocity - r_col_norm = 1.5_DP * 2 * sum(rad_frag(:)) / theta / nfrag + r_max = 1.5_DP * sum(rad_frag(:)) / PI ! We will treat the first fragment of the list as a special case. x_frag(:, 1) = -z_col_unit(:) call random_number(x_frag(:,2:nfrag)) - x_frag(:, :) = x_frag(:, :) * r_col_norm - call symba_frag_pos_com_adjust(xcom, m_frag, x_frag) - v_frag(:,:) = 0._DP + x_frag(:, :) = x_frag(:, :) * r_max + call shift_vector_to_origin(m_frag, x_frag) - call symba_frag_pos_ang_mtm(L_orb, L_spin, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag, lmerge) + do i = 1, nfrag + xb_frag(:,i) = x_frag(:,i) + xcom(:) + rmag(i) = norm2(x_frag(:, i)) + v_r_unit(:, i) = x_frag(:, i) / rmag(i) + call random_number(L_sigma(:)) ! Randomize the tangential velocity direction. This helps to ensure that the tangential velocity doesn't completely line up with the angular momentum vector, + ! otherwise we can get an ill-conditioned system + L(:) = z_col_unit(:) + 2e-3_DP * (L_sigma(:) - 0.5_DP) + L(:) = L(:) / norm2(L(:)) + call util_crossproduct(L(:), v_r_unit(:, i), v_t_unit(:, i)) + end do return - end subroutine symba_frag_pos_initialize_fragments + end subroutine set_fragment_position_vectors - subroutine symba_frag_pos_ang_mtm(L_orb, L_spin, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag, lmerge) + subroutine set_fragment_tangential_velocities() !! 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 - real(DP), dimension(:,:), intent(in) :: L_orb, L_spin !! Pre-impact position, velocity, and spin states, Ip_frag - real(DP), dimension(:), intent(in) :: m_frag, rad_frag !! Fragment masses and radii - real(DP), dimension(:,:), intent(in) :: Ip_frag, x_frag !! Fragment prinicpal moments of inertia and position vectors - real(DP), dimension(:,:), intent(out) :: v_frag, rot_frag !! Fragment velocities and spin states - logical, intent(out) :: lmerge ! Internals integer(I4B) :: i real(DP) :: L_orb_mag - real(DP), dimension(NDIM) :: L_orb_old, L_spin_frag, L_spin_new - real(DP), parameter :: f_spin = 0.20_DP !! Fraction of pre-impact orbital angular momentum that is converted to fragment spin 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(:,:), allocatable :: v_t_unit, r_unit - real(DP), dimension(:), allocatable :: v_t_mag, rmag - real(DP), dimension(NDIM) :: L_lin_others, L_orb_others, L_unit, L_frag, L_sigma - - allocate(v_t_mag, mold=m_frag) - allocate(rmag, mold=m_frag) - allocate(r_unit, mold=v_frag) - allocate(v_t_unit, mold=v_frag) + real(DP), dimension(NDIM) :: L_lin_others, L_orb_others, L + + v_frag(:,:) = 0.0_DP - L_orb_old(:) = L_orb(:, 1) + L_orb(:, 2) - ! Divide up the pre-impact spin angular momentum equally between the various bodies by mass - L_spin_new(:) = L_spin(:,1) + L_spin(:, 2) + f_spin * L_orb_old(:) - L_spin_frag(:) = L_spin_new(:) / nfrag do i = 1, nfrag - rot_frag(:,i) = L_spin_frag(:) / (Ip_frag(:, i) * m_frag(i) * rad_frag(i)**2) + rot_frag(:,i) = L_frag_spin(:) / nfrag / (Ip_frag(:, i) * m_frag(i) * rad_frag(i)**2) end do - L_orb_old(:) = L_orb_old(:) * (1.0_DP - f_spin) - L_unit(:) = L_orb_old(:) / norm2(L_orb_old(:)) - L_orb_mag = norm2(L_orb_old(:)) + L_orb_mag = norm2(L_frag_orb(:)) ! 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 - rmag(i) = norm2(x_frag(:, i)) - r_unit(:, i) = x_frag(:, i) / rmag(i) - call random_number(L_sigma(:)) ! Randomize the tangential velocity direction. This helps to ensure that the tangential velocity doesn't completely line up with the angular momentum vector, - ! otherwise we can get an ill-conditioned system - L_frag(:) = L_unit(:) + 2e-3_DP * (L_sigma(:) - 0.5_DP) - L_frag(:) = L_frag(:) / norm2(L_frag(:)) - call util_crossproduct(L_frag(:), r_unit(:, i), v_t_unit(:, i)) 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) = m_frag(i) * v_t_unit(:, i) - call util_crossproduct(r_unit(:, i), v_t_unit(:, i), L_frag(:)) - A(4:6, i) = m_frag(i) * rmag(i) * L_frag(:) + call util_crossproduct(v_r_unit(:, i), v_t_unit(:, i), L(:)) + A(4:6, i) = m_frag(i) * rmag(i) * L(:) else !For the remining bodies, distribute the angular momentum equally amongs them - v_t_mag(i) = L_orb_mag / (m_frag(i) * rmag(i) * nfrag) + v_t_mag(i) = 0.25_DP * L_orb_mag / (m_frag(i) * rmag(i) * nfrag) v_frag(:, i) = v_t_mag(i) * v_t_unit(:, i) L_lin_others(:) = L_lin_others(:) + m_frag(i) * v_frag(:, i) - call util_crossproduct(x_frag(:, i), v_frag(:, i), L_frag(:)) - L_orb_others(:) = L_orb_others(:) + m_frag(i) * L_frag(:) + call util_crossproduct(x_frag(:, i), v_frag(:, i), L(:)) + L_orb_others(:) = L_orb_others(:) + m_frag(i) * L(:) end if end do b(1:3) = -L_lin_others(:) - b(4:6) = L_orb_old(:) - L_orb_others(:) + b(4:6) = L_frag_orb(:) - L_orb_others(:) v_t_mag(1:6) = util_solve_linear_system(A, b, 6, lmerge) if (lmerge) return do i = 1, 6 v_frag(:, i) = v_t_mag(i) * v_t_unit(:, i) end do + do i = 1, nfrag + vb_frag(:,i) = v_frag(:,i) + vcom(:) + end do + return - end subroutine symba_frag_pos_ang_mtm + end subroutine set_fragment_tangential_velocities - subroutine symba_frag_pos_kinetic_energy(m_frag, ke_target, L_target, x_frag, v_frag, lmerge) + subroutine set_fragment_radial_velocities(lmerge) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! @@ -340,194 +538,100 @@ subroutine symba_frag_pos_kinetic_energy(m_frag, ke_target, L_target, x_frag, v_ !! It takes in the initial "guess" of velocities and solve for the a scaling factor applied to the radial component wrt the !! center of mass frame needed to correct the kinetic energy of the fragments in the system barycenter frame to match that of !! the target kinetic energy required to satisfy the constraints. + use symba_frag_pos_lambda_implementation implicit none ! Arguments - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:,:), intent(in) :: x_frag !! Fragment position vectors - real(DP), intent(in) :: ke_target !! Target kinetic energy - real(DP), dimension(:), intent(in) :: L_target !! Target kinetic energy - real(DP), dimension(:,:), intent(inout) :: v_frag !! Fragment position and velocities in the center of mass frame logical, intent(out) :: lmerge - ! Internals - integer(I4B) :: i, neval - real(DP), parameter :: TOL = 1e-12_DP - real(DP), dimension(:), allocatable :: vflat + real(DP), parameter :: TOL = 1e-6_DP + real(DP), dimension(:), allocatable :: vflat logical :: lerr + integer(I4B) :: i + real(DP) :: ke_tangential + real(DP), dimension(:), allocatable :: v_r_initial + real(DP), dimension(:,:), allocatable :: v_r + + ! Set the "target" ke_orb_after (the value of the orbital kinetic energy that the fragments ought to have) + ke_tangential = 0.5_DP * sum(m_frag(:) * v_t_mag(:)**2) + ke_target = ke_family - ke_tangential + (ke_spin_before - ke_spin_after) + (pe_before - pe_after) - Qloss + if (ke_target < 0.0_DP) then + lmerge = .true. + return + end if + ! Initialize radial velocity magnitudes with a random value: + allocate(v_r_initial, source=v_r_mag) + call random_number(v_r_initial(:)) + v_r_initial(:) = v_r_initial(:) * sqrt(2 * ke_target / nfrag / m_frag(:)) ! Initialize the lambda function using a structure constructor that calls the init method ! Minimize error using the BFGS optimizer - vflat = util_minimize_bfgs(symba_frag_lambda(lambda=symba_frag_pos_ke_objective_function, v_frag=v_frag, x_frag=x_frag, m_frag=m_frag, L_target=L_target, ke_target=ke_target), & - NDIM*nfrag, reshape(v_frag,[NDIM*nfrag]), TOL, lerr) + v_r_mag(:) = util_minimize_bfgs(ke_constraint(ke_objective_function, v_r_unit, v_t_mag, v_t_unit, x_frag, m_frag, L_frag_orb, ke_target), nfrag, v_r_initial, TOL, lerr) if (lerr) then ! No solution exists for this case, so we need to indicate that this should be a merge ! This may happen due to setting the tangential velocities too high when setting the angular momentum constraint lmerge = .true. v_frag(:,:) = 0.0_DP else - v_frag(:,:) = reshape(vflat,shape(v_frag)) + lmerge = .false. + ! Shift the radial velocity vectors to align with the center of mass of the collisional system (the origin) + allocate(v_r, mold=v_frag) + do i = 1, nfrag + v_r(:, i) = v_r_mag(i) * v_r_unit(:, i) + end do + call shift_vector_to_origin(m_frag, v_r) + + ! Recombine the tangential and radial components into the final velocity vector + do i = 1, nfrag + v_r_mag(i) = dot_product(v_r(:,i), v_r_unit(:, i)) + v_frag(:, i) = v_r_mag(i) * v_r_unit(:, i) + v_t_mag(i) * v_t_unit(:, i) + end do + !call shift_vector_to_origin(m_frag, v_frag) end if + do i = 1, nfrag + vb_frag(:,i) = v_frag(:,i) + vcom(:) + end do + return - end subroutine symba_frag_pos_kinetic_energy + end subroutine set_fragment_radial_velocities - function symba_frag_pos_ke_objective_function(vflat, v_frag, x_frag, m_frag, L_target, ke_target) result(fnorm) + function ke_objective_function(v_r_mag, v_r_unit, v_t_mag, v_t_unit, x_frag, m_frag, L_target, ke_target) result(fval) ! 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) :: vflat !! Unrolled unknown velocity vectors - real(DP), dimension(:,:), intent(in) :: v_frag, x_frag !! Velocity and position vectors - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:), intent(in) :: L_target !! Target orbital momentum - real(DP), intent(in) :: ke_target !! Target kinetic energ + real(DP), dimension(:), intent(in) :: v_r_mag !! Unknown radial component of fragment velocity vector + real(DP), dimension(:), intent(in) :: v_t_mag !! Tangential component of velocity vector set previously by angular momentum constraint + real(DP), dimension(:,:), intent(in) :: v_r_unit, v_t_unit !! Radial and tangential unit vectors for each fragment + real(DP), dimension(:,:), intent(in) :: x_frag !! Velocity and position vectors + real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses + real(DP), dimension(:), intent(in) :: L_target !! Target orbital momentum + real(DP), intent(in) :: ke_target !! Target kinetic energ ! Result - real(DP) :: fnorm !! The objective function result: norm of the vector composed of the tangential momentum and energy - !! Minimizing this brings us closer to our objective + real(DP) :: fval !! The objective function result: norm of the vector composed of the tangential momentum and energy + !! Minimizing this brings us closer to our objective ! Internals - real(DP), dimension(7) :: f_vec - integer(I4B) :: i, nfrag, nsol - real(DP), dimension(NDIM) :: L - real(DP), dimension(:,:), allocatable :: v_sol + integer(I4B) :: i, nfrag, nsol + real(DP), dimension(NDIM) :: L + real(DP), dimension(:), allocatable :: v_r_mag_shift + real(DP), dimension(:,:), allocatable :: v_r nfrag = size(m_frag) - nsol = size(vflat) / NDIM - allocate(v_sol(NDIM,nsol)) - v_sol = reshape(vflat, shape(v_sol)) - - ! Linear momentum constraint - f_vec(1) = sum(m_frag(:) * [v_sol(1, 1:nsol), v_frag(1, nsol + 1:nfrag)]) - f_vec(2) = sum(m_frag(:) * [v_sol(2, 1:nsol), v_frag(2, nsol + 1:nfrag)]) - f_vec(3) = sum(m_frag(:) * [v_sol(3, 1:nsol), v_frag(3, nsol + 1:nfrag)]) - ! Angular momentum and kinetic energy constraints - f_vec(4:6) = -L_target(:) - f_vec(7) = -ke_target - do i = 1, nsol - call util_crossproduct(x_frag(:,i), v_sol(:, i), L(:)) - f_vec(4:6) = f_vec(4:6) + m_frag(i) * L(:) - f_vec(7) = f_vec(7) + 0.5_DP * m_frag(i) * dot_product(v_sol(:,i), v_sol(:, i)) - end do - do i = nsol + 1, nfrag - call util_crossproduct(x_frag(:,i), v_frag(:, i), L(:)) - f_vec(4:6) = f_vec(4:6) + m_frag(i) * L(:) - f_vec(7) = f_vec(7) + 0.5_DP * m_frag(i) * dot_product(v_frag(:,i), v_frag(:, i)) - end do - - fnorm = norm2(f_vec(:)) - - return - - end function symba_frag_pos_ke_objective_function - - subroutine symba_frag_pos_com_adjust(vec_com, m_frag, vec_frag) - !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Adjusts the position or velocity of the fragments as needed to align them with the original trajectory center of mass. - implicit none - ! Arguments - real(DP), dimension(:), intent(in) :: vec_com !! Center of mass position or velocity in the system barycenter frame - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:,:), intent(inout) :: vec_frag !! Fragment positions or velocities in the center of mass frame - - ! Internals - real(DP), dimension(NDIM) :: mvec_frag, COM_offset - integer(I4B) :: i - - mvec_frag(:) = 0.0_DP - + allocate(v_r, mold=v_r_unit) + ! In order to keep satisfying the kinetic energy constraint, we must shift the origin of the radial component of the velocities to the center of mass do i = 1, nfrag - mvec_frag = mvec_frag(:) + (vec_frag(:,i) + vec_com(:)) * m_frag(i) + v_r(:,i) = v_r_mag(i) * v_r_unit(:, i) end do - COM_offset(:) = vec_com(:) - mvec_frag(:) / mtot - do i = 1, nfrag - vec_frag(:, i) = vec_frag(:, i) + COM_offset(:) + call shift_vector_to_origin(m_frag, v_r) + allocate(v_r_mag_shift, mold=v_r_mag) + do i = 1, nfrag + v_r_mag_shift(i) = dot_product(v_r(:, i), v_r_unit(:, i)) end do - return - end subroutine symba_frag_pos_com_adjust - - subroutine symba_frag_pos_energy_calc(npl, symba_plA, lexclude, ke_orbit, ke_spin, pe, Ltot, nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag) - !! Author: David A. Minton - !! - !! Calculates total system energy, including all bodies in the symba_plA list that do not have a corresponding value of the lexclude array that is true - !! and optionally including fragments. - use module_swiftestalloc - implicit none - ! Arguments - integer(I4B), intent(inout) :: npl - type(symba_pl), intent(inout) :: symba_plA - logical, dimension(:), allocatable, intent(inout) :: lexclude - real(DP), intent(out) :: ke_orbit, ke_spin, pe - real(DP), dimension(:), intent(out) :: Ltot - integer(I4B), intent(in), optional :: nfrag - real(DP), dimension(:), intent(in), optional :: m_frag, rad_frag - real(DP), dimension(:,:), intent(in), optional :: Ip_frag, xb_frag, vb_frag, rot_frag - ! Internals - integer(I4B) :: i, npl_new, nplm - logical, dimension(:), allocatable :: ltmp - logical :: lk_plpl - real(DP) :: te - type(symba_pl) :: symba_plwksp - - ! Because we're making a copy of symba_pl with the excludes/fragments appended, we need to deallocate the - ! big k_plpl array and recreate it when we're done, otherwise we run the risk of blowing up the memory by - ! allocating two of these ginormous arrays simulteouously. This is not particularly efficient, but as this - ! subroutine should be called relatively infrequently, it shouldn't matter too much. - !if (allocated(symba_plA%helio%swiftest%k_plpl)) deallocate(symba_plA%helio%swiftest%k_plpl) - - ! Build the internal planet list out of the non-excluded bodies and optionally with fragments appended. This - ! will get passed to the energy calculation subroutine so that energy is computed exactly the same way is it - ! is in the main program. - lk_plpl = allocated(symba_plA%helio%swiftest%k_plpl) - if (lk_plpl) deallocate(symba_plA%helio%swiftest%k_plpl) - if (present(nfrag)) then ! Temporarily expand the planet list to feed it into symba_energy - npl_new = npl + nfrag - else - npl_new = npl - end if - call symba_pl_allocate(symba_plwksp, npl_new) - - ! Copy over old data - symba_plwksp%helio%swiftest%id(1:npl) = symba_plA%helio%swiftest%id(1:npl) - symba_plwksp%helio%swiftest%status(1:npl) = symba_plA%helio%swiftest%status(1:npl) - symba_plwksp%helio%swiftest%mass(1:npl) = symba_plA%helio%swiftest%mass(1:npl) / mscale - symba_plwksp%helio%swiftest%radius(1:npl) = symba_plA%helio%swiftest%radius(1:npl) / rscale - symba_plwksp%helio%swiftest%xh(:,1:npl) = symba_plA%helio%swiftest%xh(:,1:npl) / rscale - symba_plwksp%helio%swiftest%vh(:,1:npl) = symba_plA%helio%swiftest%vh(:,1:npl) / vscale - symba_plwksp%helio%swiftest%rhill(1:npl) = symba_plA%helio%swiftest%rhill(1:npl) / rscale - symba_plwksp%helio%swiftest%xb(:,1:npl) = symba_plA%helio%swiftest%xb(:,1:npl) / rscale - symba_plwksp%helio%swiftest%vb(:,1:npl) = symba_plA%helio%swiftest%vb(:,1:npl) / vscale - symba_plwksp%helio%swiftest%rot(:,1:npl) = symba_plA%helio%swiftest%rot(:,1:npl) * tscale - symba_plwksp%helio%swiftest%Ip(:,1:npl) = symba_plA%helio%swiftest%Ip(:,1:npl) - - if (present(nfrag)) then ! Append the fragments if they are included - symba_plwksp%helio%swiftest%Ip(:,npl+1:npl_new) = Ip_frag(:,:) - symba_plwksp%helio%swiftest%mass(npl+1:npl_new) = m_frag(:) - symba_plwksp%helio%swiftest%radius(npl+1:npl_new) = rad_frag(:) - symba_plwksp%helio%swiftest%xb(:,npl+1:npl_new) = xb_frag(:,:) - symba_plwksp%helio%swiftest%vb(:,npl+1:npl_new) = vb_frag(:,:) - symba_plwksp%helio%swiftest%rot(:,npl+1:npl_new) = rot_frag(:,:) - symba_plwksp%helio%swiftest%status(npl+1:npl_new) = COLLISION - call coord_b2h(npl_new, symba_plwksp%helio%swiftest) - allocate(ltmp(npl_new)) - ltmp(1:npl) = lexclude(1:npl) - ltmp(npl+1:npl_new) = .false. - call move_alloc(ltmp, lexclude) - end if + fval = 0.5_DP * sum(m_frag(:) * v_r_mag_shift(:)**2) - ke_target - where (lexclude(1:npl)) - symba_plwksp%helio%swiftest%status(1:npl) = INACTIVE - end where - - nplm = count(symba_plwksp%helio%swiftest%mass > param%mtiny) - call util_dist_index_plpl(npl_new, nplm, symba_plwksp) - call symba_energy_eucl(npl_new, symba_plwksp, param%j2rp2, param%j4rp4, ke_orbit, ke_spin, pe, te, Ltot) + return - ! Restore the big array - deallocate(symba_plwksp%helio%swiftest%k_plpl) - nplm = count(symba_plA%helio%swiftest%mass > param%mtiny) - if (lk_plpl) call util_dist_index_plpl(npl, nplm, symba_plA) + end function ke_objective_function - return - end subroutine symba_frag_pos_energy_calc end subroutine symba_frag_pos diff --git a/src/util/util_minimize_bfgs.f90 b/src/util/util_minimize_bfgs.f90 index 08204427b..fe942109d 100644 --- a/src/util/util_minimize_bfgs.f90 +++ b/src/util/util_minimize_bfgs.f90 @@ -1,4 +1,4 @@ -function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) +function util_minimize_bfgs(f, N, x0, eps1, lerr) result(x1) !! author: David A. Minton !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !! This function implements the Broyden-Fletcher-Goldfarb-Shanno method to determine the minimum of a function of N variables. @@ -21,14 +21,16 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) integer(I4B), intent(in) :: N class(lambda_obj), intent(in) :: f real(DP), dimension(:), intent(in) :: x0 - real(DP), intent(in) :: eps + real(DP), intent(in) :: eps1 logical, intent(out) :: lerr ! Result real(DP), dimension(:), allocatable :: x1 ! Internals integer(I4B) :: i, j, k, l, conv, num, fnum - integer(I4B), parameter :: MAXLOOP = 200 !! Maximum number of loops before method is determined to have failed - real(DP), parameter :: gradeps = 1e-5_DP !! Tolerance for gradient calculations + !integer(I4B), parameter :: MAXLOOP = 10 !! Maximum number of loops before method is determined to have failed + integer(I4B) :: MAXLOOP + !real(DP), parameter :: gradeps = 1e-5_DP !! Tolerance for gradient calculations + real(DP) :: gradeps, eps !! Tolerance for gradient calculations real(DP), dimension(N) :: S !! Direction vectors real(DP), dimension(N) :: Snorm !! normalized direction real(DP), dimension(N,N) :: H !! Approximated inverse Hessian matrix @@ -39,6 +41,11 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) real(DP), dimension(N,N) :: PP, PyH, HyP real(DP) :: yHy, Py + open(unit=42,file='input.txt',status='old') + read(42,*) MAXLOOP + read(42,*) gradeps + read(42,*) eps + close(42) fnum = 0 lerr = .false. ! Initialize approximate Hessian with the identity matrix (i.e. begin with method of steepest descent) @@ -61,7 +68,7 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) end if ! normalize gradient Snorm(:) = S(:) / norm2(S) - num = fnum + minimize1D(f, x1, Snorm, N, eps, astar) + num = fnum + minimize1D(f, x1, Snorm, N, gradeps, astar) if (num == fnum) then write(*,*) "Exiting BFGS with error in minimize1D step" lerr = .true. @@ -85,7 +92,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(*,*) "Converged on tiny Py" + write(*,*) "Converged on tiny Py after ",i," iterations" return end if ! set up update @@ -104,6 +111,12 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) H(:,:) = H(:,:) + ((1._DP - yHy / Py) * PP(:,:) - PyH(:,:) - HyP(:,:)) / Py end do write(*,*) "Did not converge!" + write(*,*) "Py: ",Py + do k = 1, N + write(*,'("grad(",I0.2,") = ",ES12.4)') k, grad1(k) + end do + write(*,*) "|grad|: Min Mean Max" + write(*,*) minval(abs(grad1(:))), sum(abs(grad1(:))) / N, maxval(abs(grad1(:))) return contains