From f903e34c5cf3e24282093a188a794e167f265792 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 17 Dec 2022 07:41:40 -0500 Subject: [PATCH] More cleaning --- src/fraggle/fraggle_generate.f90 | 165 ++++++++++++++++-------------- src/fraggle/fraggle_io.f90 | 52 +++++----- src/fraggle/fraggle_regime.f90 | 14 +-- src/fraggle/fraggle_set.f90 | 167 +++++++++++++++---------------- src/fraggle/fraggle_util.f90 | 16 +-- src/modules/fraggle_classes.f90 | 71 +++++++------ src/symba/symba_collision.f90 | 16 +-- 7 files changed, 254 insertions(+), 247 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index ffc7acf7a..67b375940 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -47,7 +47,7 @@ module subroutine fraggle_generate_fragments(self, system, param, lfailure) select type(fragments => self%fragments) class is (fraggle_fragments) - associate(collision => self, impactors => self%impactors, nfrag => fragments%nbody, pl => system%pl) + associate(collision_system => self, impactors => self%impactors, nfrag => fragments%nbody, pl => system%pl) write(message,*) nfrag call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle generating " // trim(adjustl(message)) // " fragments.") @@ -69,12 +69,12 @@ module subroutine fraggle_generate_fragments(self, system, param, lfailure) lk_plpl = .false. end if - call fragments%set_natural_scale(impactors) + call collision_system%set_natural_scale() call fragments%reset() ! Calculate the initial energy of the system without the collisional family - call collision%get_energy_and_momentum(system, param, lbefore=.true.) + call collision_system%get_energy_and_momentum(system, param, lbefore=.true.) ! Start out the fragments close to the initial separation distance. This will be increased if there is any overlap or we fail to find a solution r_max_start = 1.2_DP * .mag.(impactors%rb(:,2) - impactors%rb(:,1)) @@ -92,38 +92,38 @@ module subroutine fraggle_generate_fragments(self, system, param, lfailure) lfailure = .false. call ieee_set_flag(ieee_all, .false.) ! Set all fpe flags to quiet - call fraggle_generate_pos_vec(fragments, impactors, r_max_start) - call fragments%set_coordinate_system(impactors) + call fraggle_generate_pos_vec(collision_system, r_max_start) + call collision_system%set_coordinate_system() ! Initial velocity guess will be the barycentric velocity of the colliding system so that the budgets are based on the much smaller collisional-frame velocities do concurrent (i = 1:nfrag) - fragments%vb(:, i) = fragments%vbcom(:) + fragments%vb(:, i) = impactors%vbcom(:) end do - call fragments%get_energy_and_momentum(impactors, system, param, lbefore=.false.) - call fragments%set_budgets() + call collision_system%get_energy_and_momentum(system, param, lbefore=.false.) + call collision_system%set_budgets() - call fraggle_generate_spins(fragments, impactors, f_spin, lfailure) + call fraggle_generate_spins(collision_system, f_spin, lfailure) if (lfailure) then call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed to find spins") cycle end if - call fraggle_generate_tan_vel(fragments, impactors, lfailure) + call fraggle_generate_tan_vel(collision_system, lfailure) if (lfailure) then call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed to find tangential velocities") cycle end if - call fraggle_generate_rad_vel(fragments, impactors, lfailure) + call fraggle_generate_rad_vel(collision_system, lfailure) if (lfailure) then call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed to find radial velocities") cycle end if - call fragments%get_energy_and_momentum(impactors, system, param, lbefore=.false.) - dEtot = fragments%Etot_after - fragments%Etot_before - dLmag = .mag. (fragments%Ltot_after(:) - fragments%Ltot_before(:)) + call collision_system%get_energy_and_momentum(system, param, lbefore=.false.) + dEtot = collision_system%Etot(2) - collision_system%Etot(1) + dLmag = .mag. (collision_system%Ltot(:,2) - collision_system%Ltot(:,1)) exit lfailure = ((abs(dEtot + impactors%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) @@ -134,9 +134,9 @@ module subroutine fraggle_generate_fragments(self, system, param, lfailure) cycle end if - lfailure = ((abs(dLmag) / (.mag.fragments%Ltot_before)) > FRAGGLE_LTOL) + lfailure = ((abs(dLmag) / (.mag.collision_system%Ltot(:,1))) > FRAGGLE_LTOL) if (lfailure) then - write(message,*) dLmag / (.mag.fragments%Ltot_before(:)) + write(message,*) dLmag / (.mag.collision_system%Ltot(:,1)) call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed due to high angular momentum error: " // & trim(adjustl(message))) cycle @@ -159,7 +159,7 @@ module subroutine fraggle_generate_fragments(self, system, param, lfailure) trim(adjustl(message)) // " tries") end if - call fragments%set_original_scale(impactors) + call collision_system%set_original_scale() ! Restore the big array if (lk_plpl) call pl%flatten(param) @@ -171,7 +171,7 @@ module subroutine fraggle_generate_fragments(self, system, param, lfailure) end subroutine fraggle_generate_fragments - subroutine fraggle_generate_pos_vec(fragments, impactors, r_max_start) + subroutine fraggle_generate_pos_vec(collision_system, r_max_start) !! 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 @@ -179,9 +179,8 @@ subroutine fraggle_generate_pos_vec(fragments, impactors, r_max_start) !! The initial positions do not conserve energy or momentum, so these need to be adjusted later. implicit none ! Arguments - class(fraggle_fragments), intent(inout) :: fragments !! Fraggle fragment system object - class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object - real(DP), intent(in) :: r_max_start !! Initial guess for the starting maximum radial distance of fragments + class(fraggle_system), intent(inout) :: collision_system !! Fraggle collision system object + real(DP), intent(in) :: r_max_start !! Initial guess for the starting maximum radial distance of fragments ! Internals real(DP) :: dis, rad, r_max, fdistort real(DP), dimension(NDIM) :: runit, vunit @@ -189,7 +188,7 @@ subroutine fraggle_generate_pos_vec(fragments, impactors, r_max_start) integer(I4B) :: i, j logical :: lnoncat, lhitandrun - associate(nfrag => fragments%nbody) + associate(fragments => collision_system%fragments, impactors => collision_system%impactors, nfrag => collision_system%fragments%nbody) allocate(loverlap(nfrag)) lnoncat = (impactors%regime /= COLLRESOLVE_REGIME_SUPERCATASTROPHIC) ! For non-catastrophic impacts, make the fragments act like ejecta and point away from the impact point @@ -216,8 +215,8 @@ subroutine fraggle_generate_pos_vec(fragments, impactors, r_max_start) call random_number(fragments%rc(:,3:nfrag)) loverlap(:) = .true. do while (any(loverlap(3:nfrag))) - fragments%rc(:, 1) = impactors%rb(:, 1) - fragments%rbcom(:) - fragments%rc(:, 2) = impactors%rb(:, 2) - fragments%rbcom(:) + fragments%rc(:, 1) = impactors%rb(:, 1) - impactors%rbcom(:) + fragments%rc(:, 2) = impactors%rb(:, 2) - impactors%rbcom(:) r_max = r_max + 0.1_DP * rad do i = 3, nfrag if (loverlap(i)) then @@ -225,7 +224,7 @@ subroutine fraggle_generate_pos_vec(fragments, impactors, r_max_start) fragments%rc(:,i) = 2 * (fragments%rc(:, i) - 0.5_DP) fragments%rc(:, i) = fragments%rc(:,i) + fdistort * vunit(:) fragments%rc(:, i) = r_max * fragments%rc(:, i) - fragments%rc(:, i) = fragments%rc(:, i) + (fragments%rbimp(:) - fragments%rbcom(:)) ! Shift the center of the fragment cloud to the impact point rather than the CoM + fragments%rc(:, i) = fragments%rc(:, i) + (impactors%rbimp(:) - impactors%rbcom(:)) ! Shift the center of the fragment cloud to the impact point rather than the CoM !if (lnoncat .and. dot_product(fragments%rc(:,i), runit(:)) < 0.0_DP) fragments%rc(:, i) = -fragments%rc(:, i) ! Make sure the fragment cloud points away from the impact point end if end do @@ -238,24 +237,24 @@ subroutine fraggle_generate_pos_vec(fragments, impactors, r_max_start) end do end do call fraggle_util_shift_vector_to_origin(fragments%mass, fragments%rc) - call fragments%set_coordinate_system(impactors) + call collision_system%set_coordinate_system() do concurrent(i = 1:nfrag) - fragments%rb(:,i) = fragments%rc(:,i) + fragments%rbcom(:) + fragments%rb(:,i) = fragments%rc(:,i) + impactors%rbcom(:) end do - fragments%rbcom(:) = 0.0_DP + impactors%rbcom(:) = 0.0_DP do concurrent(i = 1:nfrag) - fragments%rbcom(:) = fragments%rbcom(:) + fragments%mass(i) * fragments%rb(:,i) + impactors%rbcom(:) = impactors%rbcom(:) + fragments%mass(i) * fragments%rb(:,i) end do - fragments%rbcom(:) = fragments%rbcom(:) / fragments%mtot + impactors%rbcom(:) = impactors%rbcom(:) / fragments%mtot end associate return end subroutine fraggle_generate_pos_vec - subroutine fraggle_generate_spins(fragments, impactors, f_spin, lfailure) + subroutine fraggle_generate_spins(collision_system, f_spin, lfailure) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Calculates the spins of a collection of fragments such that they conserve angular momentum without blowing the fragment kinetic energy budget. @@ -263,19 +262,20 @@ subroutine fraggle_generate_spins(fragments, impactors, f_spin, lfailure) !! A failure will trigger a restructuring of the fragments so we will try new values of the radial position distribution. implicit none ! Arguments - class(fraggle_fragments), intent(inout) :: fragments !! Fraggle fragment system object - class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object - real(DP), intent(in) :: f_spin !! Fraction of energy or momentum that goes into spin (whichever gives the lowest kinetic energy) - logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds! + class(fraggle_system), intent(inout) :: collision_system !! Fraggle collision system object + real(DP), intent(in) :: f_spin !! Fraction of energy or momentum that goes into spin (whichever gives the lowest kinetic energy) + logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds! ! Internals real(DP), dimension(NDIM) :: L_remainder, rot_L, rot_ke, L - real(DP), dimension(NDIM,fragments%nbody) :: frot_rand ! The random rotation factor applied to fragments + real(DP), dimension(NDIM,collision_system%fragments%nbody) :: frot_rand ! The random rotation factor applied to fragments real(DP), parameter :: frot_rand_mag = 1.50_DP ! The magnitude of the rotation variation to apply to the fragments integer(I4B) :: i character(len=STRMAX) :: message real(DP) :: ke_remainder, ke - associate(nfrag => fragments%nbody) + associate(impactors => collision_system%impactors, nfrag => collision_system%fragments%nbody) + select type(fragments => collision_system%fragments) + class is (fraggle_fragments) lfailure = .false. L_remainder(:) = fragments%L_budget(:) ke_remainder = fragments%ke_budget @@ -337,13 +337,14 @@ subroutine fraggle_generate_spins(fragments, impactors, f_spin, lfailure) call io_log_one_message(FRAGGLE_LOG_OUT, "ke_remainder : " // trim(adjustl(message))) end if + end select end associate return end subroutine fraggle_generate_spins - subroutine fraggle_generate_tan_vel(fragments, impactors, lfailure) + subroutine fraggle_generate_tan_vel(collision_system, lfailure) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Adjusts the tangential velocities and spins of a collection of fragments such that they conserve angular momentum without blowing the fragment kinetic energy budget. @@ -358,8 +359,7 @@ subroutine fraggle_generate_tan_vel(fragments, impactors, lfailure) !! A failure will trigger a restructuring of the fragments so we will try new values of the radial position distribution. implicit none ! Arguments - class(fraggle_fragments), intent(inout) :: fragments !! Fraggle fragment system object - class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object + class(fraggle_system), intent(inout) :: collision_system !! Fraggle collision system object logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds ! Internals integer(I4B) :: i, try @@ -370,13 +370,15 @@ subroutine fraggle_generate_tan_vel(fragments, impactors, lfailure) integer(I4B), parameter :: MAXTRY = 100 real(DP) :: tol real(DP), dimension(:), allocatable :: v_t_initial, v_t_output - real(DP), dimension(fragments%nbody) :: kefrag, vnoise + real(DP), dimension(collision_system%fragments%nbody) :: kefrag, vnoise type(lambda_obj_err) :: objective_function real(DP), dimension(NDIM) :: L_frag_tot character(len=STRMAX) :: message real(DP) :: ke_diff - associate(nfrag => fragments%nbody) + associate(impactors => collision_system%impactors, nfrag => collision_system%fragments%nbody) + select type(fragments => collision_system%fragments) + class is (fraggle_fragments) lfailure = .false. allocate(v_t_initial, mold=fragments%v_t_mag) @@ -406,9 +408,9 @@ subroutine fraggle_generate_tan_vel(fragments, impactors, lfailure) ! Perform one final shift of the radial velocity vectors to align with the center of mass of the collisional system (the origin) fragments%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%v_r_unit(:,1:nfrag), fragments%v_t_mag(1:nfrag), & - fragments%v_t_unit(:,1:nfrag), fragments%mass(1:nfrag), fragments%vbcom(:)) + fragments%v_t_unit(:,1:nfrag), fragments%mass(1:nfrag), impactors%vbcom(:)) do concurrent (i = 1:nfrag) - fragments%vc(:,i) = fragments%vb(:,i) - fragments%vbcom(:) + fragments%vc(:,i) = fragments%vb(:,i) - impactors%vbcom(:) end do ! Now do a kinetic energy budget check to make sure we are still within the budget. @@ -429,7 +431,7 @@ subroutine fraggle_generate_tan_vel(fragments, impactors, lfailure) call io_log_one_message(FRAGGLE_LOG_OUT, "Tangential velocity failure diagnostics") call fragments%get_angular_momentum() L_frag_tot = fragments%Lspin(:) + fragments%Lorbit(:) - write(message, *) .mag.(fragments%L_budget(:) - L_frag_tot(:)) / (.mag.fragments%Ltot_before(:)) + write(message, *) .mag.(fragments%L_budget(:) - L_frag_tot(:)) / (.mag.collision_system%Ltot(:,1)) call io_log_one_message(FRAGGLE_LOG_OUT, "|L_remainder| : " // trim(adjustl(message))) write(message, *) fragments%ke_budget call io_log_one_message(FRAGGLE_LOG_OUT, "ke_budget : " // trim(adjustl(message))) @@ -440,6 +442,7 @@ subroutine fraggle_generate_tan_vel(fragments, impactors, lfailure) write(message, *) fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit call io_log_one_message(FRAGGLE_LOG_OUT, "ke_radial : " // trim(adjustl(message))) end if + end select end associate return @@ -461,7 +464,9 @@ function solve_fragment_tan_vel(lfailure, v_t_mag_input) result(v_t_mag_output) 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(NDIM) :: L_lin_others, L_orb_others, L, vtmp - + + select type(fragments => collision_system%fragments) + class is (fraggle_fragments) associate(nfrag => fragments%nbody) lfailure = .false. ! We have 6 constraint equations (2 vector constraints in 3 dimensions each) @@ -486,6 +491,7 @@ function solve_fragment_tan_vel(lfailure, v_t_mag_input) result(v_t_mag_output) v_t_mag_output(1:6) = util_solve_linear_system(A, b, 6, lfailure) if (present(v_t_mag_input)) v_t_mag_output(7:nfrag) = v_t_mag_input(:) end associate + end select return end function solve_fragment_tan_vel @@ -502,15 +508,17 @@ function tangential_objective_function(v_t_mag_input, lfailure) result(fval) real(DP) :: fval ! Internals integer(I4B) :: i - real(DP), dimension(NDIM,fragments%nbody) :: v_shift - real(DP), dimension(fragments%nbody) :: v_t_new, kearr + real(DP), dimension(NDIM,collision_system%fragments%nbody) :: v_shift + real(DP), dimension(collision_system%fragments%nbody) :: v_t_new, kearr real(DP) :: keo - - associate(nfrag => fragments%nbody) + + select type(fragments => collision_system%fragments) + class is (fraggle_fragments) + associate(impactors => collision_system%impactors, nfrag => fragments%nbody) lfailure = .false. v_t_new(:) = solve_fragment_tan_vel(v_t_mag_input=v_t_mag_input(:), lfailure=lfailure) - v_shift(:,:) = fraggle_util_vmag_to_vb(fragments%v_r_mag, fragments%v_r_unit, v_t_new, fragments%v_t_unit, fragments%mass, fragments%vbcom) + v_shift(:,:) = fraggle_util_vmag_to_vb(fragments%v_r_mag, fragments%v_r_unit, v_t_new, fragments%v_t_unit, fragments%mass, impactors%vbcom) kearr = 0.0_DP do concurrent(i = 1:nfrag) @@ -520,6 +528,7 @@ function tangential_objective_function(v_t_mag_input, lfailure) result(fval) fval = keo lfailure = .false. end associate + end select return end function tangential_objective_function @@ -527,15 +536,14 @@ end function tangential_objective_function end subroutine fraggle_generate_tan_vel - subroutine fraggle_generate_rad_vel(fragments, impactors, lfailure) + subroutine fraggle_generate_rad_vel(collision_system, lfailure) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! !! Adjust the fragment velocities to set the fragment orbital kinetic energy. This will minimize the difference between the fragment kinetic energy and the energy budget implicit none ! Arguments - class(fraggle_fragments), intent(inout) :: fragments !! Fraggle fragment system object - class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object + class(fraggle_system), intent(inout) :: collision_system !! Fraggle collision system object logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds! ! Internals real(DP), parameter :: TOL_MIN = FRAGGLE_ETOL ! This needs to be more accurate than the tangential step, as we are trying to minimize the total residual energy @@ -545,11 +553,13 @@ subroutine fraggle_generate_rad_vel(fragments, impactors, lfailure) real(DP) :: ke_radial, tol integer(I4B) :: i real(DP), dimension(:), allocatable :: v_r_initial - real(DP), dimension(fragments%nbody) :: vnoise + real(DP), dimension(collision_system%fragments%nbody) :: vnoise type(lambda_obj) :: objective_function character(len=STRMAX) :: message - associate(nfrag => fragments%nbody) + associate(impactors => collision_system%impactors, nfrag => collision_system%fragments%nbody) + select type(fragments => collision_system%fragments) + class is (fraggle_fragments) ! Set the "target" ke for the radial component allocate(v_r_initial, source=fragments%v_r_mag) @@ -580,9 +590,9 @@ subroutine fraggle_generate_rad_vel(fragments, impactors, lfailure) fragments%ke_orbit = 0.0_DP fragments%ke_spin = 0.0_DP fragments%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%v_r_unit(:,1:nfrag), & - fragments%v_t_mag(1:nfrag), fragments%v_t_unit(:,1:nfrag), fragments%mass(1:nfrag), fragments%vbcom(:)) + fragments%v_t_mag(1:nfrag), fragments%v_t_unit(:,1:nfrag), fragments%mass(1:nfrag), impactors%vbcom(:)) do i = 1, nfrag - fragments%vc(:, i) = fragments%vb(:, i) - fragments%vbcom(:) + fragments%vc(:, i) = fragments%vb(:, i) - impactors%vbcom(:) fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * norm2(fragments%vb(:, i)) fragments%ke_spin = fragments%ke_spin + fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * norm2(fragments%rot(:,i)) end do @@ -603,6 +613,7 @@ subroutine fraggle_generate_rad_vel(fragments, impactors, lfailure) call io_log_one_message(FRAGGLE_LOG_OUT, "ke_remainder : " // trim(adjustl(message))) end if + end select end associate return @@ -620,23 +631,29 @@ function radial_objective_function(v_r_mag_input) result(fval) ! Internals integer(I4B) :: i real(DP), dimension(:,:), allocatable :: v_shift - real(DP), dimension(fragments%nbody) :: kearr + real(DP), dimension(collision_system%fragments%nbody) :: kearr real(DP) :: keo, ke_radial, rotmag2, vmag2 - - allocate(v_shift, mold=fragments%vb) - v_shift(:,:) = fraggle_util_vmag_to_vb(v_r_mag_input, fragments%v_r_unit, fragments%v_t_mag, fragments%v_t_unit, fragments%mass, fragments%vbcom) - !$omp do simd firstprivate(fragments) - do i = 1,fragments%nbody - rotmag2 = fragments%rot(1,i)**2 + fragments%rot(2,i)**2 + fragments%rot(3,i)**2 - vmag2 = v_shift(1,i)**2 + v_shift(2,i)**2 + v_shift(3,i)**2 - kearr(i) = fragments%mass(i) * (fragments%Ip(3, i) * fragments%radius(i)**2 * rotmag2 + vmag2) - end do - !$omp end do simd - keo = 2 * fragments%ke_budget - sum(kearr(:)) - ke_radial = fragments%ke_budget - fragments%ke_orbit - fragments%ke_spin - ! The following ensures that fval = 0 is a local minimum, which is what the BFGS method is searching for - fval = (keo / (2 * ke_radial))**2 - + + associate(impactors => collision_system%impactors, nfrag => collision_system%fragments%nbody) + select type(fragments => collision_system%fragments) + class is (fraggle_fragments) + allocate(v_shift, mold=fragments%vb) + v_shift(:,:) = fraggle_util_vmag_to_vb(v_r_mag_input, fragments%v_r_unit, fragments%v_t_mag, fragments%v_t_unit, fragments%mass, impactors%vbcom) + !$omp do simd firstprivate(fragments) + do i = 1,fragments%nbody + rotmag2 = fragments%rot(1,i)**2 + fragments%rot(2,i)**2 + fragments%rot(3,i)**2 + vmag2 = v_shift(1,i)**2 + v_shift(2,i)**2 + v_shift(3,i)**2 + kearr(i) = fragments%mass(i) * (fragments%Ip(3, i) * fragments%radius(i)**2 * rotmag2 + vmag2) + end do + !$omp end do simd + keo = 2 * fragments%ke_budget - sum(kearr(:)) + ke_radial = fragments%ke_budget - fragments%ke_orbit - fragments%ke_spin + ! The following ensures that fval = 0 is a local minimum, which is what the BFGS method is searching for + fval = (keo / (2 * ke_radial))**2 + + end select + end associate + return end function radial_objective_function diff --git a/src/fraggle/fraggle_io.f90 b/src/fraggle/fraggle_io.f90 index c55db4173..e13c46472 100644 --- a/src/fraggle/fraggle_io.f90 +++ b/src/fraggle/fraggle_io.f90 @@ -12,40 +12,40 @@ contains - module subroutine fraggle_io_log_regime(impactors, fragments) + module subroutine fraggle_io_log_regime(collision_system) !! author: David A. Minton !! !! Writes a log of the results of the collisional regime determination implicit none ! Arguments - class(collision_impactors), intent(in) :: impactors !! Fraggle collider system object - class(fraggle_fragments), intent(in) :: fragments !! Fraggle fragment object + class(fraggle_system), intent(inout) :: collision_system !! Fraggle collision system object ! Internals character(STRMAX) :: errmsg - open(unit=LUN, file=FRAGGLE_LOG_OUT, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) - write(LUN, *, err = 667, iomsg = errmsg) - write(LUN, *) "--------------------------------------------------------------------" - write(LUN, *) " Fraggle collisional regime determination results" - write(LUN, *) "--------------------------------------------------------------------" - write(LUN, *) "True number of impactors : ",impactors%ncoll - write(LUN, *) "Index list of true impactors : ",impactors%idx(1:impactors%ncoll) - select case(impactors%regime) - case(COLLRESOLVE_REGIME_MERGE) - write(LUN, *) "Regime: Merge" - case(COLLRESOLVE_REGIME_DISRUPTION) - write(LUN, *) "Regime: Disruption" - case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - write(LUN, *) "Regime: Supercatastrophic disruption" - case(COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - write(LUN, *) "Regime: Graze and merge" - case(COLLRESOLVE_REGIME_HIT_AND_RUN) - write(LUN, *) "Regime: Hit and run" - end select - write(LUN, *) "Energy loss : ", impactors%Qloss - write(LUN, *) "--------------------------------------------------------------------" - close(LUN) - + associate(fragments => collision_system%fragments, impactors => collision_system%impactors) + open(unit=LUN, file=FRAGGLE_LOG_OUT, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) + write(LUN, *, err = 667, iomsg = errmsg) + write(LUN, *) "--------------------------------------------------------------------" + write(LUN, *) " Fraggle collisional regime determination results" + write(LUN, *) "--------------------------------------------------------------------" + write(LUN, *) "True number of impactors : ",impactors%ncoll + write(LUN, *) "Index list of true impactors : ",impactors%idx(1:impactors%ncoll) + select case(impactors%regime) + case(COLLRESOLVE_REGIME_MERGE) + write(LUN, *) "Regime: Merge" + case(COLLRESOLVE_REGIME_DISRUPTION) + write(LUN, *) "Regime: Disruption" + case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + write(LUN, *) "Regime: Supercatastrophic disruption" + case(COLLRESOLVE_REGIME_GRAZE_AND_MERGE) + write(LUN, *) "Regime: Graze and merge" + case(COLLRESOLVE_REGIME_HIT_AND_RUN) + write(LUN, *) "Regime: Hit and run" + end select + write(LUN, *) "Energy loss : ", impactors%Qloss + write(LUN, *) "--------------------------------------------------------------------" + close(LUN) + end associate return 667 continue write(*,*) "Error writing Fraggle regime information to log file: " // trim(adjustl(errmsg)) diff --git a/src/fraggle/fraggle_regime.f90 b/src/fraggle/fraggle_regime.f90 index 2873d236d..09015c66c 100644 --- a/src/fraggle/fraggle_regime.f90 +++ b/src/fraggle/fraggle_regime.f90 @@ -62,22 +62,22 @@ module subroutine encounter_regime_impactors(self, fragments, system, param) x1_si(:), x2_si(:), v1_si(:), v2_si(:), density_si(jtarg), density_si(jproj), & min_mfrag_si, impactors%regime, mlr, mslr, impactors%Qloss) - fragments%mass_dist(1) = min(max(mlr, 0.0_DP), mtot) - fragments%mass_dist(2) = min(max(mslr, 0.0_DP), mtot) - fragments%mass_dist(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot) + impactors%mass_dist(1) = min(max(mlr, 0.0_DP), mtot) + impactors%mass_dist(2) = min(max(mslr, 0.0_DP), mtot) + impactors%mass_dist(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot) ! Find the center of mass of the collisional system fragments%mtot = sum(impactors%mass(:)) - fragments%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / fragments%mtot - fragments%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / fragments%mtot + impactors%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / fragments%mtot + impactors%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / fragments%mtot ! Find the point of impact between the two bodies runit(:) = impactors%rb(:,2) - impactors%rb(:,1) runit(:) = runit(:) / (.mag. runit(:)) - fragments%rbimp(:) = impactors%rb(:,1) + impactors%radius(1) * runit(:) + impactors%rbimp(:) = impactors%rb(:,1) + impactors%radius(1) * runit(:) ! Convert quantities back to the system units and save them into the fragment system - fragments%mass_dist(:) = (fragments%mass_dist(:) / param%MU2KG) + impactors%mass_dist(:) = (impactors%mass_dist(:) / param%MU2KG) impactors%Qloss = impactors%Qloss * (param%TU2S / param%DU2M)**2 / param%MU2KG call fraggle_io_log_regime(impactors, fragments) diff --git a/src/fraggle/fraggle_set.f90 b/src/fraggle/fraggle_set.f90 index 0da950206..20770a1d0 100644 --- a/src/fraggle/fraggle_set.f90 +++ b/src/fraggle/fraggle_set.f90 @@ -11,31 +11,34 @@ use swiftest contains - module subroutine fraggle_set_budgets_fragments(self) + module subroutine fraggle_set_budgets(self) !! author: David A. Minton !! !! Sets the energy and momentum budgets of the fragments based on the collider values and the before/after values of energy and momentum implicit none ! Arguments - class(collision_system), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_system), intent(inout) :: self !! Fraggle collision system object ! Internals real(DP) :: dEtot real(DP), dimension(NDIM) :: dL - associate(impactors => self%impactors, fragments => self%fragments) + associate(impactors => self%impactors) + select type(fragments => self%fragments) + class is (fraggle_fragments) dEtot = self%Etot(2) - self%Etot(1) dL(:) = self%Ltot(:,2) - self%Ltot(:,1) fragments%L_budget(:) = -dL(:) - fragments%ke_budget = -(dEtot - 0.5_DP * fragments%mtot * dot_product(fragments%vbcom(:), fragments%vbcom(:))) - impactors%Qloss + fragments%ke_budget = -(dEtot - 0.5_DP * fragments%mtot * dot_product(impactors%vbcom(:), impactors%vbcom(:))) - impactors%Qloss + end select end associate return - end subroutine fraggle_set_budgets_fragments + end subroutine fraggle_set_budgets - module subroutine fraggle_set_mass_dist_fragments(self, impactors, param) + module subroutine fraggle_set_mass_dist(self, param) !! author: David A. Minton !! !! Sets the mass of fragments based on the mass distribution returned by the regime calculation. @@ -43,9 +46,8 @@ module subroutine fraggle_set_mass_dist_fragments(self, impactors, param) !! implicit none ! Arguments - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object - class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters + class(fraggle_system), intent(inout) :: self !! Fraggle collision system object + class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters ! Internals integer(I4B) :: i, jproj, jtarg, nfrag, istart real(DP), dimension(2) :: volume @@ -59,7 +61,9 @@ module subroutine fraggle_set_mass_dist_fragments(self, impactors, param) integer(I4B), parameter :: iMslr = 2 integer(I4B), parameter :: iMrem = 3 - associate(fragments => self) + associate(impactors => self%impactors) + select type(fragments => self%fragments) + class is (fraggle_fragments) ! Get mass weighted mean of Ip and density volume(1:2) = 4._DP / 3._DP * PI * impactors%radius(1:2)**3 Ip_avg(:) = (impactors%mass(1) * impactors%Ip(:,1) + impactors%mass(2) * impactors%Ip(:,2)) / fragments%mtot @@ -92,9 +96,9 @@ module subroutine fraggle_set_mass_dist_fragments(self, impactors, param) end select i = iMrem - mremaining = fragments%mass_dist(iMrem) + mremaining = impactors%mass_dist(iMrem) do while (i <= nfrag) - mfrag = (1 + i - iMslr)**(-3._DP / BETA) * fragments%mass_dist(iMslr) + mfrag = (1 + i - iMslr)**(-3._DP / BETA) * impactors%mass_dist(iMslr) if (mremaining - mfrag < 0.0_DP) exit mremaining = mremaining - mfrag i = i + 1 @@ -104,9 +108,9 @@ module subroutine fraggle_set_mass_dist_fragments(self, impactors, param) call fragments%setup(nfrag, param) case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) call fragments%setup(1, param) - fragments%mass(1) = fragments%mass_dist(1) + fragments%mass(1) = impactors%mass_dist(1) fragments%radius(1) = impactors%radius(jtarg) - fragments%density(1) = fragments%mass_dist(1) / volume(jtarg) + fragments%density(1) = impactors%mass_dist(1) / volume(jtarg) if (param%lrotation) fragments%Ip(:, 1) = impactors%Ip(:,1) return case default @@ -114,13 +118,13 @@ module subroutine fraggle_set_mass_dist_fragments(self, impactors, param) end select ! Make the first two bins the same as the Mlr and Mslr values that came from encounter_regime - fragments%mass(1) = fragments%mass_dist(iMlr) - fragments%mass(2) = fragments%mass_dist(iMslr) + fragments%mass(1) = impactors%mass_dist(iMlr) + fragments%mass(2) = impactors%mass_dist(iMslr) ! Distribute the remaining mass the 3:nfrag bodies following the model SFD given by slope BETA - mremaining = fragments%mass_dist(iMrem) + mremaining = impactors%mass_dist(iMrem) do i = iMrem, nfrag - mfrag = (1 + i - iMslr)**(-3._DP / BETA) * fragments%mass_dist(iMslr) + mfrag = (1 + i - iMslr)**(-3._DP / BETA) * impactors%mass_dist(iMslr) fragments%mass(i) = mfrag mremaining = mremaining - mfrag end do @@ -142,7 +146,7 @@ module subroutine fraggle_set_mass_dist_fragments(self, impactors, param) select case(impactors%regime) case(COLLRESOLVE_REGIME_HIT_AND_RUN) ! The hit and run case always preserves the largest body intact, so there is no need to recompute the physical properties of the first fragment fragments%radius(1) = impactors%radius(jtarg) - fragments%density(1) = fragments%mass_dist(iMlr) / volume(jtarg) + fragments%density(1) = impactors%mass_dist(iMlr) / volume(jtarg) fragments%Ip(:, 1) = impactors%Ip(:,1) istart = 2 case default @@ -153,55 +157,53 @@ module subroutine fraggle_set_mass_dist_fragments(self, impactors, param) do i = istart, nfrag fragments%Ip(:, i) = Ip_avg(:) end do - + end select end associate return - end subroutine fraggle_set_mass_dist_fragments - + end subroutine fraggle_set_mass_dist - module subroutine fraggle_set_natural_scale_factors(self, impactors) + module subroutine fraggle_set_natural_scale_factors(self) !! 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 ! Arguments - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object + class(fraggle_system), intent(inout) :: self !! Fraggle collision system object ! Internals integer(I4B) :: i - associate(fragments => self) + associate(collision_system => self, fragments => self%fragments, impactors => self%impactors) ! Set scale factors - fragments%Escale = 0.5_DP * (impactors%mass(1) * dot_product(impactors%vb(:,1), impactors%vb(:,1)) & - + impactors%mass(2) * dot_product(impactors%vb(:,2), impactors%vb(:,2))) - fragments%dscale = sum(impactors%radius(:)) - fragments%mscale = fragments%mtot - fragments%vscale = sqrt(fragments%Escale / fragments%mscale) - fragments%tscale = fragments%dscale / fragments%vscale - fragments%Lscale = fragments%mscale * fragments%dscale * fragments%vscale + collision_system%Escale = 0.5_DP * ( impactors%mass(1) * dot_product(impactors%vb(:,1), impactors%vb(:,1)) & + + impactors%mass(2) * dot_product(impactors%vb(:,2), impactors%vb(:,2))) + collision_system%dscale = sum(impactors%radius(:)) + collision_system%mscale = fragments%mtot + collision_system%vscale = sqrt(collision_system%Escale / collision_system%mscale) + collision_system%tscale = collision_system%dscale / collision_system%vscale + collision_system%Lscale = collision_system%mscale * collision_system%dscale * collision_system%vscale ! Scale all dimensioned quantities of impactors and fragments - fragments%rbcom(:) = fragments%rbcom(:) / fragments%dscale - fragments%vbcom(:) = fragments%vbcom(:) / fragments%vscale - fragments%rbimp(:) = fragments%rbimp(:) / fragments%dscale - impactors%rb(:,:) = impactors%rb(:,:) / fragments%dscale - impactors%vb(:,:) = impactors%vb(:,:) / fragments%vscale - impactors%mass(:) = impactors%mass(:) / fragments%mscale - impactors%radius(:) = impactors%radius(:) / fragments%dscale - impactors%Lspin(:,:) = impactors%Lspin(:,:) / fragments%Lscale - impactors%Lorbit(:,:) = impactors%Lorbit(:,:) / fragments%Lscale + impactors%rbcom(:) = impactors%rbcom(:) / collision_system%dscale + impactors%vbcom(:) = impactors%vbcom(:) / collision_system%vscale + impactors%rbimp(:) = impactors%rbimp(:) / collision_system%dscale + impactors%rb(:,:) = impactors%rb(:,:) / collision_system%dscale + impactors%vb(:,:) = impactors%vb(:,:) / collision_system%vscale + impactors%mass(:) = impactors%mass(:) / collision_system%mscale + impactors%radius(:) = impactors%radius(:) / collision_system%dscale + impactors%Lspin(:,:) = impactors%Lspin(:,:) / collision_system%Lscale + impactors%Lorbit(:,:) = impactors%Lorbit(:,:) / collision_system%Lscale do i = 1, 2 impactors%rot(:,i) = impactors%Lspin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i)) end do - fragments%mtot = fragments%mtot / fragments%mscale - fragments%mass = fragments%mass / fragments%mscale - fragments%radius = fragments%radius / fragments%dscale - impactors%Qloss = impactors%Qloss / fragments%Escale + fragments%mtot = fragments%mtot / collision_system%mscale + fragments%mass = fragments%mass / collision_system%mscale + fragments%radius = fragments%radius / collision_system%dscale + impactors%Qloss = impactors%Qloss / collision_system%Escale end associate return @@ -215,8 +217,7 @@ module subroutine fraggle_set_original_scale_factors(self) use, intrinsic :: ieee_exceptions implicit none ! Arguments - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object + class(fraggle_system), intent(inout) :: self !! Fraggle fragment system object ! Internals integer(I4B) :: i logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes @@ -224,58 +225,50 @@ module subroutine fraggle_set_original_scale_factors(self) call ieee_get_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily call ieee_set_halting_mode(IEEE_ALL,.false.) - associate(fragments => self) + associate(collision_system => self, fragments => self%fragments, impactors => self%impactors) ! Restore scale factors - fragments%rbcom(:) = fragments%rbcom(:) * fragments%dscale - fragments%vbcom(:) = fragments%vbcom(:) * fragments%vscale - fragments%rbimp(:) = fragments%rbimp(:) * fragments%dscale + impactors%rbcom(:) = impactors%rbcom(:) * collision_system%dscale + impactors%vbcom(:) = impactors%vbcom(:) * collision_system%vscale + impactors%rbimp(:) = impactors%rbimp(:) * collision_system%dscale - impactors%mass = impactors%mass * fragments%mscale - impactors%radius = impactors%radius * fragments%dscale - impactors%rb = impactors%rb * fragments%dscale - impactors%vb = impactors%vb * fragments%vscale - impactors%Lspin = impactors%Lspin * fragments%Lscale + impactors%mass = impactors%mass * collision_system%mscale + impactors%radius = impactors%radius * collision_system%dscale + impactors%rb = impactors%rb * collision_system%dscale + impactors%vb = impactors%vb * collision_system%vscale + impactors%Lspin = impactors%Lspin * collision_system%Lscale do i = 1, 2 impactors%rot(:,i) = impactors%Lspin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i)) end do - fragments%mtot = fragments%mtot * fragments%mscale - fragments%mass = fragments%mass * fragments%mscale - fragments%radius = fragments%radius * fragments%dscale - fragments%rot = fragments%rot / fragments%tscale - fragments%rc = fragments%rc * fragments%dscale - fragments%vc = fragments%vc * fragments%vscale + fragments%mtot = fragments%mtot * collision_system%mscale + fragments%mass = fragments%mass * collision_system%mscale + fragments%radius = fragments%radius * collision_system%dscale + fragments%rot = fragments%rot / collision_system%tscale + fragments%rc = fragments%rc * collision_system%dscale + fragments%vc = fragments%vc * collision_system%vscale do i = 1, fragments%nbody - fragments%rb(:, i) = fragments%rc(:, i) + fragments%rbcom(:) - fragments%vb(:, i) = fragments%vc(:, i) + fragments%vbcom(:) + fragments%rb(:, i) = fragments%rc(:, i) + impactors%rbcom(:) + fragments%vb(:, i) = fragments%vc(:, i) + impactors%vbcom(:) end do - impactors%Qloss = impactors%Qloss * fragments%Escale + impactors%Qloss = impactors%Qloss * collision_system%Escale - fragments%Lorbit_before(:) = fragments%Lorbit_before * fragments%Lscale - fragments%Lspin_before(:) = fragments%Lspin_before * fragments%Lscale - fragments%Ltot_before(:) = fragments%Ltot_before * fragments%Lscale - fragments%ke_orbit_before = fragments%ke_orbit_before * fragments%Escale - fragments%ke_spin_before = fragments%ke_spin_before * fragments%Escale - fragments%pe_before = fragments%pe_before * fragments%Escale - fragments%Etot_before = fragments%Etot_before * fragments%Escale - - fragments%Lorbit_after(:) = fragments%Lorbit_after * fragments%Lscale - fragments%Lspin_after(:) = fragments%Lspin_after * fragments%Lscale - fragments%Ltot_after(:) = fragments%Ltot_after * fragments%Lscale - fragments%ke_orbit_after = fragments%ke_orbit_after * fragments%Escale - fragments%ke_spin_after = fragments%ke_spin_after * fragments%Escale - fragments%pe_after = fragments%pe_after * fragments%Escale - fragments%Etot_after = fragments%Etot_after * fragments%Escale + collision_system%Lorbit(:,:) = collision_system%Lorbit(:,:) * collision_system%Lscale + collision_system%Lspin(:,:) = collision_system%Lspin(:,:) * collision_system%Lscale + collision_system%Ltot(:,:) = collision_system%Ltot(:,:) * collision_system%Lscale + collision_system%ke_orbit(:) = collision_system%ke_orbit(:) * collision_system%Escale + collision_system%ke_spin(:) = collision_system%ke_spin(:) * collision_system%Escale + collision_system%pe(:) = collision_system%pe(:) * collision_system%Escale + collision_system%Etot(:) = collision_system%Etot(:) * collision_system%Escale - fragments%mscale = 1.0_DP - fragments%dscale = 1.0_DP - fragments%vscale = 1.0_DP - fragments%tscale = 1.0_DP - fragments%Lscale = 1.0_DP - fragments%Escale = 1.0_DP + collision_system%mscale = 1.0_DP + collision_system%dscale = 1.0_DP + collision_system%vscale = 1.0_DP + collision_system%tscale = 1.0_DP + collision_system%Lscale = 1.0_DP + collision_system%Escale = 1.0_DP end associate call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 4d1a203fb..d550b7166 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -83,22 +83,22 @@ module subroutine fraggle_util_get_angular_momentum(self) end subroutine fraggle_util_get_angular_momentum - module subroutine fraggle_util_construct_temporary_system(fragments, system, param, tmpsys, tmpparam) + module subroutine fraggle_util_construct_temporary_system(collision_system, nbody_system, param, tmpsys, tmpparam) !! Author: David A. Minton !! !! Constructs a temporary internal system consisting of active bodies and additional fragments. This internal temporary system is used to calculate system energy with and without fragments implicit none ! Arguments - class(fraggle_fragments), intent(in) :: fragments !! Fraggle fragment system object - class(swiftest_nbody_system), intent(in) :: system !! Original swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters - class(swiftest_nbody_system), allocatable, intent(out) :: tmpsys !! Output temporary swiftest nbody system object - class(swiftest_parameters), allocatable, intent(out) :: tmpparam !! Output temporary configuration run parameters + class(fraggle_system), intent(inout) :: collision_system !! Fraggle collision system object + class(swiftest_nbody_system), intent(in) :: nbody_system !! Original swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(swiftest_nbody_system), allocatable, intent(out) :: tmpsys !! Output temporary swiftest nbody system object + class(swiftest_parameters), allocatable, intent(out) :: tmpparam !! Output temporary configuration run parameters ! Internals logical, dimension(:), allocatable :: linclude integer(I4B) :: npl_tot - associate(nfrag => fragments%nbody, pl => system%pl, npl => system%pl%nbody, cb => system%cb) + associate(fragments => collision_system%fragments, nfrag => collision_system%fragments%nbody, pl => nbody_system%pl, npl => nbody_system%pl%nbody, cb => nbody_system%cb) ! Set up a new system based on the original if (allocated(tmpparam)) deallocate(tmpparam) if (allocated(tmpsys)) deallocate(tmpsys) @@ -123,7 +123,7 @@ module subroutine fraggle_util_construct_temporary_system(fragments, system, par call tmpsys%pl%fill(pl, linclude) ! Scale the temporary system to the natural units of the current Fraggle calculation - call tmpsys%rescale(tmpparam, fragments%mscale, fragments%dscale, fragments%tscale) + call tmpsys%rescale(tmpparam, collision_system%mscale, collision_system%dscale, collision_system%tscale) end associate diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index b8a558af8..eff70b1f8 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -34,20 +34,10 @@ module fraggle_classes real(DP), dimension(NDIM) :: Lspin !! Spin angular momentum vector of all fragments real(DP) :: ke_orbit !! Orbital kinetic energy of all fragments real(DP) :: ke_spin !! Spin kinetic energy of all fragments - real(DP) :: ke_budget !! Kinetic energy budget for computing quanities + real(DP) :: ke_budget !! Kinetic energy budget for computing fragment trajectories + real(DP), dimension(NDIM) :: L_budget !! Angular momentum budget for computing fragment trajectories - ! Scale factors used to scale dimensioned quantities to a more "natural" system where important quantities (like kinetic energy, momentum) are of order ~1 - real(DP) :: dscale = 1.0_DP !! Distance dimension scale factor - real(DP) :: mscale = 1.0_DP !! Mass scale factor - real(DP) :: tscale = 1.0_DP !! Time scale factor - real(DP) :: vscale = 1.0_DP !! Velocity scale factor (a convenience unit that is derived from dscale and tscale) - real(DP) :: Escale = 1.0_DP !! Energy scale factor (a convenience unit that is derived from dscale, tscale, and mscale) - real(DP) :: Lscale = 1.0_DP !! Angular momentum scale factor (a convenience unit that is derived from dscale, tscale, and mscale) contains - procedure :: set_budgets => fraggle_set_budgets_fragments !! Sets the energy and momentum budgets of the fragments based on the collider value - procedure :: set_mass_dist => fraggle_set_mass_dist_fragments !! Sets the distribution of mass among the fragments depending on the regime type - procedure :: set_natural_scale => fraggle_set_natural_scale_factors !! Scales dimenional quantities to ~O(1) with respect to the collisional system. - procedure :: set_original_scale => fraggle_set_original_scale_factors !! Restores dimenional quantities back to the original system units procedure :: setup => fraggle_setup_fragments !! Allocates arrays for n fragments in a Fraggle system. Passing n = 0 deallocates all arrays. procedure :: reset => fraggle_setup_reset_fragments !! Resets all position and velocity-dependent fragment quantities in order to do a fresh calculation (does not reset mass, radius, or other values that get set prior to the call to fraggle_generate) procedure :: get_angular_momentum => fraggle_util_get_angular_momentum !! Calcualtes the current angular momentum of the fragments @@ -58,9 +48,20 @@ module fraggle_classes end type fraggle_fragments type, extends(collision_system) :: fraggle_system + ! Scale factors used to scale dimensioned quantities to a more "natural" system where important quantities (like kinetic energy, momentum) are of order ~1 + real(DP) :: dscale = 1.0_DP !! Distance dimension scale factor + real(DP) :: mscale = 1.0_DP !! Mass scale factor + real(DP) :: tscale = 1.0_DP !! Time scale factor + real(DP) :: vscale = 1.0_DP !! Velocity scale factor (a convenience unit that is derived from dscale and tscale) + real(DP) :: Escale = 1.0_DP !! Energy scale factor (a convenience unit that is derived from dscale, tscale, and mscale) + real(DP) :: Lscale = 1.0_DP !! Angular momentum scale factor (a convenience unit that is derived from dscale, tscale, and mscale) contains - procedure :: generate_fragments => fraggle_generate_fragments !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. - final :: fraggle_util_final_system !! Finalizer will deallocate all allocatables + procedure :: generate_fragments => fraggle_generate_fragments !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. + procedure :: set_budgets => fraggle_set_budgets !! Sets the energy and momentum budgets of the fragments based on the collider value + procedure :: set_mass_dist => fraggle_set_mass_dist !! Sets the distribution of mass among the fragments depending on the regime type + procedure :: set_natural_scale => fraggle_set_natural_scale_factors !! Scales dimenional quantities to ~O(1) with respect to the collisional system. + procedure :: set_original_scale => fraggle_set_original_scale_factors !! Restores dimenional quantities back to the original system units + final :: fraggle_util_final_system !! Finalizer will deallocate all allocatables end type fraggle_system @@ -74,34 +75,30 @@ module subroutine fraggle_generate_fragments(self, system, param, lfailure) logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? end subroutine fraggle_generate_fragments - module subroutine fraggle_io_log_regime(impactors, fragments) + module subroutine fraggle_io_log_regime(collision_system) implicit none - class(collision_impactors), intent(in) :: impactors - class(fraggle_fragments), intent(in) :: fragments + class(fraggle_system), intent(inout) :: collision_system !! Fraggle collision system object end subroutine fraggle_io_log_regime - module subroutine fraggle_set_budgets_fragments(self) + module subroutine fraggle_set_budgets(self) implicit none - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - end subroutine fraggle_set_budgets_fragments + class(fraggle_system), intent(inout) :: self !! Fraggle collision system object + end subroutine fraggle_set_budgets - module subroutine fraggle_set_mass_dist_fragments(self, impactors, param) + module subroutine fraggle_set_mass_dist(self, param) implicit none - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object - class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters - end subroutine fraggle_set_mass_dist_fragments + class(fraggle_system), intent(inout) :: self !! Fraggle collision system object + class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters + end subroutine fraggle_set_mass_dist - module subroutine fraggle_set_natural_scale_factors(self, impactors) + module subroutine fraggle_set_natural_scale_factors(self) implicit none - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object + class(fraggle_system), intent(inout) :: self !! Fraggle collision system object end subroutine fraggle_set_natural_scale_factors - module subroutine fraggle_set_original_scale_factors(self, impactors) + module subroutine fraggle_set_original_scale_factors(self) implicit none - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object + class(fraggle_system), intent(inout) :: self !! Fraggle collision system object end subroutine fraggle_set_original_scale_factors module subroutine fraggle_setup_fragments(self, n, param) @@ -130,14 +127,14 @@ module subroutine fraggle_util_get_angular_momentum(self) class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object end subroutine fraggle_util_get_angular_momentum - module subroutine fraggle_util_construct_temporary_system(fragments, system, param, tmpsys, tmpparam) + module subroutine fraggle_util_construct_temporary_system(collision_system, nbody_system, param, tmpsys, tmpparam) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none - class(fraggle_fragments), intent(in) :: fragments !! Fraggle fragment system object - class(swiftest_nbody_system), intent(in) :: system !! Original swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters - class(swiftest_nbody_system), allocatable, intent(out) :: tmpsys !! Output temporary swiftest nbody system object - class(swiftest_parameters), allocatable, intent(out) :: tmpparam !! Output temporary configuration run parameters + class(fraggle_system), intent(inout) :: collision_system !! Fraggle collision system object + class(swiftest_nbody_system), intent(in) :: nbody_system !! Original swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(swiftest_nbody_system), allocatable, intent(out) :: tmpsys !! Output temporary swiftest nbody system object + class(swiftest_parameters), allocatable, intent(out) :: tmpparam !! Output temporary configuration run parameters end subroutine fraggle_util_construct_temporary_system module subroutine fraggle_util_dealloc_fragments(self) diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 16fd0d7bd..663e40b55 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -118,7 +118,7 @@ module function symba_collision_casehitandrun(system, param, t) result(status) jproj = 1 end if - if (fragments%mass_dist(2) > 0.9_DP * impactors%mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched + if (impactors%mass_dist(2) > 0.9_DP * impactors%mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched call io_log_one_message(FRAGGLE_LOG_OUT, "Pure hit and run. No new fragments generated.") nfrag = 0 lpure = .true. @@ -202,8 +202,8 @@ module function symba_collision_casemerge(system, param, t) result(status) ibiggest = impactors%idx(maxloc(pl%Gmass(impactors%idx(:)), dim=1)) fragments%id(1) = pl%id(ibiggest) - fragments%rb(:,1) = fragments%rbcom(:) - fragments%vb(:,1) = fragments%vbcom(:) + fragments%rb(:,1) = impactors%rbcom(:) + fragments%vb(:,1) = impactors%vbcom(:) if (param%lrotation) then ! Conserve angular momentum by putting pre-impact orbital momentum into spin of the new body @@ -939,11 +939,11 @@ subroutine symba_resolve_collision(plplcollision_list , system, param, t) associate(fragments => system%fragments, impactors => system%impactors) impactors%regime = COLLRESOLVE_REGIME_MERGE fragments%mtot = sum(impactors%mass(:)) - fragments%mass_dist(1) = fragments%mtot - fragments%mass_dist(2) = 0.0_DP - fragments%mass_dist(3) = 0.0_DP - fragments%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / fragments%mtot - fragments%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / fragments%mtot + impactors%mass_dist(1) = fragments%mtot + impactors%mass_dist(2) = 0.0_DP + impactors%mass_dist(3) = 0.0_DP + impactors%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / fragments%mtot + impactors%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / fragments%mtot end associate end if