diff --git a/src/main/swiftest_symba.f90 b/src/main/swiftest_symba.f90 index 35ea166b2..f173f2988 100644 --- a/src/main/swiftest_symba.f90 +++ b/src/main/swiftest_symba.f90 @@ -152,35 +152,38 @@ program swiftest_symba do while ((t < tstop) .and. ((ntp0 == 0) .or. (ntp > 0))) call symba_step_eucl(t, dt, param,npl,ntp,symba_plA, symba_tpA, nplplenc, npltpenc,& plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list, mergesub_list) - - ldiscard_pl = .false. - ldiscard_tp = .false. - lfrag_add = .false. + call symba_discard_pl(t, npl, ntp, symba_plA, symba_tpA, rmin, rmax, rmaxu, qmin, qmin_coord, qmin_alo, qmin_ahi, ldiscard_pl) - call symba_discard_tp(t, npl, ntp, symba_plA, symba_tpA, dt, rmin, rmax, rmaxu, qmin, qmin_coord, & - qmin_alo, qmin_ahi, param%lrhill_present, ldiscard_tp) - call symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, mergeadd_list, nmergeadd, param) - ldiscard_pl = ldiscard_pl .or. lfrag_add - - ! These next two blocks should be streamlined/improved but right now we treat discards separately from collisions so it has to be this way - if (ldiscard_pl .or. ldiscard_tp) then + if (ldiscard_pl) then if (param%lenergy) call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Eorbit_before, Ltot) - call symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, nmergeadd, mergeadd_list, discard_plA, & - discard_tpA, ldiscard_pl, ldiscard_tp, mtiny, param) - call io_discard_write_symba(t, mtiny, npl, nsppl, nsptp, nmergesub, symba_plA, & - discard_plA%helio%swiftest, discard_tpA%helio%swiftest, mergeadd_list, mergesub_list, discard_file, param%lbig_discard) - nmergeadd = 0 - nmergesub = 0 - nsppl = 0 - nsptp = 0 - nplm = count(symba_plA%helio%swiftest%mass(1:npl) > mtiny) - - if (param%lenergy) then + call symba_rearray_pl(t, npl, symba_plA, nmergeadd, mergeadd_list, discard_plA, param) + if (param%lenergy) then call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_after, ke_spin_after, pe_after, Eorbit_after, Ltot) Ecollision = Eorbit_before - Eorbit_after ! Energy change resulting in this collisional event Total running energy offset from collision in this step symba_plA%helio%swiftest%Ecollisions = symba_plA%helio%swiftest%Ecollisions + Ecollision end if - !if (ntp > 0) call util_dist_index_pltp(nplm, ntp, symba_plA, symba_tpA) + end if + + lfrag_add = any(plplenc_list%status(1:nplplenc) == COLLISION) + if (lfrag_add) then + if (param%lenergy) call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Eorbit_before, Ltot) + call symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, nmergeadd, mergeadd_list, discard_plA, param) + if (param%lenergy) then + call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_after, ke_spin_after, pe_after, Eorbit_after, Ltot) + Ecollision = Eorbit_before - Eorbit_after ! Energy change resulting in this collisional event Total running energy offset from collision in this step + symba_plA%helio%swiftest%Ecollisions = symba_plA%helio%swiftest%Ecollisions + Ecollision + !write(*,*) 'Dke_orbit: ',ke_orbit_before - ke_orbit_after + !write(*,*) 'Dke_spin: ',ke_spin_before - ke_spin_after + !write(*,*) 'Dpe: ',pe_before - pe_after + end if + end if + + call symba_discard_tp(t, npl, ntp, symba_plA, symba_tpA, dt, rmin, rmax, rmaxu, qmin, qmin_coord, & + qmin_alo, qmin_ahi, param%lrhill_present, ldiscard_tp) + + if (ldiscard_pl .or. ldiscard_tp .or. lfrag_add) then + call io_discard_write_symba(t, mtiny, npl, nsppl, nsptp, nmergesub, symba_plA, & + discard_plA%helio%swiftest, discard_tpA%helio%swiftest, mergeadd_list, mergesub_list, discard_file, param%lbig_discard) end if iloop = iloop + 1 diff --git a/src/modules/module_interfaces.f90 b/src/modules/module_interfaces.f90 index 7f0ea3e52..14780cb80 100644 --- a/src/modules/module_interfaces.f90 +++ b/src/modules/module_interfaces.f90 @@ -796,20 +796,21 @@ END SUBROUTINE symba_chk END INTERFACE INTERFACE - SUBROUTINE symba_collision(t, symba_plA, nplplenc, plplenc_list, ldiscard, mergeadd_list, nmergeadd, param) - USE swiftest_globals - USE swiftest_data_structures - USE module_helio - USE module_symba - IMPLICIT NONE - real(DP), intent(in) :: t - integer(I4B), intent(inout) :: nplplenc, nmergeadd - type(symba_pl) :: symba_plA - type(symba_plplenc), intent(inout) :: plplenc_list - type(symba_merger), intent(inout) :: mergeadd_list - logical, intent(inout) :: ldiscard - type(user_input_parameters),intent(inout) :: param - END SUBROUTINE symba_collision + subroutine symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, nmergeadd, mergeadd_list, discard_plA, param) + use swiftest_globals + use swiftest_data_structures + use module_helio + use module_symba + implicit none + real(DP), intent(in) :: t + type(symba_pl) :: symba_plA + integer(I4B), intent(inout) :: nplplenc, nmergeadd + type(symba_plplenc), intent(inout) :: plplenc_list + logical, intent(inout) :: lfrag_add + type(symba_merger), intent(inout) :: mergeadd_list + type(symba_pl), intent(inout) :: discard_plA + type(user_input_parameters),intent(inout) :: param + end subroutine symba_collision END INTERFACE INTERFACE @@ -1074,24 +1075,19 @@ END SUBROUTINE symba_peri END INTERFACE INTERFACE - subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, nmergeadd, mergeadd_list, discard_plA, & - discard_tpA, ldiscard_pl, ldiscard_tp, mtiny, param) + subroutine symba_rearray_pl(t, npl, symba_plA, nmergeadd, mergeadd_list, discard_plA, param) use swiftest_globals use swiftest_data_structures use module_symba implicit none real(DP), intent(in) :: t - integer(I4B), intent(inout) :: npl, nplm, ntp, nsppl, nsptp + integer(I4B), intent(inout) :: npl integer(I4B), intent(in) :: nmergeadd type(symba_pl), intent(inout) :: symba_plA - type(symba_tp), intent(inout) :: symba_tpA - type(symba_tp), intent(inout) :: discard_tpA type(symba_pl), intent(inout) :: discard_plA type(symba_merger), intent(inout) :: mergeadd_list - logical(LGT), intent(in) :: ldiscard_pl, ldiscard_tp - real(DP), intent(in) :: mtiny type(user_input_parameters), intent(in) :: param - end subroutine symba_rearray + end subroutine symba_rearray_pl END INTERFACE diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index e4912975c..24bde9954 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -1,4 +1,4 @@ -subroutine symba_collision (t, symba_plA, nplplenc, plplenc_list, ldiscard, mergeadd_list, nmergeadd, param) +subroutine symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, nmergeadd, mergeadd_list, discard_plA, param) !! author: Jennifer L.L. Pouplin, Carlisle A. wishard, and David A. Minton !! !! Check for merger between planets in SyMBA. If the user has turned on the FRAGMENTATION feature, it will call the @@ -12,15 +12,16 @@ subroutine symba_collision (t, symba_plA, nplplenc, plplenc_list, ldiscard, merg use module_symba use module_interfaces, EXCEPT_THIS_ONE => symba_collision implicit none - + ! Arguments real(DP), intent(in) :: t - integer(I4B), intent(inout) :: nplplenc, nmergeadd type(symba_pl) :: symba_plA + integer(I4B), intent(inout) :: nplplenc, nmergeadd type(symba_plplenc), intent(inout) :: plplenc_list + logical, intent(inout) :: lfrag_add type(symba_merger), intent(inout) :: mergeadd_list - logical, intent(inout) :: ldiscard + type(symba_pl), intent(inout) :: discard_plA type(user_input_parameters),intent(inout) :: param - + ! internals integer(I4B), parameter :: NRES = 3 !! Number of collisional product results integer(I4B) :: i, j, k, index_enc, index_coll, jtarg, jproj real(DP), dimension(NRES) :: mass_res @@ -45,19 +46,11 @@ subroutine symba_collision (t, symba_plA, nplplenc, plplenc_list, ldiscard, merg logical, dimension(nplplenc) :: lplpl_collision logical, dimension(:), allocatable :: lplpl_unique_parent integer(I4B), dimension(:), pointer :: plparent - - ! TESTING - logical, save :: lfirst = .true. - real(DP), save :: Minitial - real(DP) :: Msystem, Madd, Mdiscard + logical :: ldiscard ! First determine the collisional regime for each colliding pair associate(npl => symba_plA%helio%swiftest%nbody, xbpl => symba_plA%helio%swiftest%xb, statpl => symba_plA%helio%swiftest%status, idpl => symba_plA%helio%swiftest%id, & idx1 => plplenc_list%index1, idx2 => plplenc_list%index2, plparent => symba_plA%kin%parent) - if (lfirst) then - Minitial = sum(symba_plA%helio%swiftest%mass(1:npl)) - lfirst = .false. - end if lplpl_collision(:) = plplenc_list%status(1:nplplenc) == COLLISION ldiscard = any(lplpl_collision) if (.not.ldiscard) return @@ -311,6 +304,8 @@ subroutine symba_collision (t, symba_plA, nplplenc, plplenc_list, ldiscard, merg end do deallocate(family) + call symba_rearray_pl(t, npl, symba_plA, nmergeadd, mergeadd_list, discard_plA, param) + end do end associate diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index 1e98ea6a6..2e418aa7f 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -110,7 +110,15 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad !write(*, "(' ---------------------------------------------------------------------------')") !write(*,fmtlabel) ' Q_loss |',-Qloss / abs(Etot_before) !write(*, "(' ---------------------------------------------------------------------------')") - +! + ! Put any residual angular momentum into the spin before constraining energy + !L_residual(:) = Ltot_before(:) - Ltot_after(:) + !L_spin_frag(:) = L_residual(:) / nfrag + !do i = 1, nfrag + ! ke_spin_after = ke_spin_after - m_frag(i) * Ip_frag(3,i) * rad_frag(i)**2 * norm2(rot_frag(:,i)) + ! rot_frag(:,i) = rot_frag(:,i) + L_spin_frag(:) / (Ip_frag(3, i) * m_frag(i) * rad_frag(i)**2) + ! ke_spin_after = ke_spin_after + m_frag(i) * Ip_frag(3,i) * rad_frag(i)**2 * norm2(rot_frag(:,i)) + !end do ! 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 @@ -132,24 +140,17 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad Lmag_after = norm2(Ltot_after(:)) lmerge = lmerge .or. ((Etot_after - Etot_before) / abs(Etot_before) > 0._DP) - if (.not.lmerge) then - L_residual(:) = Ltot_before(:) - Ltot_after(:) - L_spin_frag(:) = L_residual(:) / nfrag - do i = 1, nfrag - rot_frag(:,i) = rot_frag(:,i) + L_spin_frag(:) / (Ip_frag(3, i) * m_frag(i) * rad_frag(i)**2) - end do - end if -! write(*, "(' ---------------------------------------------------------------------------')") -! write(*, "(' | T_orb T_spin T pe Etot Ltot')") -! 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), & -! (Lmag_after - Lma_before) / Lmag_before -! write(*, "(' ---------------------------------------------------------------------------')") + !write(*, "(' ---------------------------------------------------------------------------')") + !write(*, "(' | T_orb T_spin T pe Etot Ltot')") + !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), & + ! (Lmag_after - Lmag_before) / Lmag_before + !write(*, "(' ---------------------------------------------------------------------------')") !**************************************************************************************************************** end associate @@ -229,8 +230,10 @@ subroutine symba_frag_pos_initialize_fragments(nfrag, xcom, vcom, x, v, L_orb, L x_frag(:, i) = r_frag_norm * (frag_vec(1) * x_col_unit(:) + frag_vec(2) * y_col_unit(:)) end do call symba_frag_pos_com_adjust(xcom, m_frag, x_frag) + v_frag(:,:) = 0.0_DP call symba_frag_pos_ang_mtm(nfrag, xcom, vcom, L_orb, L_spin, mass, radius, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag) + call symba_frag_pos_com_adjust(vcom, m_frag, v_frag) return end subroutine symba_frag_pos_initialize_fragments @@ -272,13 +275,12 @@ subroutine symba_frag_pos_ang_mtm(nfrag, xcom, vcom, L_orb, L_spin, mass, radius rmag = norm2(x_frag(:,i)) v_r_unit(:) = x_frag(:,i) / rmag call util_crossproduct(h_unit(:), v_r_unit(:), v_t_unit(:)) ! make a unit vector in the tangential velocity direction wrt the angular momentum plane - v_frag(:,i) = L_orb_frag_mag / (m_frag(i) * rmag) * v_t_unit(:) ! Distribute the angular momentum equally amongst the fragments + v_frag(:,i) = v_frag(:,i) + L_orb_frag_mag / (m_frag(i) * rmag) * v_t_unit(:) ! Distribute the angular momentum equally amongst the fragments end do return end subroutine symba_frag_pos_ang_mtm - subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, m_frag, x_frag, v_frag, ke_target, lmerge) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! @@ -334,7 +336,6 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, m_frag, x_frag, v_fr lmerge = .true. end if - ! Shift the fragments into the system barycenter frame do i = 1, nfrag v_frag(:,i) = v_corrected * mtot / m_frag(i) * v_r_unit(:, i) + v_t(:, i) @@ -444,11 +445,13 @@ subroutine symba_frag_pos_energy_calc(npl, symba_plA, lexclude, ke_orbit, ke_spi 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 + elsewhere + symba_plwksp%helio%swiftest%status(1:npl_new) = ACTIVE end where - nplm = count(symba_plwksp%helio%swiftest%mass > param%mtiny) + nplm = count((symba_plwksp%helio%swiftest%mass(:) > param%mtiny) .and. .not.lexclude(:)) 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) diff --git a/src/symba/symba_rearray.f90 b/src/symba/symba_rearray.f90 index 1789b607d..11c834d3d 100644 --- a/src/symba/symba_rearray.f90 +++ b/src/symba/symba_rearray.f90 @@ -1,5 +1,4 @@ -subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, nmergeadd, mergeadd_list, discard_plA,& - discard_tpA, ldiscard_pl, ldiscard_tp, mtiny, param) +subroutine symba_rearray_pl(t, npl, symba_plA, nmergeadd, mergeadd_list, discard_plA, param) !! Author: the Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! !! Clean up tp and pl arrays to remove discarded bodies and add new bodies @@ -8,29 +7,25 @@ subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, use module_swiftestalloc use module_helio use module_symba - use module_interfaces, EXCEPT_THIS_ONE => symba_rearray + use module_interfaces, EXCEPT_THIS_ONE => symba_rearray_pl implicit none ! Arguments - real(DP), intent(in) :: t - integer(I4B), intent(inout) :: npl, nplm, ntp, nsppl, nsptp - integer(I4B), intent(in) :: nmergeadd + real(DP), intent(in) :: t + integer(I4B), intent(inout) :: npl + integer(I4B), intent(in) :: nmergeadd type(symba_pl), intent(inout) :: symba_plA - type(symba_tp), intent(inout) :: symba_tpA - type(symba_tp), intent(inout) :: discard_tpA type(symba_pl), intent(inout) :: discard_plA type(symba_merger), intent(inout) :: mergeadd_list - logical, intent(in) :: ldiscard_pl, ldiscard_tp - real(DP), intent(in) :: mtiny type(user_input_parameters), intent(in) :: param ! Internals - integer(I4B) :: i, j, nkpl, nktp, ntot, dlo + integer(I4B) :: i, j, nsppl, nkpl, ntot, dlo, nplm logical, dimension(:), allocatable :: discard_l_pl, add_l_pl - logical, dimension(ntp) :: discard_l_tp - logical :: lescape + logical :: lescape, ldiscard_pl - if (any(symba_plA%helio%swiftest%status(1:npl) /= ACTIVE)) then + ldiscard_pl = any(symba_plA%helio%swiftest%status(1:npl) /= ACTIVE) + if (ldiscard_pl) then ! Deal with the central body/system discards if there are any do i = 1, npl @@ -151,57 +146,10 @@ subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, call util_hills(npl, symba_plA%helio%swiftest) - nplm = count(symba_plA%helio%swiftest%mass > mtiny) + nplm = count(symba_plA%helio%swiftest%mass > param%mtiny) CALL util_dist_index_plpl(npl, nplm, symba_plA) - end if - if (ldiscard_tp) then - nktp = 0 - nsptp = 0 - - discard_l_tp(1:ntp) = (symba_tpA%helio%swiftest%status(1:ntp) /= ACTIVE) - nsptp = count(discard_l_tp) - nktp = ntp - nsptp - - dlo = 1 - if (allocated(discard_tpA%helio%swiftest%id)) then ! We alredy made a discard list in this step, so we need to append to it - nsptp = nsptp + discard_tpA%helio%swiftest%nbody - dlo = dlo + discard_tpA%helio%swiftest%nbody - !call util_resize_tp(discard_plA, nsppl, discard_tplA%helio%swiftest%nbody) !TODO: Implement this. Probably best to wait until this gets OOFed - else - call symba_tp_allocate(discard_tpA, nsptp) - end if - - discard_tpA%helio%swiftest%id(1:nsptp) = pack(symba_tpA%helio%swiftest%id(1:ntp), discard_l_tp) - discard_tpA%helio%swiftest%status(1:nsptp) = pack(symba_tpA%helio%swiftest%status(1:ntp), discard_l_tp) - discard_tpA%helio%swiftest%isperi(1:nsptp) = pack(symba_tpA%helio%swiftest%isperi(1:ntp), discard_l_tp) - discard_tpA%helio%swiftest%peri(1:nsptp) = pack(symba_tpA%helio%swiftest%peri(1:ntp), discard_l_tp) - discard_tpA%helio%swiftest%atp(1:nsptp) = pack(symba_tpA%helio%swiftest%atp(1:ntp), discard_l_tp) - do i = 1, NDIM - discard_tpA%helio%swiftest%xh(i,1:nsptp) = pack(symba_tpA%helio%swiftest%xh(i,1:ntp), discard_l_tp) - discard_tpA%helio%swiftest%vh(i,1:nsptp) = pack(symba_tpA%helio%swiftest%vh(i,1:ntp), discard_l_tp) - discard_tpA%helio%swiftest%xb(i,1:nsptp) = pack(symba_tpA%helio%swiftest%xb(i,1:ntp), discard_l_tp) - discard_tpA%helio%swiftest%vb(i,1:nsptp) = pack(symba_tpA%helio%swiftest%vb(i,1:ntp), discard_l_tp) - end do - - symba_tpA%helio%swiftest%id(1:nktp) = pack(symba_tpA%helio%swiftest%id(1:ntp), .not. discard_l_tp) - symba_tpA%helio%swiftest%status(1:nktp) = pack(symba_tpA%helio%swiftest%status(1:ntp), .not. discard_l_tp) - symba_tpA%helio%swiftest%isperi(1:nktp) = pack(symba_tpA%helio%swiftest%isperi(1:ntp), .not. discard_l_tp) - symba_tpA%helio%swiftest%peri(1:nktp) = pack(symba_tpA%helio%swiftest%peri(1:ntp), .not. discard_l_tp) - symba_tpA%helio%swiftest%atp(1:nktp) = pack(symba_tpA%helio%swiftest%atp(1:ntp), .not. discard_l_tp) - do i = 1, NDIM - symba_tpA%helio%swiftest%xh(i,1:nktp) = pack(symba_tpA%helio%swiftest%xh(i,1:ntp), .not. discard_l_tp) - symba_tpA%helio%swiftest%vh(i,1:nktp) = pack(symba_tpA%helio%swiftest%vh(i,1:ntp), .not. discard_l_tp) - symba_tpA%helio%swiftest%xb(i,1:nktp) = pack(symba_tpA%helio%swiftest%xb(i,1:ntp), .not. discard_l_tp) - symba_tpA%helio%swiftest%vb(i,1:nktp) = pack(symba_tpA%helio%swiftest%vb(i,1:ntp), .not. discard_l_tp) - symba_tpA%helio%ah(i,1:nktp) = pack(symba_tpA%helio%ah(i,1:ntp), .not. discard_l_tp) - end do - ntp = nktp - symba_tpA%helio%swiftest%nbody = ntp - nplm = count(symba_plA%helio%swiftest%mass>mtiny) - call coord_b2h_tp(ntp, symba_tpA%helio%swiftest, symba_plA%helio%swiftest) - end if -end subroutine symba_rearray +end subroutine symba_rearray_pl