From cf54b5cb8e9ee6cb888054cfec2929c7a29881a7 Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 20 May 2021 05:14:32 -0400 Subject: [PATCH] Removed some obsolete code and made nfrag an argument to symba_frag_pos in preparation for testing allowing this to be changed --- src/modules/module_interfaces.f90 | 8 ++-- src/symba/symba_casedisruption.f90 | 2 +- src/symba/symba_casehitandrun.f90 | 2 +- src/symba/symba_casesupercatastrophic.f90 | 4 +- src/symba/symba_frag_pos.f90 | 50 +++++++++++------------ 5 files changed, 33 insertions(+), 33 deletions(-) diff --git a/src/modules/module_interfaces.f90 b/src/modules/module_interfaces.f90 index a4625d1b3..ebfd60a61 100644 --- a/src/modules/module_interfaces.f90 +++ b/src/modules/module_interfaces.f90 @@ -921,7 +921,7 @@ END SUBROUTINE symba_energy INTERFACE 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) + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lmerge) use swiftest_globals USE swiftest_data_structures USE module_symba @@ -929,11 +929,13 @@ 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(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) :: mass, radius + integer(I4B), intent(inout) :: nfrag + real(DP), dimension(:), intent(inout) :: m_frag, rad_frag real(DP), dimension(:,:), intent(inout) :: Ip_frag real(DP), dimension(:,:), intent(out) :: xb_frag, vb_frag, rot_frag + real(DP), intent(inout) :: Qloss logical, intent(out) :: lmerge end subroutine symba_frag_pos END INTERFACE diff --git a/src/symba/symba_casedisruption.f90 b/src/symba/symba_casedisruption.f90 index eb679daa6..e666aad99 100644 --- a/src/symba/symba_casedisruption.f90 +++ b/src/symba/symba_casedisruption.f90 @@ -72,7 +72,7 @@ function symba_casedisruption (symba_plA, family, nmergeadd, mergeadd_list, x, v ! Put the fragments on the circle surrounding the center of mass of the system call 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) + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lmerge) if (lmerge) then write(*,*) 'Should have been a merge instead.' diff --git a/src/symba/symba_casehitandrun.f90 b/src/symba/symba_casehitandrun.f90 index d32c8d730..dda8ae3bd 100644 --- a/src/symba/symba_casehitandrun.f90 +++ b/src/symba/symba_casehitandrun.f90 @@ -80,7 +80,7 @@ function symba_casehitandrun(symba_plA, family, nmergeadd, mergeadd_list, id, x, ! Put the fragments on the circle surrounding the center of mass of the system call 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, lpure, Qloss) + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lpure) if (lpure) then write(*,*) 'Should have been a pure hit and run instead' nfrag = 0 diff --git a/src/symba/symba_casesupercatastrophic.f90 b/src/symba/symba_casesupercatastrophic.f90 index 6b1b77180..04d242ff0 100644 --- a/src/symba/symba_casesupercatastrophic.f90 +++ b/src/symba/symba_casesupercatastrophic.f90 @@ -68,11 +68,11 @@ function symba_casesupercatastrophic (symba_plA, family, nmergeadd, mergeadd_lis ! Put the fragments on the circle surrounding the center of mass of the system call 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) + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lmerge) if (lmerge) then write(*,*) 'Should have been a merge instead.' - status = symba_casemerge (symba_plA, family, nmergeadd, mergeadd_list, x, v, mass, radius, L_spin, Ip, param) + status = symba_casemerge(symba_plA, family, nmergeadd, mergeadd_list, x, v, mass, radius, L_spin, Ip, param) else write(*,'("Generating ",I2.0," fragments")') nfrag status = SUPERCATASTROPHIC diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index a1d8d900b..b2aa0e551 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -10,7 +10,7 @@ module symba_frag_pos_lambda_implementation 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), dimension(NDIM) :: vcom real(DP) :: ke_target contains generic :: init => ke_objective_func_init @@ -23,7 +23,7 @@ module symba_frag_pos_lambda_implementation 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) + function abstract_objective_func(v_r_mag, v_r_unit, v_t_mag, v_t_unit, x_frag, m_frag, vcom, 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 @@ -32,14 +32,14 @@ function abstract_objective_func(v_r_mag, v_r_unit, v_t_mag, v_t_unit, x_frag, m 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), dimension(:), intent(in) :: vcom !! Barycentric velocity of collisional system center of mass 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) + type(ke_constraint) function ke_objective_func_init(lambda, v_r_unit, v_t_mag, v_t_unit, x_frag, m_frag, vcom, ke_target) implicit none ! Arguments procedure(abstract_objective_func) :: lambda !! The lambda function @@ -48,7 +48,7 @@ type(ke_constraint) function ke_objective_func_init(lambda, v_r_unit, v_t_mag, v 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), dimension(:), intent(in) :: vcom !! Barycentric velocity of collisional system center of mass real(DP), intent(in) :: ke_target !! Target kinetic energ ! Internals associate(self => ke_objective_func_init) @@ -58,7 +58,7 @@ type(ke_constraint) function ke_objective_func_init(lambda, v_r_unit, v_t_mag, v 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%vcom(:) = vcom(:) self%ke_target = ke_target end associate return @@ -84,7 +84,7 @@ function ke_objective_func_eval(self, x) result(fval) 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) + 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%vcom, self%ke_target) else error stop "KE Objective function was not initialized." end if @@ -92,8 +92,8 @@ 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) +subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radius, & + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lmerge) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Places the collision fragments on a circle oriented with a plane defined @@ -109,18 +109,20 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad type(user_input_parameters), intent(in) :: param type(symba_pl), intent(inout) :: symba_plA integer(I4B), dimension(:), intent(in) :: family - 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) :: mass, radius + integer(I4B), intent(inout) :: nfrag + real(DP), dimension(:), intent(inout) :: 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? + real(DP), intent(inout) :: Qloss ! Internals real(DP), parameter :: f_spin = 0.20_DP !! Fraction of pre-impact orbital angular momentum that is converted to fragment spin real(DP) :: mscale = 1.0_DP, rscale = 1.0_DP, vscale = 1.0_DP, tscale = 1.0_DP, Lscale = 1.0_DP, Escale = 1.0_DP! Scale factors that reduce quantities to O(~1) in the collisional system real(DP) :: mtot real(DP), dimension(NDIM) :: xcom, vcom - integer(I4B) :: i, j, nfrag, fam_size + integer(I4B) :: ii, fam_size logical, dimension(:), allocatable :: lexclude real(DP), dimension(NDIM, 2) :: rot, L_orb real(DP), dimension(:,:), allocatable :: x_frag, v_frag, v_r_unit, v_t_unit @@ -139,12 +141,11 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad ! Temporary until we get random number stuff settled for the whole project call random_seed(size=nseed) allocate(seed(nseed)) - seed = [(12*i, i=1, nseed)] + seed = [(12*ii, ii=1, nseed)] call random_seed(put=seed) allocate(x_frag, source=xb_frag) allocate(v_frag, source=vb_frag) - nfrag = size(m_frag) mtot = sum(m_frag) associate(npl => symba_plA%helio%swiftest%nbody, status => symba_plA%helio%swiftest%status) @@ -311,11 +312,11 @@ subroutine define_coordinate_system() L_orb(:, :) = 0.0_DP ! Compute orbital angular momentum of pre-impact system - do j = 1, 2 - xc(:) = x(:, j) - xcom(:) - vc(:) = v(:, j) - vcom(:) + do i = 1, 2 + xc(:) = x(:, i) - xcom(:) + vc(:) = v(:, i) - vcom(:) call util_crossproduct(xc(:), vc(:), x_cross_v(:)) - L_orb(:, j) = mass(j) * x_cross_v(:) + L_orb(:, i) = mass(i) * x_cross_v(:) end do ! Compute orbital angular momentum of pre-impact system. This will be the normal vector to the collision fragment plane @@ -343,6 +344,7 @@ 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) fam_size = size(family) @@ -460,11 +462,8 @@ subroutine shift_vector_to_origin(m_frag, vec_frag) ! Internals real(DP), dimension(NDIM) :: mvec_frag, COM_offset - integer(I4B) :: i, nfrag - real(DP) :: mtot + integer(I4B) :: i - nfrag = size(m_frag) - mtot = sum(m_frag) mvec_frag(:) = 0.0_DP do i = 1, nfrag @@ -656,7 +655,7 @@ subroutine set_fragment_radial_velocities(lmerge) return end subroutine set_fragment_radial_velocities - 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) + function ke_objective_function(v_r_mag, v_r_unit, v_t_mag, v_t_unit, x_frag, m_frag, vcom, 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 @@ -665,17 +664,16 @@ function ke_objective_function(v_r_mag, v_r_unit, v_t_mag, v_t_unit, x_frag, m_f 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), dimension(:), intent(in) :: vcom !! Barycentric velocity of collisional system center of mass real(DP), intent(in) :: ke_target !! Target kinetic energ ! 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, nfrag, nsol + integer(I4B) :: i real(DP), dimension(NDIM) :: L real(DP), dimension(:,:), allocatable :: v_shift - nfrag = size(m_frag) allocate(v_shift, mold=v_r_unit) ! Make sure the velocity magnitude stays positive do i = 1, nfrag