From 0a30daf509a52014e2aa4d341b33cf0b17d70f83 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 2 Sep 2021 17:34:03 -0400 Subject: [PATCH] Major restructuring of the fragmentation model. I know call it Fraggle and it is fully object oriented. --- Makefile | 7 +- Makefile.Defines | 4 +- src/fraggle/fraggle_generate.f90 | 525 +++++++++++ src/fraggle/fraggle_placeholder.f90 | 44 + src/fraggle/fraggle_regime.f90 | 356 ++++++++ src/fraggle/fraggle_set.f90 | 295 +++++++ src/fraggle/fraggle_setup.f90 | 78 ++ src/fraggle/fraggle_util.f90 | 316 +++++++ src/fragmentation/fragmentation.f90 | 1266 --------------------------- src/helio/helio_setup.f90 | 2 + src/helio/helio_util.f90 | 21 - src/modules/fraggle_classes.f90 | 250 ++++++ src/modules/helio_classes.f90 | 9 +- src/modules/swiftest.f90 | 1 + src/modules/swiftest_classes.f90 | 29 +- src/modules/symba_classes.f90 | 75 +- src/modules/whm_classes.f90 | 1 - src/setup/setup.f90 | 4 +- src/symba/symba_collision.f90 | 598 +++++-------- src/symba/symba_util.f90 | 3 +- src/whm/whm_setup.f90 | 2 + src/whm/whm_util.f90 | 16 - 22 files changed, 2124 insertions(+), 1778 deletions(-) create mode 100644 src/fraggle/fraggle_generate.f90 create mode 100644 src/fraggle/fraggle_placeholder.f90 create mode 100644 src/fraggle/fraggle_regime.f90 create mode 100644 src/fraggle/fraggle_set.f90 create mode 100644 src/fraggle/fraggle_setup.f90 create mode 100644 src/fraggle/fraggle_util.f90 delete mode 100644 src/fragmentation/fragmentation.f90 delete mode 100644 src/helio/helio_util.f90 create mode 100644 src/modules/fraggle_classes.f90 diff --git a/Makefile b/Makefile index 8d8955eb4..fb555345f 100644 --- a/Makefile +++ b/Makefile @@ -45,9 +45,10 @@ #****************************************************************************** SWIFTEST_MODULES = swiftest_globals.f90 \ + swiftest_operators.f90 \ lambda_function.f90\ swiftest_classes.f90 \ - swiftest_operators.f90 \ + fraggle_classes.f90 \ whm_classes.f90 \ rmvs_classes.f90 \ helio_classes.f90 \ @@ -91,7 +92,7 @@ lib: ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ ln -s $(SWIFTEST_HOME)/Makefile .; \ make libdir - cd $(SWIFTEST_HOME)/src/fragmentation; \ + cd $(SWIFTEST_HOME)/src/fraggle; \ rm -f Makefile.Defines Makefile; \ ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ ln -s $(SWIFTEST_HOME)/Makefile .; \ @@ -193,7 +194,7 @@ clean: cd $(SWIFTEST_HOME)/src/modules; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/discard; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/drift; rm -f Makefile.Defines Makefile *.gc* - cd $(SWIFTEST_HOME)/src/fragmentation; rm -f Makefile.Defines Makefile *.gc* + cd $(SWIFTEST_HOME)/src/fraggle; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/gr; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/helio; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/io; rm -f Makefile.Defines Makefile *.gc* diff --git a/Makefile.Defines b/Makefile.Defines index e48cb9d07..36de2cdb7 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -67,13 +67,13 @@ GWARNINGS = -Wall -Warray-bounds -Wimplicit-interface -Wextra -Warray-temporari GPRODUCTION = -O3 -ffree-line-length-none $(GPAR) #FFLAGS = $(IDEBUG) $(HEAPARR) $(SIMDVEC) $(PAR) -FFLAGS = $(IPRODUCTION) -g -traceback $(OPTREPORT) +FFLAGS = $(IPRODUCTION) $(OPTREPORT) FORTRAN = ifort #AR = xiar #FORTRAN = gfortran #FFLAGS = $(GDEBUG) $(GMEM) $(GPAR) -#FFLAGS = $(GPRODUCTION) -g -fcheck=all -Wall -fbacktrace +#FFLAGS = $(GPRODUCTION) -g -fbacktrace #-fcheck=all #-Wall AR = ar # DO NOT include in CFLAGS the "-c" option to compile object only diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 new file mode 100644 index 000000000..01fc3cf18 --- /dev/null +++ b/src/fraggle/fraggle_generate.f90 @@ -0,0 +1,525 @@ +submodule(fraggle_classes) s_fraggle_generate + use swiftest + + integer(I4B), parameter :: NFRAG_MIN = 7 !! The minimum allowable number of fragments (set to 6 because that's how many unknowns are needed in the tangential velocity calculation) + real(DP), parameter :: F_SPIN_FIRST = 0.05_DP !! The initial try value of the fraction of energy or momenum in spin (whichever has the lowest kinetic energy) + real(DP), parameter :: FRAGGLE_LTOL = 10 * epsilon(1.0_DP) + real(DP), parameter :: FRAGGLE_ETOL = 1e-8_DP + +contains + + module subroutine fraggle_generate_fragments(self, colliders, system, param, lfailure) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. + use, intrinsic :: ieee_exceptions + implicit none + ! Arguments + class(fraggle_fragments), intent(inout) :: self !! Fraggle system object the outputs will be the fragmentation + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object containing the two-body equivalent values of the colliding bodies + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? + ! Internals + integer(I4B) :: i + integer(I4B) :: try + real(DP) :: r_max_start, f_spin, dEtot, dLmag + integer(I4B), parameter :: MAXTRY = 3000 + character(len=*), parameter :: fmtlabel = "(A14,10(ES11.4,1X,:))" + logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes + logical :: lk_plpl + + call ieee_get_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily + fpe_quiet_modes(:) = .false. + call ieee_set_halting_mode(IEEE_ALL,fpe_quiet_modes) + + associate(frag => self, nfrag => self%nbody, pl => system%pl) + + if (nfrag < NFRAG_MIN) then + write(*,*) "Fraggle needs at least ",NFRAG_MIN," fragments, but only ",nfrag," were given." + lfailure = .true. + return + end if + f_spin = F_SPIN_FIRST + + lk_plpl = allocated(pl%k_plpl) + if (lk_plpl) deallocate(pl%k_plpl) + + call frag%set_natural_scale(colliders) + + call frag%reset() + + ! Calculate the initial energy of the system without the collisional family + call frag%get_energy_and_momentum(colliders, 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 * norm2(colliders%xb(:,2) - colliders%xb(:,1)) + try = 1 + lfailure = .false. + do while (try < MAXTRY) + lfailure = .false. + + call fraggle_generate_pos_vec(frag, colliders, r_max_start) + call frag%set_coordinate_system(colliders) + + ! Initial velocity guess will be the barycentric velocity of the colliding system. + ! This ensures that our budgets are based on the residual velocities needed in the + ! collisional frame + do concurrent (i = 1:nfrag) + frag%vb(:, i) = frag%vbcom(:) + end do + + call frag%get_energy_and_momentum(colliders, system, param, lbefore=.false.) + call frag%set_budgets(colliders) + + call fraggle_generate_spins(frag, colliders, f_spin, lfailure) + if (.not.lfailure) then + call fraggle_generate_tan_vel(frag, colliders, lfailure) + if (.not. lfailure) then + call fraggle_generate_rad_vel(frag, colliders, lfailure) + ! if (lfailure) write(*,*) 'Failed to find radial velocities' + if (.not.lfailure) then + call frag%get_energy_and_momentum(colliders, system, param, lbefore=.false.) + dEtot = frag%Etot_after - frag%Etot_before + dLmag = .mag. (frag%Ltot_after(:) - frag%Ltot_before(:)) + + if ((abs(dEtot + frag%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) then + ! write(*,*) 'Failed due to high energy error: ',dEtot, abs(dEtot + frag%Qloss) / FRAGGLE_ETOL + lfailure = .true. + else if ((abs(dLmag) / (.mag.frag%Ltot_before)) > FRAGGLE_LTOL) then + ! write(*,*) 'Failed due to high angular momentum error: ', dLmag / (.mag.frag%Ltot_before(:)) + lfailure = .true. + end if + end if + end if + end if + + if (.not.lfailure) exit + call frag%restructure(colliders, try, f_spin, r_max_start) + call frag%reset() + try = try + 1 + end do + + ! write(*, "(' -------------------------------------------------------------------------------------')") + ! write(*, "(' Final diagnostic')") + ! write(*, "(' -------------------------------------------------------------------------------------')") + if (lfailure) then + write(*,*) "Fraggle fragment generation failed after: ",try," tries" + else + write(*,*) "Fraggle fragment generation succeeded after: ",try," tries" + ! write(*, "(' dL_tot should be very small' )") + ! write(*,fmtlabel) ' dL_tot |', (.mag.(frag%Ltot_after(:) - frag%Ltot_before(:))) / (.mag.frag%Ltot_before(:)) + ! write(*, "(' dE_tot should be negative and equal to Qloss' )") + ! write(*,fmtlabel) ' dE_tot |', (frag%Etot_after - frag%Etot_before) / abs(frag%Etot_before) + ! write(*,fmtlabel) ' Qloss |', -frag%Qloss / abs(frag%Etot_before) + ! write(*,fmtlabel) ' dE - Qloss |', (frag%Etot_after - frag%Etot_before + frag%Qloss) / abs(frag%Etot_before) + end if + ! write(*, "(' -------------------------------------------------------------------------------------')") + call frag%set_original_scale(colliders) + + ! Restore the big array + if (lk_plpl) call pl%index(param) + end associate + + call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily + + return + end subroutine fraggle_generate_fragments + + + subroutine fraggle_generate_pos_vec(frag, colliders, 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 + !! 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 + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragment system object + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider 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 + logical, dimension(:), allocatable :: loverlap + integer(I4B) :: i, j + + associate(nfrag => frag%nbody) + allocate(loverlap(nfrag)) + + ! Place the fragments into a region that is big enough that we should usually not have overlapping bodies + ! An overlapping bodies will collide in the next time step, so it's not a major problem if they do (it just slows the run down) + r_max = r_max_start + rad = sum(colliders%radius(:)) + + ! We will treat the first two fragments of the list as special cases. They get initialized the maximum distances apart along the original impactor distance vector. + ! This is done because in a regular disruption, the first body is the largest, the second the second largest, and the rest are smaller equal-mass fragments. + + call random_number(frag%x_coll(:,3:nfrag)) + loverlap(:) = .true. + do while (any(loverlap(3:nfrag))) + frag%x_coll(:, 1) = colliders%xb(:, 1) - frag%xbcom(:) + frag%x_coll(:, 2) = colliders%xb(:, 2) - frag%xbcom(:) + r_max = r_max + 0.1_DP * rad + do i = 3, nfrag + if (loverlap(i)) then + call random_number(frag%x_coll(:,i)) + frag%x_coll(:, i) = 2 * (frag%x_coll(:, i) - 0.5_DP) * r_max + end if + end do + loverlap(:) = .false. + do j = 1, nfrag + do i = j + 1, nfrag + dis = norm2(frag%x_coll(:,j) - frag%x_coll(:,i)) + loverlap(i) = loverlap(i) .or. (dis <= (frag%radius(i) + frag%radius(j))) + end do + end do + end do + call fraggle_util_shift_vector_to_origin(frag%mass, frag%x_coll) + call frag%set_coordinate_system(colliders) + + do i = 1, nfrag + frag%xb(:,i) = frag%x_coll(:,i) + frag%xbcom(:) + end do + + frag%xbcom(:) = 0.0_DP + do i = 1, nfrag + frag%xbcom(:) = frag%xbcom(:) + frag%mass(i) * frag%xb(:,i) + end do + frag%xbcom(:) = frag%xbcom(:) / frag%mtot + end associate + + return + end subroutine fraggle_generate_pos_vec + + + subroutine fraggle_generate_spins(frag, colliders, 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. + !! + !! 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) :: frag !! Fraggle fragment system object + class(fraggle_colliders), intent(in) :: colliders !! 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! + ! Internals + real(DP), dimension(NDIM) :: L_remainder, rot_L, rot_ke + integer(I4B) :: i + + associate(nfrag => frag%nbody) + lfailure = .false. + + ! Start the first two bodies with the same rotation as the original two impactors, then distribute the remaining angular momentum among the rest + frag%rot(:,1:2) = colliders%rot(:, :) + frag%rot(:,3:nfrag) = 0.0_DP + call frag%get_ang_mtm() + L_remainder(:) = frag%L_budget(:) - frag%L_spin(:) + + frag%ke_spin = 0.0_DP + do i = 1, nfrag + ! Convert a fraction (f_spin) of either the remaining angular momentum or kinetic energy budget into spin, whichever gives the smaller rotation so as not to blow any budgets + rot_ke(:) = sqrt(2 * f_spin * frag%ke_budget / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i))) * L_remainder(:) / norm2(L_remainder(:)) + rot_L(:) = f_spin * L_remainder(:) / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i)) + if (norm2(rot_ke) < norm2(rot_L)) then + frag%rot(:,i) = frag%rot(:, i) + rot_ke(:) + else + frag%rot(:, i) = frag%rot(:, i) + rot_L(:) + end if + frag%ke_spin = frag%ke_spin + frag%mass(i) * frag%Ip(3, i) * frag%radius(i)**2 * dot_product(frag%rot(:, i), frag%rot(:, i)) + end do + frag%ke_spin = 0.5_DP * frag%ke_spin + + lfailure = ((frag%ke_budget - frag%ke_spin - frag%ke_orbit) < 0.0_DP) + end associate + + return + end subroutine fraggle_generate_spins + + + subroutine fraggle_generate_tan_vel(frag, colliders, 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. + !! This procedure works in several stages, with a goal to solve the angular and linear momentum constraints on the fragments, while still leaving a positive balance of + !! our fragment kinetic energy (frag%ke_budget) that we can put into the radial velocity distribution. + !! + !! The first thing we'll try to do is solve for the tangential velocities of the first 6 fragments, using angular and linear momentum as constraints and an initial + !! tangential velocity distribution for the remaining bodies (if there are any) that distributes their angular momentum equally between them. + !! If that doesn't work and we blow our kinetic energy budget, we will attempt to find a tangential velocity distribution that minimizes the kinetic energy while + !! conserving momentum. + !! + !! 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) :: frag !! Fraggle fragment system object + class(fraggle_colliders), intent(in) :: colliders !! Fraggle collider system object + logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds + ! Internals + integer(I4B) :: i + real(DP), parameter :: TOL_MIN = 1e-1_DP ! This doesn't have to be very accurate, as we really just want a tangential velocity distribution with less kinetic energy than our initial guess. + real(DP), parameter :: TOL_INIT = 1e-14_DP + integer(I4B), parameter :: MAXLOOP = 10 + real(DP) :: tol + real(DP), dimension(:), allocatable :: v_t_initial + real(DP), dimension(frag%nbody) :: kefrag + type(lambda_obj) :: spinfunc + type(lambda_obj_err) :: objective_function + real(DP), dimension(NDIM) :: Li, L_remainder + + associate(nfrag => frag%nbody) + lfailure = .false. + + ! write(*,*) '***************************************************' + ! write(*,*) 'Original dis : ',norm2(x(:,2) - x(:,1)) + ! write(*,*) 'r_max : ',r_max + ! write(*,*) 'f_spin : ',f_spin + ! write(*,*) '***************************************************' + ! write(*,*) 'Energy balance so far: ' + ! write(*,*) 'frag%ke_budget : ',frag%ke_budget + ! write(*,*) 'ke_orbit_before: ',ke_orbit_before + ! write(*,*) 'ke_orbit_after : ',ke_orbit_after + ! write(*,*) 'ke_spin_before : ',ke_spin_before + ! write(*,*) 'ke_spin_after : ',ke_spin_after + ! write(*,*) 'pe_before : ',pe_before + ! write(*,*) 'pe_after : ',pe_after + ! write(*,*) 'Qloss : ',Qloss + ! write(*,*) '***************************************************' + + allocate(v_t_initial, mold=frag%v_t_mag) + v_t_initial(:) = 0.0_DP + frag%v_coll(:,:) = 0.0_DP + + ! Next we will solve for the tangential component of the velocities that both conserves linear momentum and uses the remaining angular momentum not used in spin. + ! This will be done using a linear solver that solves for the tangential velocities of the first 6 fragments, constrained by the linear and angular momentum vectors, + ! which is embedded in a non-linear minimizer that will adjust the tangential velocities of the remaining i>6 fragments to minimize kinetic energy for a given momentum solution + ! The initial conditions fed to the minimizer for the fragments will be the remaining angular momentum distributed between the fragments. + call frag%get_ang_mtm() + L_remainder(:) = frag%L_budget(:) - frag%L_spin(:) + do i = 1, nfrag + v_t_initial(i) = norm2(L_remainder(:)) / ((nfrag - i + 1) * frag%mass(i) * norm2(frag%x_coll(:,i))) + Li(:) = frag%mass(i) * (frag%x_coll(:,i) .cross. (v_t_initial(i) * frag%v_t_unit(:, i))) + L_remainder(:) = L_remainder(:) - Li(:) + end do + + ! Find the local kinetic energy minimum for the system that conserves linear and angular momentum + objective_function = lambda_obj(tangential_objective_function, lfailure) + + tol = TOL_INIT + do while(tol < TOL_MIN) + frag%v_t_mag(7:nfrag) = util_minimize_bfgs(objective_function, nfrag-6, v_t_initial(7:nfrag), tol, MAXLOOP, lfailure) + ! Now that the KE-minimized values of the i>6 fragments are found, calculate the momentum-conserving solution for tangential velociteis + v_t_initial(7:nfrag) = frag%v_t_mag(7:nfrag) + if (.not.lfailure) exit + tol = tol * 2_DP ! Keep increasing the tolerance until we converge on a solution + end do + frag%v_t_mag(1:nfrag) = solve_fragment_tan_vel(v_t_mag_input=v_t_initial(7:nfrag), lfailure=lfailure) + + ! Perform one final shift of the radial velocity vectors to align with the center of mass of the collisional system (the origin) + frag%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(frag%v_r_mag(1:nfrag), frag%v_r_unit(:,1:nfrag), frag%v_t_mag(1:nfrag), frag%v_t_unit(:,1:nfrag), frag%mass(1:nfrag), frag%vbcom(:)) + do concurrent (i = 1:nfrag) + frag%v_coll(:,i) = frag%vb(:,i) - frag%vbcom(:) + end do + + ! Now do a kinetic energy budget check to make sure we are still within the budget. + kefrag = 0.0_DP + do concurrent(i = 1:nfrag) + kefrag(i) = frag%mass(i) * dot_product(frag%vb(:, i), frag%vb(:, i)) + end do + frag%ke_orbit = 0.5_DP * sum(kefrag(:)) + + ! If we are over the energy budget, flag this as a failure so we can try again + lfailure = ((frag%ke_budget - frag%ke_spin - frag%ke_orbit) < 0.0_DP) + ! write(*,*) 'Tangential' + ! write(*,*) 'Failure? ',lfailure + ! call calculate_fragment_ang_mtm() + ! write(*,*) '|L_remainder| : ',.mag.(frag%L_budget(:) - L_frag_tot(:)) / Lmag_before + ! write(*,*) 'frag%ke_budget: ',frag%ke_budget + ! write(*,*) 'frag%ke_spin : ',frag%ke_spin + ! write(*,*) 'ke_tangential : ',frag%ke_orbit + ! write(*,*) 'ke_radial : ',frag%ke_budget - frag%ke_spin - frag%ke_orbit + end associate + + return + + contains + function solve_fragment_tan_vel(lfailure, v_t_mag_input) result(v_t_mag_output) + !! 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 + logical, intent(out) :: lfailure !! Error flag + real(DP), dimension(:), optional, intent(in) :: v_t_mag_input !! Unknown tangential velocities for fragments 7:nfrag + ! Internals + integer(I4B) :: i + ! Result + real(DP), dimension(:), allocatable :: 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 + + associate(nfrag => frag%nbody) + lfailure = .false. + ! 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 + 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) = frag%mass(i) * frag%v_t_unit(:, i) + A(4:6, i) = frag%mass(i) * frag%rmag(i) * (frag%v_r_unit(:, i) .cross. frag%v_t_unit(:, i)) + else if (present(v_t_mag_input)) then + vtmp(:) = v_t_mag_input(i - 6) * frag%v_t_unit(:, i) + L_lin_others(:) = L_lin_others(:) + frag%mass(i) * vtmp(:) + L(:) = frag%mass(i) * (frag%x_coll(:, i) .cross. vtmp(:)) + L_orb_others(:) = L_orb_others(:) + L(:) + end if + end do + b(1:3) = -L_lin_others(:) + b(4:6) = frag%L_budget(:) - frag%L_spin(:) - L_orb_others(:) + allocate(v_t_mag_output(nfrag)) + 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 + return + end function solve_fragment_tan_vel + + + function tangential_objective_function(v_t_mag_input, lfailure) result(fval) + !! Author: David A. Minton + !! + !! 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) :: v_t_mag_input !! Unknown tangential component of velocity vector set previously by angular momentum constraint + logical, intent(out) :: lfailure !! Error flag + ! Result + real(DP) :: fval + ! Internals + integer(I4B) :: i + real(DP), dimension(NDIM,frag%nbody) :: v_shift + real(DP), dimension(frag%nbody) :: v_t_new, kearr + real(DP) :: keo + + associate(nfrag => frag%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(frag%v_r_mag, frag%v_r_unit, v_t_new, frag%v_t_unit, frag%mass, frag%vbcom) + + kearr = 0.0_DP + do concurrent(i = 1:nfrag) + kearr(i) = frag%mass(i) * dot_product(v_shift(:, i), v_shift(:, i)) + end do + keo = 0.5_DP * sum(kearr(:)) + fval = keo + lfailure = .false. + end associate + + return + end function tangential_objective_function + + end subroutine fraggle_generate_tan_vel + + + subroutine fraggle_generate_rad_vel(frag, colliders, 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) :: frag !! Fraggle fragment system object + class(fraggle_colliders), intent(in) :: colliders !! Fraggle collider 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 + real(DP), parameter :: TOL_INIT = 1e-14_DP + integer(I4B), parameter :: MAXLOOP = 100 + real(DP) :: ke_radial, tol + integer(I4B) :: i, j + real(DP), dimension(:), allocatable :: v_r_initial, v_r_sigma + real(DP), dimension(:,:), allocatable :: v_r + type(lambda_obj) :: objective_function + + associate(nfrag => frag%nbody) + ! Set the "target" ke for the radial component + ke_radial = frag%ke_budget - frag%ke_spin - frag%ke_orbit + + allocate(v_r_initial, source=frag%v_r_mag) + ! Initialize radial velocity magnitudes with a random value that is approximately 10% of that found by distributing the kinetic energy equally + allocate(v_r_sigma, source=frag%v_r_mag) + call random_number(v_r_sigma(1:nfrag)) + v_r_sigma(1:nfrag) = sqrt(1.0_DP + 2 * (v_r_sigma(1:nfrag) - 0.5_DP) * 1e-4_DP) + v_r_initial(1:nfrag) = v_r_sigma(1:nfrag) * sqrt(abs(2 * ke_radial) / (frag%mass(1:nfrag) * nfrag)) + + ! Initialize the lambda function using a structure constructor that calls the init method + ! Minimize the ke objective function using the BFGS optimizer + objective_function = lambda_obj(radial_objective_function) + tol = TOL_INIT + do while(tol < TOL_MIN) + frag%v_r_mag = util_minimize_bfgs(objective_function, nfrag, v_r_initial, tol, MAXLOOP, lfailure) + if (.not.lfailure) exit + tol = tol * 2 ! Keep increasing the tolerance until we converge on a solution + v_r_initial(:) = frag%v_r_mag(:) + end do + + ! Shift the radial velocity vectors to align with the center of mass of the collisional system (the origin) + frag%ke_orbit = 0.0_DP + frag%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(frag%v_r_mag(1:nfrag), frag%v_r_unit(:,1:nfrag), frag%v_t_mag(1:nfrag), frag%v_t_unit(:,1:nfrag), frag%mass(1:nfrag), frag%vbcom(:)) + do i = 1, nfrag + frag%v_coll(:, i) = frag%vb(:, i) - frag%vbcom(:) + frag%ke_orbit = frag%ke_orbit + frag%mass(i) * dot_product(frag%vb(:, i), frag%vb(:, i)) + end do + frag%ke_orbit = 0.5_DP * frag%ke_orbit + + ! write(*,*) 'Radial' + ! write(*,*) 'Failure? ',lfailure + ! write(*,*) 'frag%ke_budget: ',frag%ke_budget + ! write(*,*) 'frag%ke_spin : ',frag%ke_spin + ! write(*,*) 'frag%ke_orbit : ',frag%ke_orbit + ! write(*,*) 'ke_remainder : ',frag%ke_budget - (frag%ke_orbit + frag%ke_spin) + lfailure = .false. + + end associate + return + + contains + function radial_objective_function(v_r_mag_input) result(fval) + !! Author: David A. Minton + !! + !! 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) :: v_r_mag_input !! Unknown radial component of fragment velocity vector + ! 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 + real(DP), dimension(:,:), allocatable :: v_shift + real(DP), dimension(frag%nbody) :: kearr + real(DP) :: keo, ke_radial + + associate(nfrag => frag%nbody) + allocate(v_shift, mold=frag%vb) + v_shift(:,:) = fraggle_util_vmag_to_vb(v_r_mag_input, frag%v_r_unit, frag%v_t_mag, frag%v_t_unit, frag%mass, frag%vbcom) + do concurrent(i = 1:nfrag) + kearr(i) = frag%mass(i) * (frag%Ip(3, i) * frag%radius(i)**2 * dot_product(frag%rot(:, i), frag%rot(:, i)) + dot_product(v_shift(:, i), v_shift(:, i))) + end do + keo = 2 * frag%ke_budget - sum(kearr(:)) + ke_radial = frag%ke_budget - frag%ke_orbit - frag%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 associate + + return + end function radial_objective_function + + end subroutine fraggle_generate_rad_vel + +end submodule s_fraggle_generate \ No newline at end of file diff --git a/src/fraggle/fraggle_placeholder.f90 b/src/fraggle/fraggle_placeholder.f90 new file mode 100644 index 000000000..bbf08bb04 --- /dev/null +++ b/src/fraggle/fraggle_placeholder.f90 @@ -0,0 +1,44 @@ +submodule(fraggle_classes) s_fraggle_placeholder + use swiftest + +contains + + !> The following interfaces are placeholders intended to satisfy the required abstract methods given by the parent class + module subroutine fraggle_placeholder_accel(self, system, param, t, lbeg) + implicit none + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current simulation time + logical, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step + write(*,*) "The type-bound procedure 'accel' is not defined for type fraggle_fragments" + return + end subroutine fraggle_placeholder_accel + + module subroutine fraggle_placeholder_kick(self, system, param, t, dt, lbeg) + implicit none + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system objec + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current time + real(DP), intent(in) :: dt !! Stepsize + logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. + + write(*,*) "The type-bound procedure 'kick' is not defined for type fraggle_fragments" + return + end subroutine fraggle_placeholder_kick + + module subroutine fraggle_placeholder_step(self, system, param, t, dt) + implicit none + class(fraggle_fragments), intent(inout) :: self !! Swiftest body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Simulation time + real(DP), intent(in) :: dt !! Current stepsize + + write(*,*) "The type-bound procedure 'step' is not defined for type fraggle_fragments" + return + end subroutine fraggle_placeholder_step + + +end submodule s_fraggle_placeholder \ No newline at end of file diff --git a/src/fraggle/fraggle_regime.f90 b/src/fraggle/fraggle_regime.f90 new file mode 100644 index 000000000..ab786917e --- /dev/null +++ b/src/fraggle/fraggle_regime.f90 @@ -0,0 +1,356 @@ +submodule(fraggle_classes) s_fraggle_regime + use swiftest + +contains + + module subroutine fraggle_regime_colliders(self, frag, system, param) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Determine which fragmentation regime the set of colliders will be. This subroutine is a wrapper for the non-polymorphic raggle_regime_collresolve subroutine. + !! It converts to SI units prior to calling + implicit none + ! Arguments + class(fraggle_colliders), intent(inout) :: self !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragment system object + class(swiftest_nbody_system), intent(in) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters + ! Internals + integer(I4B) :: jtarg, jproj, regime + real(DP), dimension(2) :: radius_si, mass_si, density_si + real(DP) :: min_mfrag_si, Mcb_si + real(DP), dimension(NDIM) :: x1_si, v1_si, x2_si, v2_si + real(DP) :: mlr, mslr, mtot, dentot, msys, msys_new, Qloss, impact_parameter + + associate(colliders => self) + ! Convert all quantities to SI units and determine which of the pair is the projectile vs. target before sending them to the regime determination subroutine + if (colliders%mass(1) > colliders%mass(2)) then + jtarg = 1 + jproj = 2 + else + jtarg = 2 + jproj = 1 + end if + mass_si(:) = colliders%mass([jtarg, jproj]) * param%MU2KG !! The two-body equivalent masses of the collider system + radius_si(:) = colliders%radius([jtarg, jproj]) * param%DU2M !! The two-body equivalent radii of the collider system + density_si(:) = mass_si(:) / (4.0_DP / 3._DP * PI * radius_si(:)**3) !! The two-body equivalent density of the collider system + x1_si(:) = colliders%xb(:,jtarg) * param%DU2M !! The first body of the two-body equivalent position vector the collider system + v1_si(:) = colliders%vb(:,jtarg) * param%DU2M / param%TU2S !! The first body of the two-body equivalent velocity vector the collider system + x2_si(:) = colliders%xb(:,jproj) * param%DU2M !! The second body of the two-body equivalent position vector the collider system + v2_si(:) = colliders%vb(:,jproj) * param%DU2M / param%TU2S !! The second body of the two-body equivalent velocity vector the collider system + Mcb_si = system%cb%mass * param%MU2KG !! The central body mass of the system + select type(param) + class is (symba_parameters) + min_mfrag_si = (param%min_GMfrag / param%GU) * param%MU2KG !! The minimum fragment mass to generate. Collider systems that would otherwise generate less massive fragments than this value will be forced to merge instead + class default + min_mfrag_si = 0.0_DP + end select + + mtot = sum(mass_si(:)) + dentot = sum(mass_si(:) * density_si(:)) / mtot + + !! Use the positions and velocities of the parents from indside the step (at collision) to calculate the collisional regime + call fraggle_regime_collresolve(Mcb_si, mass_si(jtarg), mass_si(jproj), radius_si(jtarg), radius_si(jproj), x1_si(:), x2_si(:),& + v1_si(:), v2_si(:), density_si(jtarg), density_si(jproj), min_mfrag_si, frag%regime, mlr, mslr, frag%Qloss) + + frag%mass_dist(1) = min(max(mlr, 0.0_DP), mtot) + frag%mass_dist(2) = min(max(mslr, 0.0_DP), mtot) + frag%mass_dist(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot) + + ! Convert quantities back to the system units and save them into the fragment system + frag%mass_dist(:) = (frag%mass_dist(:) / param%MU2KG) + frag%Qloss = frag%Qloss * (param%TU2S / param%DU2M)**2 / param%MU2KG + frag%mtot = sum(colliders%mass(:)) + frag%xbcom(:) = (colliders%mass(1) * colliders%xb(:,1) + colliders%mass(2) * colliders%xb(:,2)) / frag%mtot + frag%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / frag%mtot + end associate + + return + end subroutine fraggle_regime_colliders + + + subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, min_mfrag, regime, Mlr, Mslr, Qloss) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Determine the collisional regime of two colliding bodies. + !! Current version requires all values to be converted to SI units prior to calling the function + !! References: + !! Kokubo, E., Genda, H., 2010. Formation of Terrestrial Planets from Protoplanets Under a Realistic Accretion + !! Condition. ApJL 714, L21. https://doi.org/10.1088/2041-8205/714/1/L21 + !! Leinhardt, Z.M., Stewart, S.T., 2012. Collisions between Gravity-dominated Bodies. I. Outcome Regimes and Scaling + !! Laws 745, 79. https://doi.org/10.1088/0004-637X/745/1/79 + !! Mustill, A.J., Davies, M.B., Johansen, A., 2018. The dynamical evolution of transiting planetary systems including + !! a realistic collision prescription. Mon Not R Astron Soc 478, 2896–2908. https://doi.org/10.1093/mnras/sty1273 + !! Rufu, R., Aharonson, O., 2019. Impact Dynamics of Moons Within a Planetary Potential. J. Geophys. Res. Planets 124, + !! 1008–1019. https://doi.org/10.1029/2018JE005798 + !! Stewart, S.T., Leinhardt, Z.M., 2012. Collisions between Gravity-dominated Bodies. II. The Diversity of Impact + !! Outcomes during the End Stage of Planet Formation. ApJ 751, 32. https://doi.org/10.1088/0004-637X/751/1/32 + !! + implicit none + ! Arguments + real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, min_mfrag + real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 + integer(I4B), intent(out) :: regime + real(DP), intent(out) :: Mlr, Mslr + real(DP), intent(out) :: Qloss !! The residual energy after the collision + ! Constants + integer(I4B), parameter :: N1 = 1 !number of objects with mass equal to the largest remnant from LS12 + integer(I4B), parameter :: N2 = 2 !number of objects with mass larger than second largest remnant from LS12 + real(DP), parameter :: DENSITY1 = 1000.0_DP !standard density parameter from LS12 [kg/m3] + real(DP), parameter :: MU_BAR = 0.37_DP !0.385#0.37#0.3333# 3.978 # 1/3 material parameter for hydrodynamic planet-size bodies (LS12) + real(DP), parameter :: BETA = 2.85_DP !slope of sfd for remnants from LS12 2.85 + real(DP), parameter :: C1 = 2.43_DP !! Kokubo & Genda (2010) eq. (3) + real(DP), parameter :: C2 = -0.0408_DP !! Kokubo & Genda (2010) eq. (3) + real(DP), parameter :: C3 = 1.86_DP !! Kokubo & Genda (2010) eq. (3) + real(DP), parameter :: C4 = 1.08_DP !! Kokubo & Genda (2010) eq. (3) + real(DP), parameter :: CRUFU = 2.0_DP - 3 * MU_BAR ! central potential variable from Rufu and Aharonson (2019) + real(DP), parameter :: SUPERCAT_QRATIO = 1.8_DP ! See Section 4.1 of LS12 + ! Internals + real(DP) :: a1, alpha, aint, b, bcrit, c_star, egy, zeta, l, lint, mu, phi, theta + real(DP) :: Qr, Qrd_pstar, Qr_erosion, Qr_supercat + real(DP) :: Vhr, Verosion, Vescp, Vhill, Vimp, Vsupercat + real(DP) :: Mint, Mtot + real(DP) :: Rp, rhill + real(DP) :: Mresidual + real(DP) :: U_binding + + Vimp = norm2(vb2(:) - vb1(:)) + b = calc_b(xh2, vb2, xh1, vb1) + l = (rad1 + rad2) * (1 - b) + egy = 0.5_DP * dot_product(vb1, vb1) - GC * Mcb / norm2(xh1) + a1 = - GC * Mcb / 2.0_DP / egy + Mtot = m1 + m2 + mu = (m1 * m2) / Mtot + if (l < 2 * rad2) then + !calculate Mint + phi = 2 * acos((l - rad2) / rad2) + aint = rad2**2 * (PI - (phi - sin(phi)) / 2.0_DP) + lint = 2 * sqrt(rad2**2 - (rad2 - l / 2.0_DP) ** 2) + Mint = aint * lint ![kg] + alpha = (l**2) * (3 * rad2 - l) / (4 * (rad2**3)) + else + alpha = 1.0_DP + Mint = m2 + end if + Rp = (3 * (m1 / den1 + alpha * m2 / den2) / (4 * PI))**(1.0_DP/3.0_DP) ! (Mustill et al. 2018) + c_star = calc_c_star(Rp) + + !calculate Vescp + Vescp = sqrt(2 * GC * Mtot / Rp) !Mustill et al. 2018 eq 6 + + !calculate rhill + rhill = a1 * (m1 / 3.0_DP / (Mcb + m1))**(1.0_DP/3.0_DP) + + !calculate Vhill + if ((rad2 + rad1) < rhill) then + Vhill = sqrt(2 * GC * m1 * ((rhill**2 - rhill * (rad1 + rad2)) / & + (rhill**2 - 0.5_DP * (rad1 + rad2)**2)) / (rad1 + rad2)) + else + Vhill = Vescp + end if + + !calculate Qr_pstar + Qrd_pstar = calc_Qrd_pstar(m1, m2, alpha, c_star) * (Vhill / Vescp)**CRUFU !Rufu and Aharaonson eq (3) + + !calculate Verosion + Qr_erosion = 2 * (1.0_DP - m1 / Mtot) * Qrd_pstar + Verosion = (2 * Qr_erosion * Mtot / mu)** (1.0_DP / 2.0_DP) + Qr = mu*(Vimp**2) / Mtot / 2.0_DP + + !calculate mass largest remnant Mlr + Mlr = (1.0_DP - Qr / Qrd_pstar / 2.0_DP) * Mtot ! [kg] # LS12 eq (5) + + !calculate Vsupercat + Qr_supercat = SUPERCAT_QRATIO * Qrd_pstar ! See LS12 Section 4.1 + Vsupercat = sqrt(2 * Qr_supercat * Mtot / mu) + + !calculate Vhr + zeta = (m1 - m2) / Mtot + theta = 1.0_DP - b + Vhr = Vescp * (C1 * zeta**2 * theta**(2.5_DP) + C2 * zeta**2 + C3 * theta**(2.5_DP) + C4) ! Kokubo & Genda (2010) eq. (3) + bcrit = rad1 / (rad1 + rad2) + Qloss = 0.0_DP + U_binding = (3.0_DP * Mtot) / (5.0_DP * Rp) ! LS12 eq. 27 + + if ((m1 < min_mfrag).or.(m2 < min_mfrag)) then + regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime + Mlr = Mtot + Mslr = 0.0_DP + Qloss = 0.0_DP + write(*,*) "FORCE MERGE" + else + if( Vimp < Vescp) then + regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime + Mlr = Mtot + Mslr = 0.0_DP + Qloss = 0.0_DP + else if (Vimp < Verosion) then + if (b < bcrit) then + regime = COLLRESOLVE_REGIME_MERGE !partial accretion regime" + Mlr = Mtot + Mslr = 0.0_DP + Qloss = 0.0_DP + else if ((b > bcrit) .and. (Vimp < Vhr)) then + regime = COLLRESOLVE_REGIME_MERGE ! graze and merge + Mlr = Mtot + Mslr = 0.0_DP + Qloss = 0.0_DP + else + Mlr = m1 + Mslr = calc_Qrd_rev(m2, m1, Mint, den1, den2, Vimp, c_star) + regime = COLLRESOLVE_REGIME_HIT_AND_RUN !hit and run + Qloss = (c_star + 1.0_DP) * U_binding ! Qr + end if + else if (Vimp > Verosion .and. Vimp < Vsupercat) then + if (m2 < 0.001_DP * m1) then + regime = COLLRESOLVE_REGIME_MERGE !cratering regime" + Mlr = Mtot + Mslr = 0.0_DP + Qloss = 0.0_DP + else + Mslr = Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA) ! LS12 eq (37) + regime = COLLRESOLVE_REGIME_DISRUPTION !disruption + Qloss = (c_star + 1.0_DP) * U_binding ! Qr - Qr_erosion + end if + else if (Vimp > Vsupercat) then + Mlr = Mtot * 0.1_DP * (Qr / (Qrd_pstar * SUPERCAT_QRATIO))**(-1.5_DP) !LS12 eq (44) + Mslr = Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA) !LS12 eq (37) + regime = COLLRESOLVE_REGIME_SUPERCATASTROPHIC ! supercatastrophic + Qloss = (c_star + 1.0_DP) * U_binding ! Qr - Qr_supercat + else + write(*,*) "Error no regime found in symba_regime" + end if + end if + + Mresidual = Mtot - Mlr - Mslr + if (Mresidual < 0.0_DP) then ! prevents final masses from going negative + Mlr = Mlr + Mresidual + end if + + return + + ! Internal functions + contains + function calc_Qrd_pstar(Mtarg, Mp, alpha, c_star) result(Qrd_pstar) + !! author: Jennifer L.L. Pouplin and Carlisle A. Wishard + !! + !! Calculates the corrected Q* for oblique impacts. See Eq. (15) of LS12. + !! Reference: + !! Leinhardt, Z.M., Stewart, S.T., 2012. Collisions between Gravity-dominated Bodies. I. Outcome Regimes and Scaling + !! Laws 745, 79. https://doi.org/10.1088/0004-637X/745/1/79 + !! + implicit none + ! Arguments + real(DP),intent(in) :: Mtarg, Mp, alpha, c_star + ! Result + real(DP) :: Qrd_pstar + ! Internals + real(DP) :: Qrd_star1, mu_alpha, mu, Qrd_star + + ! calc mu, mu_alpha + mu = (Mtarg * Mp) / (Mtarg + Mp) ! [kg] + mu_alpha = (Mtarg * alpha * Mp) / (Mtarg + alpha * Mp) ! [kg] + ! calc Qrd_star1 + Qrd_star1 = (c_star * 4 * PI * DENSITY1 * GC * Rp**2) / 5.0_DP + ! calc Qrd_star + Qrd_star = Qrd_star1 * (((Mp / Mtarg + 1.0_DP)**2) / (4 * Mp / Mtarg))**(2.0_DP / (3.0_DP * MU_BAR) - 1.0_DP) !(eq 23) + ! calc Qrd_pstar, v_pstar + Qrd_pstar = ((mu / mu_alpha)**(2.0_DP - 3.0_DP * MU_BAR / 2.0_DP)) * Qrd_star ! (eq 15) + + return + end function calc_Qrd_pstar + + function calc_Qrd_rev(Mp, Mtarg, Mint, den1, den2, Vimp, c_star) result(Mslr) + !! author: Jennifer L.L. Pouplin and Carlisle A. Wishard + !! + !! Calculates mass of second largest fragment. + !! + implicit none + ! Arguments + real(DP),intent(in) :: Mp, Mtarg, Mint, den1, den2, Vimp, c_star + ! Result + real(DP) :: Mslr + ! Internals + real(DP) :: mtot_rev, mu_rev, gamma_rev, Qrd_star1, Qrd_star, mu_alpha_rev + real(DP) :: Qrd_pstar, Rc1, Qr_rev, Qrd_pstar_rev, Qr_supercat_rev + + ! calc Mslr, Rc1, mu, gammalr + mtot_rev = Mint + Mp + Rc1 = (3 * (Mint / den1 + Mp / den2) / (4 * PI))**(1.0_DP/3.0_DP) ! [m] Mustill et al 2018 + mu_rev = (Mint * Mp) / mtot_rev ! [kg] eq 49 LS12 + mu_alpha_rev = (Mtarg * alpha * Mp) / (Mtarg + alpha * Mp) + gamma_rev = Mint / Mp ! eq 50 LS12 + !calc Qr_rev + Qr_rev = mu_rev * (Vimp**2) / (2 * mtot_rev) + ! calc Qrd_star1, v_star1 + Qrd_star1 = (c_star * 4 * PI * mtot_rev * GC ) / Rc1 / 5.0_DP + ! calc Qrd_pstar_rev + Qrd_star = Qrd_star1 * (((gamma_rev + 1.0_DP)**2) / (4 * gamma_rev)) ** (2.0_DP / (3.0_DP * MU_BAR) - 1.0_DP) !(eq 52) + Qrd_pstar = Qrd_star * ((mu_rev / mu_alpha_rev)**(2.0_DP - 3.0_DP * MU_BAR / 2.0_DP)) + Qrd_pstar_rev = Qrd_pstar * (Vhill / Vescp)**CRUFU !Rufu and Aharaonson eq (3) + !calc Qr_supercat_rev + Qr_supercat_rev = 1.8_DP * Qrd_pstar_rev + if (Qr_rev > Qr_supercat_rev ) then + Mslr = mtot_rev * (0.1_DP * ((Qr_rev / (Qrd_pstar_rev * 1.8_DP))**(-1.5_DP))) !eq (44) + else if ( Qr_rev < Qrd_pstar_rev ) then + Mslr = Mp + else + Mslr = (1.0_DP - Qr_rev / Qrd_pstar_rev / 2.0_DP) * (mtot_rev) ! [kg] #(eq 5) + end if + + if ( Mslr > Mp ) Mslr = Mp !check conservation of mass + + return + end function calc_Qrd_rev + + function calc_b(proj_pos, proj_vel, targ_pos, targ_vel) result(sintheta) + !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Calculates the impact factor b = sin(theta), where theta is the angle between the relative velocity + !! and distance vectors of the target and projectile bodies. See Fig. 2 of Leinhardt and Stewart (2012) + !! + implicit none + !! Arguments + real(DP), dimension(:), intent(in) :: proj_pos, proj_vel, targ_pos, targ_vel + !! Result + real(DP) :: sintheta + !! Internals + real(DP), dimension(NDIM) :: imp_vel, distance, x_cross_v + + imp_vel(:) = proj_vel(:) - targ_vel(:) + distance(:) = proj_pos(:) - targ_pos(:) + x_cross_v(:) = distance(:) .cross. imp_vel(:) + sintheta = norm2(x_cross_v(:)) / norm2(distance(:)) / norm2(imp_vel(:)) + return + end function calc_b + + + function calc_c_star(Rc1) result(c_star) + !! author: David A. Minton + !! + !! Calculates c_star as a function of impact equivalent radius. It inteRpolates between 5 for ~1 km sized bodies to + !! 1.8 for ~10000 km sized bodies. See LS12 Fig. 4 for details. + !! + implicit none + !! Arguments + real(DP), intent(in) :: Rc1 + !! Result + real(DP) :: c_star + !! Internals + real(DP), parameter :: loR = 1.0e3_DP ! Lower bound of inteRpolation size (m) + real(DP), parameter :: hiR = 1.0e7_DP ! Upper bound of inteRpolation size (m) + real(DP), parameter :: loval = 5.0_DP ! Value of C* at lower bound + real(DP), parameter :: hival = 1.9_DP ! Value of C* at upper bound + + if (Rc1 < loR) then + c_star = loval + else if (Rc1 < hiR) then + c_star = loval + (hival - loval) * log(Rc1 / loR) / log(hiR /loR) + else + c_star = hival + end if + return + end function calc_c_star + + end subroutine fraggle_regime_collresolve + +end submodule s_fraggle_regime \ No newline at end of file diff --git a/src/fraggle/fraggle_set.f90 b/src/fraggle/fraggle_set.f90 new file mode 100644 index 000000000..b0b2a69be --- /dev/null +++ b/src/fraggle/fraggle_set.f90 @@ -0,0 +1,295 @@ +submodule(fraggle_classes) s_fraggle_set + use swiftest +contains + + module subroutine fraggle_set_budgets_fragments(self, colliders) + !! 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(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + ! Internals + real(DP) :: dEtot + real(DP), dimension(NDIM) :: dL + + associate(frag => self) + + dEtot = frag%Etot_after - frag%Etot_before + dL(:) = frag%Ltot_after(:) - frag%Ltot_before(:) + + frag%L_budget(:) = -dL(:) + frag%ke_budget = -(dEtot - 0.5_DP * frag%mtot * dot_product(frag%vbcom(:), frag%vbcom(:))) - frag%Qloss + + end associate + return + end subroutine fraggle_set_budgets_fragments + + + module subroutine fraggle_set_mass_dist_fragments(self, colliders) + !! author: David A. Minton + !! + !! Sets the mass of fragments based on the mass distribution returned by the regime calculation. + !! This subroutine must be run after the the setup rourtine has been run on the fragments + implicit none + ! Arguments + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + ! Internals + integer(I4B) :: i, istart, jproj, jtarg + real(DP), dimension(2) :: vol + real(DP) :: min_frag_mass + real(DP), dimension(NDIM) :: Ip_avg + + associate(frag => self, nfrag => self%nbody) + ! Get mass weighted mean of Ip and density + vol(1:2) = 4._DP / 3._DP * PI * colliders%radius(1:2)**3 + frag%density(1:nfrag) = frag%mtot / sum(vol(:)) + Ip_avg(:) = (colliders%mass(1) * colliders%Ip(:,1) + colliders%mass(2) * colliders%Ip(:,2)) / frag%mtot + do i = 1, nfrag + frag%Ip(:, i) = Ip_avg(:) + end do + + select case(frag%regime) + case(COLLRESOLVE_REGIME_DISRUPTION) + ! Distribute the mass among fragments, with a branch to check for the size of the second largest fragment + frag%mass(1) = frag%mass_dist(1) + if (frag%mass_dist(2) > frag%mass_dist(1) / 3._DP) then + frag%mass(2) = frag%mass_dist(2) + istart = 3 + else + istart = 2 + end if + ! Distribute remaining mass among the remaining bodies + do i = istart, nfrag + frag%mass(i) = (frag%mtot - sum(frag%mass(1:istart - 1))) / (nfrag - istart + 1) + end do + frag%radius(1:nfrag) = (3 * frag%mass(1:nfrag) / (4 * PI * frag%density(1:nfrag)))**(1.0_DP / 3.0_DP) + case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + ! If we are adding the first and largest fragment (lr), check to see if its mass is SMALLER than an equal distribution of + ! mass between all fragments. If so, we will just distribute the mass equally between the fragments + min_frag_mass = frag%mtot / nfrag + if (frag%mass_dist(1) < min_frag_mass) then + frag%mass(:) = min_frag_mass + else + frag%mass(1) = frag%mass_dist(1) + frag%mass(2:nfrag) = (frag%mtot - frag%mass_dist(1)) / (nfrag - 1) + end if + ! Distribute any residual mass if there is any and set the radius + frag%mass(nfrag) = frag%mass(nfrag) + (frag%mtot - sum(frag%mass(:))) + frag%radius(1:nfrag) = (3 * frag%mass(1:nfrag) / (4 * PI * frag%density(1:nfrag)))**(1.0_DP / 3.0_DP) + case(COLLRESOLVE_REGIME_HIT_AND_RUN) + if (colliders%mass(1) > colliders%mass(2)) then + jtarg = 1 + jproj = 2 + else + jtarg = 2 + jproj = 1 + end if + frag%mass(1) = frag%mass_dist(1) + frag%radius(1) = colliders%radius(jtarg) + frag%density(1) = frag%mass_dist(1) / vol(jtarg) + + frag%mass(2:nfrag) = (frag%mtot - frag%mass(1)) / (nfrag - 1) + frag%mass(nfrag) = frag%mass(nfrag) + (frag%mtot - sum(frag%mass(:))) + frag%radius(2:nfrag) = (3 * frag%mass(2:nfrag) / (4 * PI * frag%density(2:nfrag)))**(1.0_DP / 3.0_DP) + case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) + frag%mass(1) = frag%mtot + frag%radius(1) = (3 * sum(vol(:)) / (4 * PI))**(1._DP / 3._DP) + case default + write(*,*) "fraggle_set_mass_dist_fragments error: Unrecognized regime code",frag%regime + end select + + end associate + + return + end subroutine fraggle_set_mass_dist_fragments + + + module subroutine fraggle_set_coordinate_system(self, colliders) + !! author: David A. Minton + !! + !! Defines the collisional coordinate system, including the unit vectors of both the system and individual fragments. + implicit none + ! Arguments + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + ! Internals + integer(I4B) :: i + real(DP), dimension(NDIM) :: x_cross_v, delta_r, delta_v, Ltot + real(DP) :: r_col_norm, v_col_norm + real(DP), dimension(NDIM, self%nbody) :: L_sigma + + associate(frag => self, nfrag => self%nbody) + delta_v(:) = colliders%vb(:, 2) - colliders%vb(:, 1) + v_col_norm = .mag. delta_v(:) + delta_r(:) = colliders%xb(:, 2) - colliders%xb(:, 1) + r_col_norm = .mag. delta_r(:) + + ! 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. + Ltot = colliders%L_orbit(:,1) + colliders%L_orbit(:,2) + colliders%L_spin(:,1) + colliders%L_spin(:,2) + frag%y_coll_unit(:) = delta_r(:) / r_col_norm + frag%z_coll_unit(:) = Ltot(:) / (.mag. Ltot(:)) + ! The cross product of the y- by z-axis will give us the x-axis + frag%x_coll_unit(:) = frag%y_coll_unit(:) .cross. frag%z_coll_unit(:) + + if (.not.any(frag%x_coll(:,:) > 0.0_DP)) return + frag%rmag(:) = .mag. frag%x_coll(:,:) + + 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 + do concurrent(i = 1:nfrag, frag%rmag(i) > 0.0_DP) + frag%v_r_unit(:, i) = frag%x_coll(:, i) / frag%rmag(i) + frag%v_n_unit(:, i) = frag%z_coll_unit(:) + 2e-1_DP * (L_sigma(:,i) - 0.5_DP) + frag%v_n_unit(:, i) = frag%v_n_unit(:, i) / (.mag. frag%v_n_unit(:, i)) + frag%v_t_unit(:, i) = frag%v_n_unit(:, i) .cross. frag%v_r_unit(:, i) + frag%v_t_unit(:, i) = frag%v_t_unit(:, i) / (.mag. frag%v_t_unit(:, i)) + end do + end associate + + return + end subroutine fraggle_set_coordinate_system + + + ! module subroutine symba_set_collresolve_colliders(self, cb, pl, idx) + ! !! author: David A. Minton + ! !! + ! !! Calculate the two-body equivalent values given a set of input collider indices + ! use swiftest_classes, only : swiftest_nbody_system + ! implicit none + ! ! Arguments + ! class(fraggle_colliders), intent(inout) :: self !! Fraggle collider object + ! class(symba_cb), intent(in) :: cb !! Swiftest central body object system object + ! class(symba_pl), intent(in) :: pl !! Swiftest central body object system object + ! integer(I4B), dimension(:), intent(in) :: idx !! Index array of bodies from the pl object to use to calculate a "two-body equivalent" collisional pair + ! ! Internals + ! real(DP), dimension(NDIM, 2) :: mxc, vc + ! real(DP), dimension(NDIM) :: vcom, xcom + + ! associate(colliders => self) + + ! ! Compute orbital angular momentum of pre-impact system + ! xcom(:) = (colliders%mass(1) * colliders%xb(:, 1) + colliders%mass(2) * colliders%xb(:, 2)) / sum(colliders%mass(:)) + ! vcom(:) = (colliders%mass(1) * colliders%vb(:, 1) + colliders%mass(2) * colliders%vb(:, 2)) / sum(colliders%mass(:)) + ! mxc(:, 1) = colliders%mass(1) * (colliders%xb(:, 1) - xcom(:)) + ! mxc(:, 2) = colliders%mass(2) * (colliders%xb(:, 2) - xcom(:)) + ! vc(:, 1) = colliders%vb(:, 1) - vcom(:) + ! vc(:, 2) = colliders%vb(:, 2) - vcom(:) + + ! colliders%L_orbit(:,:) = mxc(:,:) .cross. vc(:,:) + + ! end associate + + ! return + ! end subroutine symbe_set_collresolve_colliders + + + module subroutine fraggle_set_natural_scale_factors(self, colliders) + !! 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(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + ! Internals + integer(I4B) :: i + + associate(frag => self) + ! Find the center of mass of the collisional system + frag%mtot = sum(colliders%mass(:)) + frag%xbcom(:) = (colliders%mass(1) * colliders%xb(:,1) + colliders%mass(2) * colliders%xb(:,2)) / frag%mtot + frag%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / frag%mtot + + ! Set scale factors + frag%Escale = 0.5_DP * (colliders%mass(1) * dot_product(colliders%vb(:,1), colliders%vb(:,1)) + colliders%mass(2) * dot_product(colliders%vb(:,2), colliders%vb(:,2))) + frag%dscale = sum(colliders%radius(:)) + frag%mscale = frag%mtot + frag%vscale = sqrt(frag%Escale / frag%mscale) + frag%tscale = frag%dscale / frag%vscale + frag%Lscale = frag%mscale * frag%dscale * frag%vscale + + ! Scale all dimensioned quantities of colliders and fragments + frag%xbcom(:) = frag%xbcom(:) / frag%dscale + frag%vbcom(:) = frag%vbcom(:) / frag%vscale + colliders%xb(:,:) = colliders%xb(:,:) / frag%dscale + colliders%vb(:,:) = colliders%vb(:,:) / frag%vscale + colliders%mass(:) = colliders%mass(:) / frag%mscale + colliders%radius(:) = colliders%radius(:) / frag%dscale + colliders%L_spin(:,:) = colliders%L_spin(:,:) / frag%Lscale + + do i = 1, 2 + colliders%rot(:,i) = colliders%L_spin(:,i) / (colliders%mass(i) * colliders%radius(i)**2 * colliders%Ip(3, i)) + end do + + frag%mtot = frag%mtot / frag%mscale + frag%mass = frag%mass / frag%mscale + frag%radius = frag%radius / frag%dscale + frag%Qloss = frag%Qloss / frag%Escale + end associate + + return + end subroutine fraggle_set_natural_scale_factors + + + module subroutine fraggle_set_original_scale_factors(self, colliders) + !! author: David A. Minton + !! + !! Restores dimenional quantities back to the system units + use, intrinsic :: ieee_exceptions + implicit none + ! Arguments + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + ! Internals + integer(I4B) :: i + logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes + + 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(frag => self) + + ! Restore scale factors + frag%xbcom(:) = frag%xbcom(:) * frag%dscale + frag%vbcom(:) = frag%vbcom(:) * frag%vscale + + colliders%mass = colliders%mass * frag%mscale + colliders%radius = colliders%radius * frag%dscale + colliders%xb = colliders%xb * frag%dscale + colliders%vb = colliders%vb * frag%vscale + colliders%L_spin = colliders%L_spin * frag%Lscale + do i = 1, 2 + colliders%rot(:,i) = colliders%L_spin(:,i) * (colliders%mass(i) * colliders%radius(i)**2 * colliders%Ip(3, i)) + end do + frag%Qloss = frag%Qloss * frag%Escale + + frag%mtot = frag%mtot * frag%mscale + frag%mass = frag%mass * frag%mscale + frag%radius = frag%radius * frag%dscale + frag%rot = frag%rot / frag%tscale + frag%x_coll = frag%x_coll * frag%dscale + frag%v_coll = frag%v_coll * frag%vscale + + do i = 1, frag%nbody + frag%xb(:, i) = frag%x_coll(:, i) + frag%xbcom(:) + frag%vb(:, i) = frag%v_coll(:, i) + frag%vbcom(:) + end do + + frag%mscale = 1.0_DP + frag%dscale = 1.0_DP + frag%vscale = 1.0_DP + frag%tscale = 1.0_DP + frag%Lscale = 1.0_DP + frag%Escale = 1.0_DP + end associate + call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) + + return + end subroutine fraggle_set_original_scale_factors + + +end submodule s_fraggle_set \ No newline at end of file diff --git a/src/fraggle/fraggle_setup.f90 b/src/fraggle/fraggle_setup.f90 new file mode 100644 index 000000000..65e107b03 --- /dev/null +++ b/src/fraggle/fraggle_setup.f90 @@ -0,0 +1,78 @@ +submodule (fraggle_classes) s_fraggle_setup + use swiftest +contains + + module subroutine fraggle_setup_reset_fragments(self) + !! author: David A. Minton + !! + !! 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) + implicit none + ! Arguments + class(fraggle_fragments), intent(inout) :: self + + self%xb(:,:) = 0.0_DP + self%vb(:,:) = 0.0_DP + self%rot(:,:) = 0.0_DP + self%x_coll(:,:) = 0.0_DP + self%v_coll(:,:) = 0.0_DP + self%v_r_unit(:,:) = 0.0_DP + self%v_t_unit(:,:) = 0.0_DP + self%v_n_unit(:,:) = 0.0_DP + + self%rmag(:) = 0.0_DP + self%rotmag(:) = 0.0_DP + self%v_r_mag(:) = 0.0_DP + self%v_t_mag(:) = 0.0_DP + + self%ke_orbit = 0.0_DP + self%ke_spin = 0.0_DP + self%L_orbit(:) = 0.0_DP + self%L_spin(:) = 0.0_DP + + return + end subroutine fraggle_setup_reset_fragments + + + module subroutine fraggle_setup_fragments(self, n, param) + !! author: David A. Minton + !! + !! Allocates arrays for n fragments in a Fraggle system. Passing n = 0 deallocates all arrays. + implicit none + ! Arguments + class(fraggle_fragments), intent(inout) :: self + integer(I4B), intent(in) :: n + class(swiftest_parameters), intent(in) :: param + + call setup_pl(self, n, param) + if (n < 0) return + + if (allocated(self%x_coll)) deallocate(self%x_coll) + if (allocated(self%v_coll)) deallocate(self%v_coll) + if (allocated(self%v_r_unit)) deallocate(self%v_r_unit) + if (allocated(self%v_t_unit)) deallocate(self%v_t_unit) + if (allocated(self%v_n_unit)) deallocate(self%v_n_unit) + if (allocated(self%rmag)) deallocate(self%rmag) + if (allocated(self%rotmag)) deallocate(self%rotmag) + if (allocated(self%v_r_mag)) deallocate(self%v_r_mag) + if (allocated(self%v_t_mag)) deallocate(self%v_t_mag) + + if (n == 0) return + + allocate(self%x_coll(NDIM,n)) + allocate(self%v_coll(NDIM,n)) + allocate(self%v_r_unit(NDIM,n)) + allocate(self%v_t_unit(NDIM,n)) + allocate(self%v_n_unit(NDIM,n)) + allocate(self%rmag(n)) + allocate(self%rotmag(n)) + allocate(self%v_r_mag(n)) + allocate(self%v_t_mag(n)) + + call self%reset() + + return + end subroutine fraggle_setup_fragments + + + +end submodule s_fraggle_setup \ No newline at end of file diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 new file mode 100644 index 000000000..df8b89c41 --- /dev/null +++ b/src/fraggle/fraggle_util.f90 @@ -0,0 +1,316 @@ +submodule(fraggle_classes) s_fraggle_util + use swiftest +contains + + module subroutine fraggle_util_add_fragments_to_system(frag, colliders, system, param) + !! Author: David A. Minton + !! + !! Adds fragments to the temporary system pl object + implicit none + ! Arguments + class(fraggle_fragments), intent(in) :: frag !! Fraggle fragment system object + class(fraggle_colliders), intent(in) :: colliders !! Fraggle collider system object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + ! Internals + integer(I4B) :: i, npl_before, npl_after + logical, dimension(:), allocatable :: lexclude + + associate(nfrag => frag%nbody, pl => system%pl, cb => system%cb) + npl_after = pl%nbody + npl_before = npl_after - nfrag + allocate(lexclude(npl_after)) + + pl%status(npl_before+1:npl_after) = ACTIVE + pl%mass(npl_before+1:npl_after) = frag%mass(1:nfrag) + pl%Gmass(npl_before+1:npl_after) = frag%mass(1:nfrag) * param%GU + pl%radius(npl_before+1:npl_after) = frag%radius(1:nfrag) + do concurrent (i = 1:nfrag) + pl%xb(:,npl_before+i) = frag%xb(:,i) + pl%vb(:,npl_before+i) = frag%vb(:,i) + pl%xh(:,npl_before+i) = frag%xb(:,i) - cb%xb(:) + pl%vh(:,npl_before+i) = frag%vb(:,i) - cb%vb(:) + end do + if (param%lrotation) then + pl%Ip(:,npl_before+1:npl_after) = frag%Ip(:,1:nfrag) + pl%rot(:,npl_before+1:npl_after) = frag%rot(:,1:nfrag) + end if + ! This will remove the colliders from the system since we've replaced them with fragments + lexclude(1:npl_after) = .false. + lexclude(colliders%idx(1:colliders%ncoll)) = .true. + where(lexclude(1:npl_after)) + pl%status(1:npl_after) = INACTIVE + elsewhere + pl%status(1:npl_after) = ACTIVE + endwhere + + end associate + + return + end subroutine fraggle_util_add_fragments_to_system + + + module subroutine fraggle_util_ang_mtm(self) + !! Author: David A. Minton + !! + !! Calcualtes the current angular momentum of the fragments + implicit none + ! Arguments + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + ! Internals + integer(I4B) :: i + + associate(frag => self, nfrag => self%nbody) + frag%L_orbit(:) = 0.0_DP + frag%L_spin(:) = 0.0_DP + + do i = 1, nfrag + frag%L_orbit(:) = frag%L_orbit(:) + frag%mass(i) * (frag%x_coll(:, i) .cross. frag%v_coll(:, i)) + frag%L_spin(:) = frag%L_spin(:) + frag%mass(i) * frag%radius(i)**2 * frag%Ip(:, i) * frag%rot(:, i) + end do + end associate + + return + end subroutine fraggle_util_ang_mtm + + + module subroutine fraggle_util_construct_temporary_system(frag, 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 + !! and optionally including fragments. + implicit none + ! Arguments + class(fraggle_fragments), intent(in) :: frag !! 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 + ! Internals + logical, dimension(:), allocatable :: linclude + integer(I4B) :: npl_tot + + associate(nfrag => frag%nbody, pl => system%pl, npl => system%pl%nbody, cb => system%cb) + ! Set up a new system based on the original + if (allocated(tmpparam)) deallocate(tmpparam) + if (allocated(tmpsys)) deallocate(tmpsys) + allocate(tmpparam, source=param) + call setup_construct_system(tmpsys, tmpparam) + + ! No test particles necessary for energy/momentum calcs + call tmpsys%tp%setup(0, param) + + ! Replace the empty central body object with a copy of the original + deallocate(tmpsys%cb) + allocate(tmpsys%cb, source=cb) + + ! Make space for the fragments + npl_tot = npl + nfrag + call tmpsys%pl%setup(npl_tot, tmpparam) + allocate(linclude(npl_tot)) + + ! Fill up the temporary system with all of the original bodies, leaving the spaces for fragments empty until we add them in later + linclude(1:npl) = .true. + linclude(npl+1:npl_tot) = .false. + call tmpsys%pl%fill(pl, linclude) + + ! Scale the temporary system to the natural units of the current Fraggle calculation + call tmpsys%rescale(tmpparam, frag%mscale, frag%dscale, frag%tscale) + + end associate + + return + end subroutine fraggle_util_construct_temporary_system + + + module subroutine fraggle_util_get_energy_momentum(self, colliders, system, param, lbefore) + !! Author: David A. Minton + !! + !! Calculates total system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) + !! This subrourtine works by building a temporary internal massive body object 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. + !! This will temporarily expand the massive body object in a temporary system object called tmpsys to feed it into symba_energy + implicit none + ! Arguments + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the system, with colliders included and fragments excluded or vice versa + ! Internals + integer(I4B) :: i, nplm + logical, dimension(:), allocatable :: lexclude + logical :: lk_plpl + logical, save :: ladd_frags + class(swiftest_nbody_system), allocatable, save :: tmpsys + class(swiftest_parameters), allocatable, save :: tmpparam + integer(I4B) :: npl_before, npl_after + + associate(frag => self, nfrag => self%nbody, pl => system%pl, cb => system%cb) + + ! Because we're making a copy of the massive body object 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. + + npl_before = pl%nbody + npl_after = npl_before + nfrag + + ! Build the exluded body logical mask + allocate(lexclude(npl_after)) + if (lbefore) then + lexclude(1:npl_before) = .false. + lexclude(npl_before+1:npl_after) = .true. + call fraggle_util_construct_temporary_system(frag, system, param, tmpsys, tmpparam) + else + lexclude(1:npl_after) = .false. + lexclude(colliders%idx(1:colliders%ncoll)) = .true. + if (.not.allocated(tmpsys)) then + write(*,*) "Error in fraggle_util_get_energy_momentum. This must be called with lbefore=.true. at least once before calling it with lbefore=.false." + call util_exit(FAILURE) + end if + call fraggle_util_add_fragments_to_system(frag, colliders, tmpsys, tmpparam) + end if + + call tmpsys%pl%index(param) + + call tmpsys%get_energy_and_momentum(param) + + + ! Calculate the current fragment energy and momentum balances + if (lbefore) then + frag%Lorbit_before(:) = tmpsys%Lorbit(:) + frag%Lspin_before(:) = tmpsys%Lspin(:) + frag%Ltot_before(:) = tmpsys%Ltot(:) + frag%ke_orbit_before = tmpsys%ke_orbit + frag%ke_spin_before = tmpsys%ke_spin + frag%pe_before = tmpsys%pe + frag%Etot_before = tmpsys%te + else + frag%Lorbit_after(:) = tmpsys%Lorbit(:) + frag%Lspin_after(:) = tmpsys%Lspin(:) + frag%Ltot_after(:) = tmpsys%Ltot(:) + frag%ke_orbit_after = tmpsys%ke_orbit + frag%ke_spin_after = tmpsys%ke_spin + frag%pe_after = tmpsys%pe + frag%Etot_after = tmpsys%te + end if + end associate + + return + end subroutine fraggle_util_get_energy_momentum + + + module subroutine fraggle_util_restructure(self, colliders, try, f_spin, r_max_start) + !! Author: David A. Minton + !! + !! Restructure the inputs after a failed attempt failed to find a set of positions and velocities that satisfy the energy and momentum constraints + implicit none + ! Arguments + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_colliders), intent(in) :: colliders !! Fraggle collider system object + integer(I4B), intent(in) :: try !! The current number of times Fraggle has tried to find a solution + real(DP), intent(inout) :: f_spin !! Fraction of energy/momentum that goes into spin. This decreases ater a failed attempt + real(DP), intent(inout) :: r_max_start !! The maximum radial distance that the position calculation starts with. This increases after a failed attempt + ! Internals + integer(I4B) :: i + real(DP), save :: ke_tot_deficit, r_max_start_old, ke_avg_deficit_old + real(DP), dimension(:), allocatable :: m_frag_new, rad_frag_new + real(DP), dimension(:,:), allocatable :: xb_frag_new, vb_frag_new, Ip_frag_new, rot_frag_new + real(DP) :: delta_r, delta_r_max, ke_avg_deficit + real(DP), parameter :: ke_avg_deficit_target = 0.0_DP + + ! Introduce a bit of noise in the radius determination so we don't just flip flop between similar failed positions + associate(frag => self) + call random_number(delta_r_max) + delta_r_max = sum(colliders%radius(:)) * (1.0_DP + 2e-1_DP * (delta_r_max - 0.5_DP)) + if (try == 1) then + ke_tot_deficit = - (frag%ke_budget - frag%ke_orbit - frag%ke_spin) + ke_avg_deficit = ke_tot_deficit + delta_r = delta_r_max + else + ! Linearly interpolate the last two failed solution ke deficits to find a new distance value to try + ke_tot_deficit = ke_tot_deficit - (frag%ke_budget - frag%ke_orbit - frag%ke_spin) + ke_avg_deficit = ke_tot_deficit / try + delta_r = (r_max_start - r_max_start_old) * (ke_avg_deficit_target - ke_avg_deficit_old) / (ke_avg_deficit - ke_avg_deficit_old) + if (abs(delta_r) > delta_r_max) delta_r = sign(delta_r_max, delta_r) + end if + r_max_start_old = r_max_start + r_max_start = r_max_start + delta_r ! The larger lever arm can help if the problem is in the angular momentum step + ke_avg_deficit_old = ke_avg_deficit + + if (f_spin > epsilon(1.0_DP)) then ! Try reducing the fraction in spin + f_spin = f_spin / 2 + else + f_spin = 0.0_DP + end if + end associate + + return + end subroutine fraggle_util_restructure + + + module subroutine fraggle_util_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 + + ! Internals + real(DP), dimension(NDIM) :: mvec_frag, COM_offset + integer(I4B) :: i, nfrag + real(DP) :: mtot + + mvec_frag(:) = 0.0_DP + mtot = sum(m_frag) + nfrag = size(m_frag) + + 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 + + return + end subroutine fraggle_util_shift_vector_to_origin + + + module function fraggle_util_vmag_to_vb(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom) result(vb) + !! Author: David A. Minton + !! + !! Converts radial and tangential velocity magnitudes into barycentric velocity + implicit none + ! Arguments + 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) :: m_frag !! Fragment masses + real(DP), dimension(:), intent(in) :: vcom !! Barycentric velocity of collisional system center of mass + ! Result + real(DP), dimension(:,:), allocatable :: vb + ! Internals + integer(I4B) :: i, nfrag + + allocate(vb, mold=v_r_unit) + ! Make sure the velocity magnitude stays positive + nfrag = size(m_frag) + do i = 1, nfrag + vb(:,i) = abs(v_r_mag(i)) * v_r_unit(:, i) + end do + ! 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 + call fraggle_util_shift_vector_to_origin(m_frag, vb) + + do i = 1, nfrag + vb(:, i) = vb(:, i) + v_t_mag(i) * v_t_unit(:, i) + vcom(:) + end do + + return + end function fraggle_util_vmag_to_vb + + +end submodule s_fraggle_util diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 deleted file mode 100644 index c62d337ca..000000000 --- a/src/fragmentation/fragmentation.f90 +++ /dev/null @@ -1,1266 +0,0 @@ -submodule(swiftest_classes) s_fragmentation - use swiftest -contains - - module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & - nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) - !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Initialize the position and velocity of fragments to conserve energy and momentum. - use, intrinsic :: ieee_exceptions - implicit none - ! Arguments - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - integer(I4B), dimension(:), intent(in) :: family !! Index of bodies involved in the collision - real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Two-body equivalent position, vector, spin momentum, and rotational inertia values for the collision - real(DP), dimension(:), intent(inout) :: mass, radius !! Two-body equivalent mass and radii for the bodies in the collision - integer(I4B), intent(inout) :: nfrag !! Number of fragments to generate - real(DP), dimension(:), allocatable, intent(inout) :: m_frag, rad_frag !! Distribution of fragment mass and radii - real(DP), dimension(:,:), allocatable, intent(inout) :: Ip_frag !! Fragment rotational inertia vectors - real(DP), dimension(:,:), allocatable, intent(inout) :: xb_frag, vb_frag, rot_frag !! Fragment barycentric position, barycentric velocity, and rotation vectors - real(DP), intent(inout) :: Qloss !! Energy lost during the collision - logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? - ! Internals - real(DP) :: mscale, dscale, vscale, tscale, Lscale, Escale ! Scale factors that reduce quantities to O(~1) in the collisional system - real(DP) :: mtot - real(DP), dimension(NDIM) :: xcom, vcom - integer(I4B) :: ii, npl_new - logical, dimension(:), allocatable :: lexclude - real(DP), dimension(NDIM, 2) :: rot, L_orb, mxc, vc - real(DP), dimension(:,:), allocatable :: x_frag, v_frag, v_r_unit, v_t_unit, v_h_unit - real(DP), dimension(:), allocatable :: rmag, rotmag, v_r_mag, v_t_mag - real(DP), dimension(NDIM) :: Ltot_before - real(DP), dimension(NDIM) :: Ltot_after - real(DP) :: Etot_before, ke_orbit_before, ke_spin_before, pe_before, Lmag_before - real(DP) :: Etot_after, ke_orbit_after, ke_spin_after, pe_after, Lmag_after, dEtot, dLmag - real(DP), dimension(NDIM) :: L_frag_tot, L_frag_orb, L_frag_spin, L_frag_budget, Lorbit_before, Lorbit_after, Lspin_before, Lspin_after, dL - real(DP) :: ke_frag_budget, ke_frag_orbit, ke_frag_spin, ke_tot_deficit, ke_avg_deficit, ke_avg_deficit_old - real(DP), dimension(NDIM) :: x_col_unit, y_col_unit, z_col_unit - character(len=*), parameter :: fmtlabel = "(A14,10(ES11.4,1X,:))" - integer(I4B) :: try - integer(I4B), parameter :: NFRAG_MIN = 7 !! The minimum allowable number of fragments (set to 6 because that's how many unknowns are needed in the tangential velocity calculation) - real(DP) :: r_max_start, r_max_start_old, r_max, f_spin - real(DP), parameter :: Ltol = 10 * epsilon(1.0_DP) - real(DP), parameter :: Etol = 1e-8_DP - integer(I4B), parameter :: MAXTRY = 3000 - logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes - - if (nfrag < NFRAG_MIN) then - write(*,*) "symba_frag_pos needs at least ",NFRAG_MIN," fragments, but only ",nfrag," were given." - lfailure = .true. - return - end if - - call ieee_get_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily - fpe_quiet_modes(:) = .false. - call ieee_set_halting_mode(IEEE_ALL,fpe_quiet_modes) - - allocate(x_frag, source=xb_frag) - allocate(v_frag, source=vb_frag) - allocate(rmag(nfrag)) - allocate(rotmag(nfrag)) - allocate(v_r_mag(nfrag)) - allocate(v_t_mag(nfrag)) - allocate(v_r_unit(NDIM,nfrag)) - allocate(v_t_unit(NDIM,nfrag)) - allocate(v_h_unit(NDIM,nfrag)) - - rmag(:) = 0.0_DP - rotmag(:) = 0.0_DP - v_r_mag(:) = 0.0_DP - v_t_mag(:) = 0.0_DP - v_r_unit(:,:) = 0.0_DP - v_t_unit(:,:) = 0.0_DP - v_h_unit(:,:) = 0.0_DP - - associate(pl => system%pl, npl => system%pl%nbody) - npl_new = npl + nfrag - allocate(lexclude(npl_new)) - lexclude(1:npl) = pl%status(1:npl) == INACTIVE - lexclude(npl+1:npl_new) = .true. - end associate - - call set_scale_factors() - - ! Compute orbital angular momentum of pre-impact system - mxc(:, 1) = mass(1) * (x(:, 1) - xcom(:)) - mxc(:, 2) = mass(2) * (x(:, 2) - xcom(:)) - vc(:, 1) = v(:, 1) - vcom(:) - vc(:, 2) = v(:, 2) - vcom(:) - L_orb(:,:) = mxc(:,:) .cross. vc(:,:) - - ! Compute orbital angular momentum of pre-impact system. We'll use this to start the coordinate system, but it will get updated as we divide up the angular momentum - L_frag_orb(:) = L_orb(:, 1) + L_orb(:, 2) - L_frag_spin(:) = L_spin(:, 1) + L_spin(:, 2) - L_frag_budget(:) = L_frag_orb(:) + L_frag_spin(:) - f_spin = 0.05_DP - - call reset_fragments() - call define_coordinate_system() - - ! Calculate the initial energy of the system without the collisional family - call calculate_system_energy(linclude_fragments=.false.) - - r_max_start = 1 * norm2(x(:,2) - x(:,1)) - try = 1 - lfailure = .false. - ke_tot_deficit = 0.0_DP - do while (try < MAXTRY) - lfailure = .false. - - call set_fragment_position_vectors() - - do concurrent (ii = 1:nfrag) - vb_frag(:, ii) = vcom(:) - end do - - call calculate_system_energy(linclude_fragments=.true.) - L_frag_budget(:) = -dL(:) - ! The ke constraints are calcualted in the collision frame, so undo the barycentric velocity component - ke_frag_budget = -(dEtot - 0.5_DP * mtot * dot_product(vcom(:), vcom(:))) - Qloss - - call set_fragment_spin(lfailure) - if (.not.lfailure) call set_fragment_tan_vel(lfailure) - - if (lfailure) then - ! write(*,*) 'Failed to find tangential velocities' - else - call set_fragment_radial_velocities(lfailure) - ! if (lfailure) write(*,*) 'Failed to find radial velocities' - if (.not.lfailure) then - call calculate_system_energy(linclude_fragments=.true.) - if ((abs(dEtot + Qloss) > Etol) .or. (dEtot > 0.0_DP)) then - ! write(*,*) 'Failed due to high energy error: ',dEtot, abs(dEtot + Qloss) / Etol - lfailure = .true. - else if (abs(dLmag) / Lmag_before > Ltol) then - ! write(*,*) 'Failed due to high angular momentum error: ', dLmag / Lmag_before - lfailure = .true. - end if - end if - end if - - if (.not.lfailure) exit - call restructure_failed_fragments() - call reset_fragments() - try = try + 1 - end do - call restore_scale_factors() - - ! write(*, "(' -------------------------------------------------------------------------------------')") - ! write(*, "(' Final diagnostic')") - ! write(*, "(' -------------------------------------------------------------------------------------')") - ! call calculate_system_energy(linclude_fragments=.true.) - if (lfailure) then - write(*,*) "symba_frag_pos failed after: ",try," tries" - do ii = 1, nfrag - vb_frag(:, ii) = vcom(:) - end do - else - write(*,*) "symba_frag_pos succeeded after: ",try," tries" - ! write(*, "(' dL_tot should be very small' )") - ! write(*,fmtlabel) ' dL_tot |', dLmag / Lmag_before - ! write(*, "(' dE_tot should be negative and equal to Qloss' )") - ! write(*,fmtlabel) ' dE_tot |', dEtot / abs(Etot_before) - ! write(*,fmtlabel) ' Qloss |', -Qloss / abs(Etot_before) - ! write(*,fmtlabel) ' dE - Qloss |', (Etot_after - Etot_before + Qloss) / abs(Etot_before) - end if - ! write(*, "(' -------------------------------------------------------------------------------------')") - - call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily - - 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 - integer(I4B) :: i - - ! Find the center of mass of the collisional system - mtot = sum(mass(:)) - xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot - vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot - - ! Set scale factors - Escale = 0.5_DP * (mass(1) * dot_product(v(:,1), v(:,1)) + mass(2) * dot_product(v(:,2), v(:,2))) - dscale = sum(radius(:)) - mscale = mtot - vscale = sqrt(Escale / mscale) - tscale = dscale / vscale - Lscale = mscale * dscale * vscale - - xcom(:) = xcom(:) / dscale - vcom(:) = vcom(:) / vscale - - mtot = mtot / mscale - mass = mass / mscale - radius = radius / dscale - x = x / dscale - v = v / vscale - L_spin = L_spin / Lscale - do i = 1, 2 - rot(:,i) = L_spin(:,i) / (mass(i) * radius(i)**2 * Ip(3, i)) - end do - - m_frag = m_frag / mscale - rad_frag = rad_frag / dscale - Qloss = Qloss / Escale - - return - end subroutine set_scale_factors - - - subroutine restore_scale_factors() - !! author: David A. Minton - !! - !! Restores dimenional quantities back to the system units - implicit none - integer(I4B) :: i - - call ieee_set_halting_mode(IEEE_ALL,.false.) - ! Restore scale factors - xcom(:) = xcom(:) * dscale - vcom(:) = vcom(:) * vscale - - mtot = mtot * mscale - mass = mass * mscale - radius = radius * dscale - x = x * dscale - v = v * vscale - L_spin = L_spin * Lscale - do i = 1, 2 - rot(:,i) = L_spin(:,i) * (mass(i) * radius(i)**2 * Ip(3, i)) - end do - - m_frag = m_frag * mscale - rad_frag = rad_frag * dscale - rot_frag = rot_frag / tscale - x_frag = x_frag * dscale - v_frag = v_frag * vscale - Qloss = Qloss * Escale - - do i = 1, nfrag - xb_frag(:, i) = x_frag(:, i) + xcom(:) - vb_frag(:, i) = v_frag(:, i) + vcom(:) - end do - - Etot_before = Etot_before * Escale - pe_before = pe_before * Escale - ke_spin_before = ke_spin_before * Escale - ke_orbit_before = ke_orbit_before * Escale - Ltot_before = Ltot_before * Lscale - Lmag_before = Lmag_before * Lscale - Etot_after = Etot_after * Escale - pe_after = pe_after * Escale - ke_spin_after = ke_spin_after * Escale - ke_orbit_after = ke_orbit_after * Escale - Ltot_after = Ltot_after * Lscale - Lmag_after = Lmag_after * Lscale - - dL(:) = Ltot_after(:) - Ltot_before(:) - dLmag = .mag.dL(:) - dEtot = Etot_after - Etot_before - - mscale = 1.0_DP - dscale = 1.0_DP - vscale = 1.0_DP - tscale = 1.0_DP - Lscale = 1.0_DP - Escale = 1.0_DP - - return - end subroutine restore_scale_factors - - subroutine reset_fragments() - !! author: David A. Minton - !! - !! Resets all tracked fragment quantities in order to do a fresh calculation - !! Initialize the fragments with 0 velocity and spin so we can divide up the balance between the tangential, radial, and spin components while conserving momentum - implicit none - - xb_frag(:,:) = 0.0_DP - vb_frag(:,:) = 0.0_DP - x_frag(:,:) = 0.0_DP - v_frag(:,:) = 0.0_DP - rot_frag(:,:) = 0.0_DP - v_t_mag(:) = 0.0_DP - v_r_mag(:) = 0.0_DP - ke_frag_orbit = 0.0_DP - ke_frag_spin = 0.0_DP - L_frag_orb(:) = 0.0_DP - L_frag_spin(:) = 0.0_DP - - return - end subroutine reset_fragments - - - 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, delta_r, delta_v - real(DP) :: r_col_norm, v_col_norm - real(DP), dimension(NDIM, nfrag) :: L_sigma - - delta_v(:) = v(:, 2) - v(:, 1) - v_col_norm = .mag. delta_v(:) - delta_r(:) = x(:, 2) - x(:, 1) - r_col_norm = .mag. delta_r(:) - - ! 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_budget(:) - L_frag_spin(:)) / (.mag. (L_frag_budget(:) - L_frag_spin(:))) - ! The cross product of the y- by z-axis will give us the x-axis - x_col_unit(:) = y_col_unit(:) .cross. z_col_unit(:) - - if (.not.any(x_frag(:,:) > 0.0_DP)) return - rmag(:) = .mag. x_frag(:,:) - - 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 - do concurrent(i = 1:nfrag, rmag(i) > 0.0_DP) - v_r_unit(:, i) = x_frag(:, i) / rmag(i) - v_h_unit(:, i) = z_col_unit(:) + 2e-1_DP * (L_sigma(:,i) - 0.5_DP) - v_h_unit(:, i) = v_h_unit(:, i) / (.mag. v_h_unit(:, i)) - v_t_unit(:, i) = v_h_unit(:, i) .cross. v_r_unit(:, i) - v_t_unit(:, i) = v_t_unit(:, i) / (.mag. v_t_unit(:, i)) - end do - - return - end subroutine define_coordinate_system - - - subroutine construct_temporary_system(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 - !! and optionally including fragments. - implicit none - ! Arguments - class(swiftest_nbody_system), allocatable, intent(inout) :: tmpsys - class(swiftest_parameters), allocatable, intent(inout) :: tmpparam - ! Internals - logical, dimension(:), allocatable :: lexclude_tmp - - associate(pl => system%pl, npl => system%pl%nbody, cb => system%cb) - if (size(lexclude) /= npl + nfrag) then - allocate(lexclude_tmp(npl + nfrag)) - lexclude_tmp(1:npl) = lexclude(1:npl) - call move_alloc(lexclude_tmp, lexclude) - end if - where (pl%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 - lexclude(npl+1:(npl + nfrag)) = .true. - allocate(tmpparam, source=param) - call setup_construct_system(tmpsys, tmpparam) - call tmpsys%tp%setup(0, param) - deallocate(tmpsys%cb) - allocate(tmpsys%cb, source=cb) - call tmpsys%pl%setup(npl + nfrag, tmpparam) - call tmpsys%pl%fill(pl, .not.lexclude) - call tmpsys%rescale(tmpparam, mscale, dscale, tscale) - - end associate - - return - end subroutine construct_temporary_system - - - subroutine add_fragments_to_tmpsys(tmpsys, tmpparam) - !! Author: David A. Minton - !! - !! Adds fragments to the temporary system pl object - implicit none - ! Arguments - class(swiftest_nbody_system), intent(inout) :: tmpsys - class(swiftest_parameters), intent(inout) :: tmpparam - ! Internals - integer(I4B) :: i - class(swiftest_pl), allocatable :: pl_discards - logical, dimension(:), allocatable :: lexclude_tmp - - associate(pl => system%pl, npl => system%pl%nbody) - npl_new = npl + nfrag - - tmpsys%pl%mass(npl+1:npl_new) = m_frag(1:nfrag) - tmpsys%pl%Gmass(npl+1:npl_new) = m_frag(1:nfrag) * tmpparam%GU - tmpsys%pl%radius(npl+1:npl_new) = rad_frag(1:nfrag) - do concurrent (i = 1:nfrag) - tmpsys%pl%xb(:,npl+i) = xb_frag(:,i) - tmpsys%pl%vb(:,npl+i) = vb_frag(:,i) - tmpsys%pl%xh(:,npl+i) = xb_frag(:,i) - tmpsys%cb%xb(:) - tmpsys%pl%vh(:,npl+i) = vb_frag(:,i) - tmpsys%cb%vb(:) - end do - if (tmpparam%lrotation) then - tmpsys%pl%Ip(:,npl+1:npl_new) = Ip_frag(:,1:nfrag) - tmpsys%pl%rot(:,npl+1:npl_new) = rot_frag(:,1:nfrag) - end if - ! Disable the collisional family for subsequent energy calculations and coordinate shifts - lexclude(family(:)) = .true. - lexclude(npl+1:npl_new) = .false. - where(lexclude(1:npl_new)) - tmpsys%pl%status(1:npl_new) = INACTIVE - elsewhere - tmpsys%pl%status(1:npl_new) = ACTIVE - end where - allocate(pl_discards, mold=tmpsys%pl) - call tmpsys%pl%spill(pl_discards, lspill_list=lexclude(1:npl_new), ldestructive=.true.) - npl_new = count(.not.lexclude(:)) - - if (size(lexclude) /= npl_new) then - allocate(lexclude_tmp(npl_new)) - call move_alloc(lexclude_tmp, lexclude) - end if - lexclude(1:npl_new) = .false. - - end associate - - return - end subroutine add_fragments_to_tmpsys - - - subroutine calculate_system_energy(linclude_fragments) - !! Author: David A. Minton - !! - !! Calculates total system energy, including all bodies in the pl list that do not have a corresponding value of the lexclude array that is true - !! and optionally including fragments. - implicit none - ! Arguments - logical, intent(in) :: linclude_fragments - ! Internals - integer(I4B) :: i, nplm - logical, dimension(:), allocatable :: lexclude_tmp - logical :: lk_plpl - class(swiftest_nbody_system), allocatable :: tmpsys - class(swiftest_parameters), allocatable :: tmpparam - - ! 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. This will temporarily expand the planet list in a temporary system object called tmpsys to feed it into symba_energy - associate(pl => system%pl, npl => system%pl%nbody, cb => system%cb) - - ! 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. - lk_plpl = allocated(pl%k_plpl) - if (lk_plpl) deallocate(pl%k_plpl) - - call construct_temporary_system(tmpsys, tmpparam) - if (linclude_fragments) call add_fragments_to_tmpsys(tmpsys, tmpparam) - - call tmpsys%pl%index(param) - - call tmpsys%get_energy_and_momentum(param) - - ! Restore the big array - deallocate(tmpsys%pl%k_plpl) - - if (lk_plpl) call pl%index(param) - - ! Calculate the current fragment energy and momentum balances - if (linclude_fragments) then - Lorbit_after(:) = tmpsys%Lorbit - Lspin_after(:) = tmpsys%Lspin - Ltot_after(:) = tmpsys%Lorbit(:) + tmpsys%Lspin(:) - Lmag_after = norm2(Ltot_after(:)) - ke_orbit_after = tmpsys%ke_orbit - ke_spin_after = tmpsys%ke_spin - pe_after = tmpsys%pe - Etot_after = tmpsys%te - dEtot = Etot_after - Etot_before - dL(:) = Ltot_after(:) - Ltot_before(:) - dLmag = .mag.dL(:) - else - Lorbit_before(:) = tmpsys%Lorbit - Lspin_before(:) = tmpsys%Lspin - Ltot_before(:) = tmpsys%Lorbit(:) + tmpsys%Lspin(:) - Lmag_before = norm2(Ltot_before(:)) - ke_orbit_before = tmpsys%ke_orbit - ke_spin_before = tmpsys%ke_spin - pe_before = tmpsys%pe - Etot_before = tmpsys%te - end if - end associate - - return - end subroutine calculate_system_energy - - - subroutine calculate_fragment_ang_mtm() - !! Author: David A. Minton - !! - !! Calcualtes the current angular momentum of the fragments - implicit none - integer(I4B) :: i - - L_frag_orb(:) = 0.0_DP - L_frag_spin(:) = 0.0_DP - - do i = 1, nfrag - L_frag_orb(:) = L_frag_orb(:) + m_frag(i) * (x_frag(:, i) .cross. v_frag(:, i)) - L_frag_spin(:) = L_frag_spin(:) + m_frag(i) * rad_frag(i)**2 * Ip_frag(:, i) * rot_frag(:, i) - end do - - L_frag_tot(:) = L_frag_orb(:) + L_frag_spin(:) - - return - end subroutine calculate_fragment_ang_mtm - - - 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 - - ! Internals - real(DP), dimension(NDIM) :: mvec_frag, COM_offset - integer(I4B) :: i - - mvec_frag(:) = 0.0_DP - - 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 - - return - end subroutine shift_vector_to_origin - - - 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 - real(DP) :: dis, rad - logical, dimension(:), allocatable :: loverlap - integer(I4B) :: i, j - - allocate(loverlap(nfrag)) - - ! Place the fragments into a region that is big enough that we should usually not have overlapping bodies - ! An overlapping bodies will collide in the next time step, so it's not a major problem if they do (it just slows the run down) - r_max = r_max_start - rad = sum(radius(:)) - - ! We will treat the first two fragments of the list as special cases. They get initialized the maximum distances apart along the original impactor distance vector. - ! This is done because in a regular disruption, the first body is the largest, the second the second largest, and the rest are smaller equal-mass fragments. - - call random_number(x_frag(:,3:nfrag)) - loverlap(:) = .true. - do while (any(loverlap(3:nfrag))) - x_frag(:, 1) = x(:, 1) - xcom(:) - x_frag(:, 2) = x(:, 2) - xcom(:) - r_max = r_max + 0.1_DP * rad - do i = 3, nfrag - if (loverlap(i)) then - call random_number(x_frag(:,i)) - x_frag(:, i) = 2 * (x_frag(:, i) - 0.5_DP) * r_max - end if - end do - loverlap(:) = .false. - do j = 1, nfrag - do i = j + 1, nfrag - dis = norm2(x_frag(:,j) - x_frag(:,i)) - loverlap(i) = loverlap(i) .or. (dis <= (rad_frag(i) + rad_frag(j))) - end do - end do - end do - call shift_vector_to_origin(m_frag, x_frag) - call define_coordinate_system() - - do i = 1, nfrag - xb_frag(:,i) = x_frag(:,i) + xcom(:) - end do - - xcom(:) = 0.0_DP - do i = 1, nfrag - xcom(:) = xcom(:) + m_frag(i) * xb_frag(:,i) - end do - xcom(:) = xcom(:) / mtot - - return - end subroutine set_fragment_position_vectors - - - subroutine set_fragment_spin(lerr) - !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Sets the spins of a collection of fragments such that they conserve angular momentum without blowing the fragment kinetic energy budget. - !! - !! A failure will trigger a restructuring of the fragments so we will try new values of the radial position distribution. - implicit none - ! Arguments - logical, intent(out) :: lerr - ! Internals - real(DP), dimension(NDIM) :: L_remainder, rot_L, rot_ke - integer(I4B) :: i - - lerr = .false. - - ! Start the first two bodies with the same rotation as the original two impactors, then distribute the remaining angular momentum among the rest - rot_frag(:,1:2) = rot(:, :) - rot_frag(:,3:nfrag) = 0.0_DP - call calculate_fragment_ang_mtm() - L_remainder(:) = L_frag_budget(:) - L_frag_spin(:) - - ke_frag_spin = 0.0_DP - do i = 1, nfrag - ! Convert a fraction (f_spin) of either the remaining angular momentum or kinetic energy budget into spin, whichever gives the smaller rotation so as not to blow any budgets - rot_ke(:) = sqrt(2 * f_spin * ke_frag_budget / (nfrag * m_frag(i) * rad_frag(i)**2 * Ip_frag(3, i))) * L_remainder(:) / norm2(L_remainder(:)) - rot_L(:) = f_spin * L_remainder(:) / (nfrag * m_frag(i) * rad_frag(i)**2 * Ip_frag(3, i)) - if (norm2(rot_ke) < norm2(rot_L)) then - rot_frag(:,i) = rot_frag(:, i) + rot_ke(:) - else - rot_frag(:, i) = rot_frag(:, i) + rot_L(:) - end if - ke_frag_spin = ke_frag_spin + m_frag(i) * Ip_frag(3, i) * rad_frag(i)**2 * dot_product(rot_frag(:, i), rot_frag(:, i)) - end do - ke_frag_spin = 0.5_DP * ke_frag_spin - - lerr = ((ke_frag_budget - ke_frag_spin - ke_frag_orbit) < 0.0_DP) - - return - end subroutine set_fragment_spin - - - subroutine set_fragment_tan_vel(lerr) - !! 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. - !! This procedure works in several stages, with a goal to solve the angular and linear momentum constraints on the fragments, while still leaving a positive balance of - !! our fragment kinetic energy (ke_frag_budget) that we can put into the radial velocity distribution. - !! - !! The first thing we'll try to do is solve for the tangential velocities of the first 6 fragments, using angular and linear momentum as constraints and an initial - !! tangential velocity distribution for the remaining bodies (if there are any) that distributes their angular momentum equally between them. - !! If that doesn't work and we blow our kinetic energy budget, we will attempt to find a tangential velocity distribution that minimizes the kinetic energy while - !! conserving momentum. - !! - !! A failure will trigger a restructuring of the fragments so we will try new values of the radial position distribution. - implicit none - ! Arguments - logical, intent(out) :: lerr - ! Internals - integer(I4B) :: i - real(DP), parameter :: TOL_MIN = 1e-1_DP ! This doesn't have to be very accurate, as we really just want a tangential velocity distribution with less kinetic energy than our initial guess. - real(DP), parameter :: TOL_INIT = 1e-14_DP - integer(I4B), parameter :: MAXLOOP = 10 - real(DP) :: tol - real(DP), dimension(:), allocatable :: v_t_initial - real(DP), dimension(nfrag) :: kefrag - type(lambda_obj) :: spinfunc - type(lambda_obj_err) :: objective_function - real(DP), dimension(NDIM) :: Li, L_remainder - - lerr = .false. - - ! write(*,*) '***************************************************' - ! write(*,*) 'Original dis : ',norm2(x(:,2) - x(:,1)) - ! write(*,*) 'r_max : ',r_max - ! write(*,*) 'f_spin : ',f_spin - ! write(*,*) '***************************************************' - ! write(*,*) 'Energy balance so far: ' - ! write(*,*) 'ke_frag_budget : ',ke_frag_budget - ! write(*,*) 'ke_orbit_before: ',ke_orbit_before - ! write(*,*) 'ke_orbit_after : ',ke_orbit_after - ! write(*,*) 'ke_spin_before : ',ke_spin_before - ! write(*,*) 'ke_spin_after : ',ke_spin_after - ! write(*,*) 'pe_before : ',pe_before - ! write(*,*) 'pe_after : ',pe_after - ! write(*,*) 'Qloss : ',Qloss - ! write(*,*) '***************************************************' - - allocate(v_t_initial, mold=v_t_mag) - v_t_initial(:) = 0.0_DP - v_frag(:,:) = 0.0_DP - - ! Next we will solve for the tangential component of the velocities that both conserves linear momentum and uses the remaining angular momentum not used in spin. - ! This will be done using a linear solver that solves for the tangential velocities of the first 6 fragments, constrained by the linear and angular momentum vectors, - ! which is embedded in a non-linear minimizer that will adjust the tangential velocities of the remaining i>6 fragments to minimize kinetic energy for a given momentum solution - ! The initial conditions fed to the minimizer for the fragments will be the remaining angular momentum distributed between the fragments. - call calculate_fragment_ang_mtm() - call define_coordinate_system() ! Make sure that the tangential velocity components are defined properly - L_remainder(:) = L_frag_budget(:) - L_frag_spin(:) - do i = 1, nfrag - v_t_initial(i) = norm2(L_remainder(:)) / ((nfrag - i + 1) * m_frag(i) * norm2(x_frag(:,i))) - Li(:) = m_frag(i) * (x_frag(:,i) .cross. (v_t_initial(i) * v_t_unit(:, i))) - L_remainder(:) = L_remainder(:) - Li(:) - end do - - ! Find the local kinetic energy minimum for the system that conserves linear and angular momentum - objective_function = lambda_obj(tangential_objective_function, lerr) - - tol = TOL_INIT - do while(tol < TOL_MIN) - v_t_mag(7:nfrag) = util_minimize_bfgs(objective_function, nfrag-6, v_t_initial(7:nfrag), tol, MAXLOOP, lerr) - ! Now that the KE-minimized values of the i>6 fragments are found, calculate the momentum-conserving solution for tangential velociteis - v_t_initial(7:nfrag) = v_t_mag(7:nfrag) - if (.not.lerr) exit - tol = tol * 2_DP ! Keep increasing the tolerance until we converge on a solution - end do - v_t_mag(1:nfrag) = solve_fragment_tan_vel(v_t_mag_input=v_t_initial(7:nfrag), lerr=lerr) - - ! Perform one final shift of the radial velocity vectors to align with the center of mass of the collisional system (the origin) - vb_frag(:,1:nfrag) = vmag_to_vb(v_r_mag(1:nfrag), v_r_unit(:,1:nfrag), v_t_mag(1:nfrag), v_t_unit(:,1:nfrag), m_frag(1:nfrag), vcom(:)) - do concurrent (i = 1:nfrag) - v_frag(:,i) = vb_frag(:,i) - vcom(:) - end do - - ! Now do a kinetic energy budget check to make sure we are still within the budget. - kefrag = 0.0_DP - do concurrent(i = 1:nfrag) - kefrag(i) = m_frag(i) * dot_product(vb_frag(:, i), vb_frag(:, i)) - end do - ke_frag_orbit = 0.5_DP * sum(kefrag(:)) - - ! If we are over the energy budget, flag this as a failure so we can try again - lerr = ((ke_frag_budget - ke_frag_spin - ke_frag_orbit) < 0.0_DP) - ! write(*,*) 'Tangential' - ! write(*,*) 'Failure? ',lerr - ! call calculate_fragment_ang_mtm() - ! write(*,*) '|L_remainder| : ',.mag.(L_frag_budget(:) - L_frag_tot(:)) / Lmag_before - ! write(*,*) 'ke_frag_budget: ',ke_frag_budget - ! write(*,*) 'ke_frag_spin : ',ke_frag_spin - ! write(*,*) 'ke_tangential : ',ke_frag_orbit - ! write(*,*) 'ke_radial : ',ke_frag_budget - ke_frag_spin - ke_frag_orbit - - return - end subroutine set_fragment_tan_vel - - - function tangential_objective_function(v_t_mag_input, lerr) result(fval) - !! Author: David A. Minton - !! - !! 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) :: v_t_mag_input !! Unknown tangential component of velocity vector set previously by angular momentum constraint - logical, intent(out) :: lerr !! Error flag - ! Result - real(DP) :: fval - ! Internals - integer(I4B) :: i - real(DP), dimension(NDIM,nfrag) :: v_shift - real(DP), dimension(nfrag) :: v_t_new, kearr - real(DP) :: keo - - lerr = .false. - - v_t_new(:) = solve_fragment_tan_vel(v_t_mag_input=v_t_mag_input(:), lerr=lerr) - v_shift(:,:) = vmag_to_vb(v_r_mag, v_r_unit, v_t_new, v_t_unit, m_frag, vcom) - - kearr = 0.0_DP - do concurrent(i = 1:nfrag) - kearr(i) = m_frag(i) * dot_product(v_shift(:, i), v_shift(:, i)) - end do - keo = 0.5_DP * sum(kearr(:)) - fval = keo - lerr = .false. - - return - end function tangential_objective_function - - - function solve_fragment_tan_vel(lerr, v_t_mag_input) result(v_t_mag_output) - !! 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 - logical, intent(out) :: lerr !! Error flag - real(DP), dimension(:), optional, intent(in) :: v_t_mag_input !! Unknown tangential velocities for fragments 7:nfrag - ! Internals - integer(I4B) :: i - ! Result - real(DP), dimension(:), allocatable :: 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 - - lerr = .false. - ! 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 - 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) - A(4:6, i) = m_frag(i) * rmag(i) * (v_r_unit(:, i) .cross. v_t_unit(:, i)) - else if (present(v_t_mag_input)) then - vtmp(:) = v_t_mag_input(i - 6) * v_t_unit(:, i) - L_lin_others(:) = L_lin_others(:) + m_frag(i) * vtmp(:) - L(:) = m_frag(i) * (x_frag(:, i) .cross. vtmp(:)) - L_orb_others(:) = L_orb_others(:) + L(:) - end if - end do - b(1:3) = -L_lin_others(:) - b(4:6) = L_frag_budget(:) - L_frag_spin(:) - L_orb_others(:) - allocate(v_t_mag_output(nfrag)) - v_t_mag_output(1:6) = util_solve_linear_system(A, b, 6, lerr) - if (present(v_t_mag_input)) v_t_mag_output(7:nfrag) = v_t_mag_input(:) - - return - end function solve_fragment_tan_vel - - - subroutine set_fragment_radial_velocities(lerr) - !! 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 - logical, intent(out) :: lerr - ! Internals - real(DP), parameter :: TOL_MIN = Etol ! This needs to be more accurate than the tangential step, as we are trying to minimize the total residual energy - real(DP), parameter :: TOL_INIT = 1e-14_DP - integer(I4B), parameter :: MAXLOOP = 100 - real(DP) :: ke_radial, tol - integer(I4B) :: i, j - real(DP), dimension(:), allocatable :: v_r_initial, v_r_sigma - real(DP), dimension(:,:), allocatable :: v_r - type(lambda_obj) :: objective_function - - ! Set the "target" ke for the radial component - ke_radial = ke_frag_budget - ke_frag_spin - ke_frag_orbit - - allocate(v_r_initial, source=v_r_mag) - ! Initialize radial velocity magnitudes with a random value that is approximately 10% of that found by distributing the kinetic energy equally - allocate(v_r_sigma, source=v_r_mag) - call random_number(v_r_sigma(1:nfrag)) - v_r_sigma(1:nfrag) = sqrt(1.0_DP + 2 * (v_r_sigma(1:nfrag) - 0.5_DP) * 1e-4_DP) - v_r_initial(1:nfrag) = v_r_sigma(1:nfrag) * sqrt(abs(2 * ke_radial) / (m_frag(1:nfrag) * nfrag)) - - ! Initialize the lambda function using a structure constructor that calls the init method - ! Minimize the ke objective function using the BFGS optimizer - objective_function = lambda_obj(radial_objective_function) - tol = TOL_INIT - do while(tol < TOL_MIN) - v_r_mag = util_minimize_bfgs(objective_function, nfrag, v_r_initial, tol, MAXLOOP, lerr) - if (.not.lerr) exit - tol = tol * 2 ! Keep increasing the tolerance until we converge on a solution - v_r_initial(:) = v_r_mag(:) - end do - - ! Shift the radial velocity vectors to align with the center of mass of the collisional system (the origin) - ke_frag_orbit = 0.0_DP - vb_frag(:,1:nfrag) = vmag_to_vb(v_r_mag(1:nfrag), v_r_unit(:,1:nfrag), v_t_mag(1:nfrag), v_t_unit(:,1:nfrag), m_frag(1:nfrag), vcom(:)) - do i = 1, nfrag - v_frag(:, i) = vb_frag(:, i) - vcom(:) - ke_frag_orbit = ke_frag_orbit + m_frag(i) * dot_product(vb_frag(:, i), vb_frag(:, i)) - end do - ke_frag_orbit = 0.5_DP * ke_frag_orbit - - ! write(*,*) 'Radial' - ! write(*,*) 'Failure? ',lerr - ! write(*,*) 'ke_frag_budget: ',ke_frag_budget - ! write(*,*) 'ke_frag_spin : ',ke_frag_spin - ! write(*,*) 'ke_frag_orbit : ',ke_frag_orbit - ! write(*,*) 'ke_remainder : ',ke_frag_budget - (ke_frag_orbit + ke_frag_spin) - lerr = .false. - - return - end subroutine set_fragment_radial_velocities - - - function radial_objective_function(v_r_mag_input) result(fval) - !! Author: David A. Minton - !! - !! 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) :: v_r_mag_input !! Unknown radial component of fragment velocity vector - ! 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 - real(DP), dimension(:,:), allocatable :: v_shift - real(DP), dimension(nfrag) :: kearr - real(DP) :: keo, ke_radial - - allocate(v_shift, mold=vb_frag) - v_shift(:,:) = vmag_to_vb(v_r_mag_input, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom) - do concurrent(i = 1:nfrag) - kearr(i) = m_frag(i) * (Ip_frag(3, i) * rad_frag(i)**2 * dot_product(rot_frag(:, i), rot_frag(:, i)) + dot_product(v_shift(:, i), v_shift(:, i))) - end do - keo = 2 * ke_frag_budget - sum(kearr(:)) - ke_radial = ke_frag_budget - ke_frag_orbit - ke_frag_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 - - return - end function radial_objective_function - - - function vmag_to_vb(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom) result(vb) - !! Author: David A. Minton - !! - !! Converts radial and tangential velocity magnitudes into barycentric velocity - implicit none - ! Arguments - 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) :: m_frag !! Fragment masses - real(DP), dimension(:), intent(in) :: vcom !! Barycentric velocity of collisional system center of mass - ! Result - real(DP), dimension(:,:), allocatable :: vb - ! Internals - integer(I4B) :: i - - allocate(vb, mold=v_r_unit) - ! Make sure the velocity magnitude stays positive - do i = 1, nfrag - vb(:,i) = abs(v_r_mag(i)) * v_r_unit(:, i) - end do - ! 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 - call shift_vector_to_origin(m_frag, vb) - - do i = 1, nfrag - vb(:, i) = vb(:, i) + v_t_mag(i) * v_t_unit(:, i) + vcom(:) - end do - - return - end function vmag_to_vb - - - subroutine restructure_failed_fragments() - !! Author: David A. Minton - !! - !! We failed to find a set of positions and velocities that satisfy all the constraints, and so we will alter the fragments and try again. - implicit none - integer(I4B) :: i - real(DP), dimension(:), allocatable :: m_frag_new, rad_frag_new - real(DP), dimension(:,:), allocatable :: xb_frag_new, vb_frag_new, Ip_frag_new, rot_frag_new - real(DP) :: delta_r, delta_r_max - real(DP), parameter :: ke_avg_deficit_target = 0.0_DP - - ke_tot_deficit = ke_tot_deficit - (ke_frag_budget - ke_frag_orbit - ke_frag_spin) - ke_avg_deficit = ke_tot_deficit / try - ! Introduce a bit of noise in the radius determination so we don't just flip flop between similar failed positions - call random_number(delta_r_max) - delta_r_max = sum(radius(:)) * (1.0_DP + 2e-1_DP * (delta_r_max - 0.5_DP)) - if (try > 1) then - ! Linearly interpolate the last two failed solution ke deficits to find a new distance value to try - delta_r = (r_max_start - r_max_start_old) * (ke_avg_deficit_target - ke_avg_deficit_old) / (ke_avg_deficit - ke_avg_deficit_old) - if (abs(delta_r) > delta_r_max) delta_r = sign(delta_r_max, delta_r) - else - delta_r = delta_r_max - end if - r_max_start_old = r_max_start - r_max_start = r_max_start + delta_r ! The larger lever arm can help if the problem is in the angular momentum step - ke_avg_deficit_old = ke_avg_deficit - - if (f_spin > epsilon(1.0_DP)) then ! Try reducing the fraction in spin - f_spin = f_spin / 2 - else - f_spin = 0.0_DP - end if - - return - end subroutine restructure_failed_fragments - end subroutine fragmentation_initialize - - - module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, min_mfrag, Qloss) - !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Determine the collisional regime of two colliding bodies. - !! Current version requires all values to be converted to SI units prior to calling the function - !! References: - !! Kokubo, E., Genda, H., 2010. Formation of Terrestrial Planets from Protoplanets Under a Realistic Accretion - !! Condition. ApJL 714, L21. https://doi.org/10.1088/2041-8205/714/1/L21 - !! Leinhardt, Z.M., Stewart, S.T., 2012. Collisions between Gravity-dominated Bodies. I. Outcome Regimes and Scaling - !! Laws 745, 79. https://doi.org/10.1088/0004-637X/745/1/79 - !! Mustill, A.J., Davies, M.B., Johansen, A., 2018. The dynamical evolution of transiting planetary systems including - !! a realistic collision prescription. Mon Not R Astron Soc 478, 2896–2908. https://doi.org/10.1093/mnras/sty1273 - !! Rufu, R., Aharonson, O., 2019. Impact Dynamics of Moons Within a Planetary Potential. J. Geophys. Res. Planets 124, - !! 1008–1019. https://doi.org/10.1029/2018JE005798 - !! Stewart, S.T., Leinhardt, Z.M., 2012. Collisions between Gravity-dominated Bodies. II. The Diversity of Impact - !! Outcomes during the End Stage of Planet Formation. ApJ 751, 32. https://doi.org/10.1088/0004-637X/751/1/32 - !! - implicit none - ! Arguments - integer(I4B), intent(out) :: regime - real(DP), intent(out) :: Mlr, Mslr - real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, min_mfrag - real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 - real(DP), intent(out) :: Qloss !! The residual energy after the collision - ! Constants - integer(I4B), parameter :: N1 = 1 !number of objects with mass equal to the largest remnant from LS12 - integer(I4B), parameter :: N2 = 2 !number of objects with mass larger than second largest remnant from LS12 - real(DP), parameter :: DENSITY1 = 1000.0_DP !standard density parameter from LS12 [kg/m3] - real(DP), parameter :: MU_BAR = 0.37_DP !0.385#0.37#0.3333# 3.978 # 1/3 material parameter for hydrodynamic planet-size bodies (LS12) - real(DP), parameter :: BETA = 2.85_DP !slope of sfd for remnants from LS12 2.85 - real(DP), parameter :: C1 = 2.43_DP !! Kokubo & Genda (2010) eq. (3) - real(DP), parameter :: C2 = -0.0408_DP !! Kokubo & Genda (2010) eq. (3) - real(DP), parameter :: C3 = 1.86_DP !! Kokubo & Genda (2010) eq. (3) - real(DP), parameter :: C4 = 1.08_DP !! Kokubo & Genda (2010) eq. (3) - real(DP), parameter :: CRUFU = 2.0_DP - 3 * MU_BAR ! central potential variable from Rufu and Aharonson (2019) - real(DP), parameter :: SUPERCAT_QRATIO = 1.8_DP ! See Section 4.1 of LS12 - ! Internals - real(DP) :: a1, alpha, aint, b, bcrit, c_star, egy, zeta, l, lint, mu, phi, theta - real(DP) :: Qr, Qrd_pstar, Qr_erosion, Qr_supercat - real(DP) :: Vhr, Verosion, Vescp, Vhill, Vimp, Vsupercat - real(DP) :: Mint, Mtot - real(DP) :: Rp, rhill - real(DP) :: Mresidual - real(DP) :: U_binding - - Vimp = norm2(vb2(:) - vb1(:)) - b = calc_b(xh2, vb2, xh1, vb1) - l = (rad1 + rad2) * (1 - b) - egy = 0.5_DP * dot_product(vb1, vb1) - GC * Mcb / norm2(xh1) - a1 = - GC * Mcb / 2.0_DP / egy - Mtot = m1 + m2 - mu = (m1 * m2) / Mtot - if (l < 2 * rad2) then - !calculate Mint - phi = 2 * acos((l - rad2) / rad2) - aint = rad2**2 * (PI - (phi - sin(phi)) / 2.0_DP) - lint = 2 * sqrt(rad2**2 - (rad2 - l / 2.0_DP) ** 2) - Mint = aint * lint ![kg] - alpha = (l**2) * (3 * rad2 - l) / (4 * (rad2**3)) - else - alpha = 1.0_DP - Mint = m2 - end if - Rp = (3 * (m1 / den1 + alpha * m2 / den2) / (4 * PI))**(1.0_DP/3.0_DP) ! (Mustill et al. 2018) - c_star = calc_c_star(Rp) - !calculate Vescp - Vescp = sqrt(2 * GC * Mtot / Rp) !Mustill et al. 2018 eq 6 - !calculate rhill - rhill = a1 * (m1 / 3.0_DP / (Mcb + m1))**(1.0_DP/3.0_DP) - !calculate Vhill - if ((rad2 + rad1) < rhill) then - Vhill = sqrt(2 * GC * m1 * ((rhill**2 - rhill * (rad1 + rad2)) / & - (rhill**2 - 0.5_DP * (rad1 + rad2)**2)) / (rad1 + rad2)) - else - Vhill = Vescp - end if - !calculate Qr_pstar - Qrd_pstar = calc_Qrd_pstar(m1, m2, alpha, c_star) * (Vhill / Vescp)**CRUFU !Rufu and Aharaonson eq (3) - !calculate Verosion - Qr_erosion = 2 * (1.0_DP - m1 / Mtot) * Qrd_pstar - Verosion = (2 * Qr_erosion * Mtot / mu)** (1.0_DP / 2.0_DP) - Qr = mu*(Vimp**2) / Mtot / 2.0_DP - !calculate mass largest remnant Mlr - Mlr = (1.0_DP - Qr / Qrd_pstar / 2.0_DP) * Mtot ! [kg] # LS12 eq (5) - !calculate Vsupercat - Qr_supercat = SUPERCAT_QRATIO * Qrd_pstar ! See LS12 Section 4.1 - Vsupercat = sqrt(2 * Qr_supercat * Mtot / mu) - !calculate Vhr - zeta = (m1 - m2) / Mtot - theta = 1.0_DP - b - Vhr = Vescp * (C1 * zeta**2 * theta**(2.5_DP) + C2 * zeta**2 + C3 * theta**(2.5_DP) + C4) ! Kokubo & Genda (2010) eq. (3) - bcrit = rad1 / (rad1 + rad2) - Qloss = 0.0_DP - U_binding = (3.0_DP * Mtot) / (5.0_DP * Rp) ! LS12 eq. 27 - - if ((m1 < min_mfrag).or.(m2 < min_mfrag)) then - regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime - Mlr = Mtot - Mslr = 0.0_DP - Qloss = 0.0_DP - write(*,*) "FORCE MERGE" - else - if( Vimp < Vescp) then - regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime - Mlr = Mtot - Mslr = 0.0_DP - Qloss = 0.0_DP - else if (Vimp < Verosion) then - if (b < bcrit) then - regime = COLLRESOLVE_REGIME_MERGE !partial accretion regime" - Mlr = Mtot - Mslr = 0.0_DP - Qloss = 0.0_DP - else if ((b > bcrit) .and. (Vimp < Vhr)) then - regime = COLLRESOLVE_REGIME_MERGE ! graze and merge - Mlr = Mtot - Mslr = 0.0_DP - Qloss = 0.0_DP - else - Mlr = m1 - Mslr = calc_Qrd_rev(m2, m1, Mint, den1, den2, Vimp, c_star) - regime = COLLRESOLVE_REGIME_HIT_AND_RUN !hit and run - Qloss = (c_star + 1.0_DP) * U_binding ! Qr - end if - else if (Vimp > Verosion .and. Vimp < Vsupercat) then - if (m2 < 0.001_DP * m1) then - regime = COLLRESOLVE_REGIME_MERGE !cratering regime" - Mlr = Mtot - Mslr = 0.0_DP - Qloss = 0.0_DP - else - Mslr = Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA) ! LS12 eq (37) - regime = COLLRESOLVE_REGIME_DISRUPTION !disruption - Qloss = (c_star + 1.0_DP) * U_binding ! Qr - Qr_erosion - end if - else if (Vimp > Vsupercat) then - Mlr = Mtot * 0.1_DP * (Qr / (Qrd_pstar * SUPERCAT_QRATIO))**(-1.5_DP) !LS12 eq (44) - Mslr = Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA) !LS12 eq (37) - regime = COLLRESOLVE_REGIME_SUPERCATASTROPHIC ! supercatastrophic - Qloss = (c_star + 1.0_DP) * U_binding ! Qr - Qr_supercat - else - write(*,*) "Error no regime found in symba_regime" - end if - end if - Mresidual = Mtot - Mlr - Mslr - if (Mresidual < 0.0_DP) then ! prevents final masses from going negative - Mlr = Mlr + Mresidual - end if - - return - - ! Internal functions - contains - function calc_Qrd_pstar(Mtarg, Mp, alpha, c_star) result(Qrd_pstar) - !! author: Jennifer L.L. Pouplin and Carlisle A. Wishard - !! - !! Calculates the corrected Q* for oblique impacts. See Eq. (15) of LS12. - !! Reference: - !! Leinhardt, Z.M., Stewart, S.T., 2012. Collisions between Gravity-dominated Bodies. I. Outcome Regimes and Scaling - !! Laws 745, 79. https://doi.org/10.1088/0004-637X/745/1/79 - !! - implicit none - ! Arguments - real(DP),intent(in) :: Mtarg, Mp, alpha, c_star - ! Result - real(DP) :: Qrd_pstar - ! Internals - real(DP) :: Qrd_star1, mu_alpha, mu, Qrd_star - - ! calc mu, mu_alpha - mu = (Mtarg * Mp) / (Mtarg + Mp) ! [kg] - mu_alpha = (Mtarg * alpha * Mp) / (Mtarg + alpha * Mp) ! [kg] - ! calc Qrd_star1 - Qrd_star1 = (c_star * 4 * PI * DENSITY1 * GC * Rp**2) / 5.0_DP - ! calc Qrd_star - Qrd_star = Qrd_star1 * (((Mp / Mtarg + 1.0_DP)**2) / (4 * Mp / Mtarg))**(2.0_DP / (3.0_DP * MU_BAR) - 1.0_DP) !(eq 23) - ! calc Qrd_pstar, v_pstar - Qrd_pstar = ((mu / mu_alpha)**(2.0_DP - 3.0_DP * MU_BAR / 2.0_DP)) * Qrd_star ! (eq 15) - - return - end function calc_Qrd_pstar - - function calc_Qrd_rev(Mp, Mtarg, Mint, den1, den2, Vimp, c_star) result(Mslr) - !! author: Jennifer L.L. Pouplin and Carlisle A. Wishard - !! - !! Calculates mass of second largest fragment. - !! - implicit none - ! Arguments - real(DP),intent(in) :: Mp, Mtarg, Mint, den1, den2, Vimp, c_star - ! Result - real(DP) :: Mslr - ! Internals - real(DP) :: mtot_rev, mu_rev, gamma_rev, Qrd_star1, Qrd_star, mu_alpha_rev - real(DP) :: Qrd_pstar, Rc1, Qr_rev, Qrd_pstar_rev, Qr_supercat_rev - - ! calc Mslr, Rc1, mu, gammalr - mtot_rev = Mint + Mp - Rc1 = (3 * (Mint / den1 + Mp / den2) / (4 * PI))**(1.0_DP/3.0_DP) ! [m] Mustill et al 2018 - mu_rev = (Mint * Mp) / mtot_rev ! [kg] eq 49 LS12 - mu_alpha_rev = (Mtarg * alpha * Mp) / (Mtarg + alpha * Mp) - gamma_rev = Mint / Mp ! eq 50 LS12 - !calc Qr_rev - Qr_rev = mu_rev * (Vimp**2) / (2 * mtot_rev) - ! calc Qrd_star1, v_star1 - Qrd_star1 = (c_star * 4 * PI * mtot_rev * GC ) / Rc1 / 5.0_DP - ! calc Qrd_pstar_rev - Qrd_star = Qrd_star1 * (((gamma_rev + 1.0_DP)**2) / (4 * gamma_rev)) ** (2.0_DP / (3.0_DP * MU_BAR) - 1.0_DP) !(eq 52) - Qrd_pstar = Qrd_star * ((mu_rev / mu_alpha_rev)**(2.0_DP - 3.0_DP * MU_BAR / 2.0_DP)) - Qrd_pstar_rev = Qrd_pstar * (Vhill / Vescp)**CRUFU !Rufu and Aharaonson eq (3) - !calc Qr_supercat_rev - Qr_supercat_rev = 1.8_DP * Qrd_pstar_rev - if (Qr_rev > Qr_supercat_rev ) then - Mslr = mtot_rev * (0.1_DP * ((Qr_rev / (Qrd_pstar_rev * 1.8_DP))**(-1.5_DP))) !eq (44) - else if ( Qr_rev < Qrd_pstar_rev ) then - Mslr = Mp - else - Mslr = (1.0_DP - Qr_rev / Qrd_pstar_rev / 2.0_DP) * (mtot_rev) ! [kg] #(eq 5) - end if - - if ( Mslr > Mp ) Mslr = Mp !check conservation of mass - - return - end function calc_Qrd_rev - - function calc_b(proj_pos, proj_vel, targ_pos, targ_vel) result(sintheta) - !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Calculates the impact factor b = sin(theta), where theta is the angle between the relative velocity - !! and distance vectors of the target and projectile bodies. See Fig. 2 of Leinhardt and Stewart (2012) - !! - implicit none - !! Arguments - real(DP), dimension(:), intent(in) :: proj_pos, proj_vel, targ_pos, targ_vel - !! Result - real(DP) :: sintheta - !! Internals - real(DP), dimension(NDIM) :: imp_vel, distance, x_cross_v - - imp_vel(:) = proj_vel(:) - targ_vel(:) - distance(:) = proj_pos(:) - targ_pos(:) - x_cross_v(:) = distance(:) .cross. imp_vel(:) - sintheta = norm2(x_cross_v(:)) / norm2(distance(:)) / norm2(imp_vel(:)) - return - end function calc_b - - function calc_c_star(Rc1) result(c_star) - !! author: David A. Minton - !! - !! Calculates c_star as a function of impact equivalent radius. It inteRpolates between 5 for ~1 km sized bodies to - !! 1.8 for ~10000 km sized bodies. See LS12 Fig. 4 for details. - !! - implicit none - !! Arguments - real(DP), intent(in) :: Rc1 - !! Result - real(DP) :: c_star - !! Internals - real(DP), parameter :: loR = 1.0e3_DP ! Lower bound of inteRpolation size (m) - real(DP), parameter :: hiR = 1.0e7_DP ! Upper bound of inteRpolation size (m) - real(DP), parameter :: loval = 5.0_DP ! Value of C* at lower bound - real(DP), parameter :: hival = 1.9_DP ! Value of C* at upper bound - - if (Rc1 < loR) then - c_star = loval - else if (Rc1 < hiR) then - c_star = loval + (hival - loval) * log(Rc1 / loR) / log(hiR /loR) - else - c_star = hival - end if - return - end function calc_c_star - - end subroutine fragmentation_regime - -end submodule s_fragmentation \ No newline at end of file diff --git a/src/helio/helio_setup.f90 b/src/helio/helio_setup.f90 index 2a8c24019..cf5f57d8a 100644 --- a/src/helio/helio_setup.f90 +++ b/src/helio/helio_setup.f90 @@ -15,6 +15,8 @@ module subroutine helio_setup_initialize_system(self, param) call whm_setup_initialize_system(self, param) call self%pl%h2b(self%cb) call self%tp%h2b(self%cb) + call self%pl%sort("mass", ascending=.false.) + call self%pl%index(param) return end subroutine helio_setup_initialize_system diff --git a/src/helio/helio_util.f90 b/src/helio/helio_util.f90 deleted file mode 100644 index 3b5ee6116..000000000 --- a/src/helio/helio_util.f90 +++ /dev/null @@ -1,21 +0,0 @@ -submodule(helio_classes) s_helio_eucl - use swiftest -contains - - module subroutine helio_util_index_eucl_plpl(self, param) - !! author: David A. Minton - !! - !! Wrapper for the indexing method for WHM massive bodies. Sorts the massive bodies by heliocentric distance and then flattens the pl-pl upper triangular matrix - implicit none - ! Arguments - class(helio_pl), intent(inout) :: self !! Helio massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - - call self%sort("mass", ascending=.false.) - call util_index_eucl_plpl(self, param) - - return - end subroutine helio_util_index_eucl_plpl - -end submodule s_helio_eucl - diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 new file mode 100644 index 000000000..d6ea0ac32 --- /dev/null +++ b/src/modules/fraggle_classes.f90 @@ -0,0 +1,250 @@ +module fraggle_classes + !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott + !! + !! Definition of classes and methods specific to Fraggel: The Fragment Generation Model + use swiftest_globals + use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system, swiftest_cb, swiftest_pl + implicit none + public + + integer(I4B), parameter :: FRAGGLE_NMASS_DIST = 3 !! Number of mass bins returned by the regime calculation (largest fragment, second largest, and remainder) + !******************************************************************************************************************************** + ! fraggle_colliders class definitions and method interfaces + !******************************************************************************************************************************* + !> Class definition for the variables that describe the bodies involved in the collision + type :: fraggle_colliders + integer(I4B) :: ncoll !! Number of bodies involved in the collision + integer(I4B), dimension(:), allocatable :: idx !! Index of bodies involved in the collision + real(DP), dimension(NDIM,2) :: xb !! Two-body equivalent position vectors of the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: vb !! Two-body equivalent velocity vectors of the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: rot !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: L_spin !! Two-body equivalent spin angular momentum vectors of the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: L_orbit !! Two-body equivalent orbital angular momentum vectors of the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: Ip !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision + real(DP), dimension(2) :: mass !! Two-body equivalent mass of the collider bodies prior to the collision + real(DP), dimension(2) :: radius !! Two-body equivalent radii of the collider bodies prior to the collision + contains + procedure :: regime => fraggle_regime_colliders !! Determine which fragmentation regime the set of colliders will be + end type fraggle_colliders + + !******************************************************************************************************************************** + ! fraggle_fragments class definitions and method interfaces + !******************************************************************************************************************************* + !> Class definition for the variables that describe a collection of fragments by Fraggle barycentric coordinates + type, extends(swiftest_pl) :: fraggle_fragments + real(DP) :: mtot !! Total mass of fragments + real(DP) :: Qloss !! Energy lost during the collision + real(DP), dimension(FRAGGLE_NMASS_DIST) :: mass_dist !! Distribution of fragment mass determined by the regime calculation (largest fragment, second largest, and remainder) + integer(I4B) :: regime !! Collresolve regime code for this collision + + ! Values in a coordinate frame centered on the collider barycenter and collisional system unit vectors (these are used internally by the fragment generation subroutine) + real(DP), dimension(NDIM) :: xbcom !! Center of mass position vector of the collider system in system barycentric coordinates + real(DP), dimension(NDIM) :: vbcom !! Velocity vector of the center of mass of the collider system in system barycentric coordinates + real(DP), dimension(NDIM) :: x_coll_unit !! x-direction unit vector of collisional system + real(DP), dimension(NDIM) :: y_coll_unit !! y-direction unit vector of collisional system + real(DP), dimension(NDIM) :: z_coll_unit !! z-direction unit vector of collisional system + real(DP), dimension(:,:), allocatable :: x_coll !! Array of fragment position vectors in the collisional coordinate frame + real(DP), dimension(:,:), allocatable :: v_coll !! Array of fragment velocity vectors in the collisional coordinate frame + real(DP), dimension(:,:), allocatable :: v_r_unit !! Array of radial direction unit vectors of individual fragments in the collisional coordinate frame + real(DP), dimension(:,:), allocatable :: v_t_unit !! Array of tangential direction unit vectors of individual fragments in the collisional coordinate frame + real(DP), dimension(:,:), allocatable :: v_n_unit !! Array of normal direction unit vectors of individual fragments in the collisional coordinate frame + real(DP), dimension(:), allocatable :: rmag !! Array of radial distance magnitudes of individual fragments in the collisional coordinate frame + real(DP), dimension(:), allocatable :: rotmag !! Array of rotation magnitudes of individual fragments + real(DP), dimension(:), allocatable :: v_r_mag !! Array of radial direction velocity magnitudes of individual fragments + real(DP), dimension(:), allocatable :: v_t_mag !! Array of tangential direction velocity magnitudes of individual fragments + + ! Energy and momentum book-keeping variables that characterize the whole system of fragments + real(DP) :: ke_orbit !! Current orbital kinetic energy of the system of fragments in the collisional frame + real(DP) :: ke_spin !! Current spin kinetic energy of the system of fragments in the collisional frame + real(DP), dimension(NDIM) :: L_orbit !! Current orbital angular momentum of the system of fragments in the collisional frame + real(DP), dimension(NDIM) :: L_spin !! Current spin angular momentum of the system of fragments in the collisional frame + real(DP) :: ke_budget !! Total kinetic energy budget for the system of fragmens in the collisional frame + real(DP), dimension(NDIM) :: L_budget !! Total angular momentum budget for the system of fragmens in the collisional frame + + ! For the following variables, "before" refers to the *entire* n-body system in its pre-collisional state and "after" refers to the system in its post-collisional state + real(DP), dimension(NDIM) :: Lorbit_before, Lorbit_after !! Before/after orbital angular momentum + real(DP), dimension(NDIM) :: Lspin_before, Lspin_after !! Before/after spin angular momentum + real(DP), dimension(NDIM) :: Ltot_before, Ltot_after !! Before/after total system angular momentum + real(DP) :: ke_orbit_before, ke_orbit_after !! Before/after orbital kinetic energy + real(DP) :: ke_spin_before, ke_spin_after !! Before/after spin kinetic energy + real(DP) :: pe_before, pe_after !! Before/after potential energy + real(DP) :: Etot_before, Etot_after !! Before/after total system energy + + ! 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 !! Distance dimension scale factor + real(DP) :: mscale !! Mass scale factor + real(DP) :: tscale !! Time scale factor + real(DP) :: vscale !! Velocity scale factor (a convenience unit that is derived from dscale and tscale) + real(DP) :: Escale !! Energy scale factor (a convenience unit that is derived from dscale, tscale, and mscale) + real(DP) :: Lscale !! 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. + procedure :: accel => fraggle_placeholder_accel !! Placeholder subroutine to fulfill requirement for an accel method + procedure :: kick => fraggle_placeholder_kick !! Placeholder subroutine to fulfill requirement for a kick method + procedure :: step => fraggle_placeholder_step !! Placeholder subroutine to fulfill requirement for a step method + procedure :: set_budgets => fraggle_set_budgets_fragments !! Sets the energy and momentum budgets of the fragments based on the collider value + procedure :: set_coordinate_system => fraggle_set_coordinate_system !! Defines the collisional coordinate system, including the unit vectors of both the system and individual fragments. + 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_ang_mtm => fraggle_util_ang_mtm !! Calcualtes the current angular momentum of the fragments + procedure :: get_energy_and_momentum => fraggle_util_get_energy_momentum !! Calculates total system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) + procedure :: restructure => fraggle_util_restructure !! Restructure the inputs after a failed attempt failed to find a set of positions and velocities that satisfy the energy and momentum constraints + end type fraggle_fragments + + interface + module subroutine fraggle_generate_fragments(self, colliders, system, param, lfailure) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object containing the two-body equivalent values of the colliding bodies + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? + end subroutine fraggle_generate_fragments + + !> The following interfaces are placeholders intended to satisfy the required abstract methods given by the parent class + module subroutine fraggle_placeholder_accel(self, system, param, t, lbeg) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current simulation time + logical, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step + end subroutine fraggle_placeholder_accel + + module subroutine fraggle_placeholder_kick(self, system, param, t, dt, lbeg) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system objec + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current time + real(DP), intent(in) :: dt !! Stepsize + logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. + end subroutine fraggle_placeholder_kick + + module subroutine fraggle_placeholder_step(self, system, param, t, dt) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(fraggle_fragments), intent(inout) :: self !! Helio massive body particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nboody system + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current simulation time + real(DP), intent(in) :: dt !! Stepsiz + end subroutine fraggle_placeholder_step + + module subroutine fraggle_regime_colliders(self, frag, system, param) + implicit none + class(fraggle_colliders), intent(inout) :: self !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragment system object + class(swiftest_nbody_system), intent(in) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameter + end subroutine fraggle_regime_colliders + + module subroutine fraggle_set_budgets_fragments(self, colliders) + implicit none + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + end subroutine fraggle_set_budgets_fragments + + module subroutine fraggle_set_coordinate_system(self, colliders) + implicit none + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + end subroutine fraggle_set_coordinate_system + + module subroutine fraggle_set_mass_dist_fragments(self, colliders) + implicit none + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + end subroutine fraggle_set_mass_dist_fragments + + module subroutine fraggle_set_natural_scale_factors(self, colliders) + implicit none + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + end subroutine fraggle_set_natural_scale_factors + + module subroutine fraggle_set_original_scale_factors(self, colliders) + implicit none + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + end subroutine fraggle_set_original_scale_factors + + module subroutine fraggle_setup_fragments(self, n, param) + implicit none + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + integer(I4B), intent(in) :: n !! Number of fragments + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + end subroutine fraggle_setup_fragments + + module subroutine fraggle_setup_reset_fragments(self) + implicit none + class(fraggle_fragments), intent(inout) :: self + end subroutine fraggle_setup_reset_fragments + + module subroutine fraggle_util_add_fragments_to_system(frag, colliders, system, param) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(fraggle_fragments), intent(in) :: frag !! Fraggle fragment system object + class(fraggle_colliders), intent(in) :: colliders !! Fraggle collider system object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + end subroutine fraggle_util_add_fragments_to_system + + module subroutine fraggle_util_ang_mtm(self) + implicit none + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + end subroutine fraggle_util_ang_mtm + + module subroutine fraggle_util_construct_temporary_system(frag, system, param, tmpsys, tmpparam) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(fraggle_fragments), intent(in) :: frag !! 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 + end subroutine fraggle_util_construct_temporary_system + + module subroutine fraggle_util_get_energy_momentum(self, colliders, system, param, lbefore) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the system, with colliders included and fragments excluded or vice versa + end subroutine fraggle_util_get_energy_momentum + + module subroutine fraggle_util_restructure(self, colliders, try, f_spin, r_max_start) + implicit none + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_colliders), intent(in) :: colliders !! Fraggle collider system object + integer(I4B), intent(in) :: try !! The current number of times Fraggle has tried to find a solution + real(DP), intent(inout) :: f_spin !! Fraction of energy/momentum that goes into spin. This decreases ater a failed attempt + real(DP), intent(inout) :: r_max_start !! The maximum radial distance that the position calculation starts with. This increases after a failed attempt + end subroutine fraggle_util_restructure + + module subroutine fraggle_util_shift_vector_to_origin(m_frag, vec_frag) + implicit none + 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 + end subroutine + + module function fraggle_util_vmag_to_vb(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom) result(vb) + implicit none + 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) :: m_frag !! Fragment masses + real(DP), dimension(:), intent(in) :: vcom !! Barycentric velocity of collisional system center of mass + real(DP), dimension(:,:), allocatable :: vb + end function fraggle_util_vmag_to_vb + end interface + +end module fraggle_classes \ No newline at end of file diff --git a/src/modules/helio_classes.f90 b/src/modules/helio_classes.f90 index dcfcdde2e..1c93bdf6b 100644 --- a/src/modules/helio_classes.f90 +++ b/src/modules/helio_classes.f90 @@ -2,7 +2,7 @@ module helio_classes !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! !! Definition of classes and methods specific to the Democratic Heliocentric Method - !! Adapted from David E. Kaufmann's Swifter routine: helio.f90 + !! Adapted from David E. Kaufmann's Swifter routine: module_helio.f90 use swiftest_globals use swiftest_classes, only : swiftest_cb, swiftest_pl, swiftest_tp, swiftest_nbody_system use whm_classes, only : whm_nbody_system @@ -38,7 +38,6 @@ module helio_classes contains procedure :: drift => helio_drift_pl !! Method for Danby drift in Democratic Heliocentric coordinates procedure :: lindrift => helio_drift_linear_pl !! Method for linear drift of massive bodies due to barycentric momentum of Sun - procedure :: index => helio_util_index_eucl_plpl !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix procedure :: accel_gr => helio_gr_kick_getacch_pl !! Acceleration term arising from the post-Newtonian correction procedure :: gr_pos_kick => helio_gr_p4_pl !! Position kick due to p**4 term in the post-Newtonian correction procedure :: accel => helio_kick_getacch_pl !! Compute heliocentric accelerations of massive bodies @@ -214,12 +213,6 @@ module subroutine helio_step_tp(self, system, param, t, dt) real(DP), intent(in) :: dt !! Stepsizee end subroutine helio_step_tp - module subroutine helio_util_index_eucl_plpl(self, param) - use swiftest_classes, only : swiftest_parameters - implicit none - class(helio_pl), intent(inout) :: self !! Helio massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine helio_util_index_eucl_plpl end interface end module helio_classes diff --git a/src/modules/swiftest.f90 b/src/modules/swiftest.f90 index 61d45163c..a6daff9a7 100644 --- a/src/modules/swiftest.f90 +++ b/src/modules/swiftest.f90 @@ -10,6 +10,7 @@ module swiftest use rmvs_classes use helio_classes use symba_classes + use fraggle_classes use lambda_function !use advisor_annotate !$ use omp_lib diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 746781117..2694fcc25 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -573,31 +573,6 @@ module subroutine util_index_eucl_pltp(self, pl, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine - module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & - nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) - implicit none - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - integer(I4B), dimension(:), intent(in) :: family !! Index of bodies involved in the collision - real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Two-body equivalent position, vector, spin momentum, and rotational inertia values for the collision - real(DP), dimension(:), intent(inout) :: mass, radius !! Two-body equivalent mass and radii for the bodies in the collision - integer(I4B), intent(inout) :: nfrag !! Number of fragments to generate - real(DP), dimension(:), allocatable, intent(inout) :: m_frag, rad_frag !! Distribution of fragment mass and radii - real(DP), dimension(:,:), allocatable, intent(inout) :: Ip_frag !! Fragment rotational inertia vectors - real(DP), dimension(:,:), allocatable, intent(inout) :: xb_frag, vb_frag, rot_frag !! Fragment barycentric position, barycentric velocity, and rotation vectors - real(DP), intent(inout) :: Qloss !! Energy lost during the collision - logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? - end subroutine fragmentation_initialize - - module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, min_mfrag, Qloss) - implicit none - integer(I4B), intent(out) :: regime - real(DP), intent(out) :: Mlr, Mslr - real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, min_mfrag - real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 - real(DP), intent(out) :: Qloss !! Energy lost during the collision - end subroutine fragmentation_regime - module pure subroutine gr_kick_getaccb_ns_body(self, system, param) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest generic body object @@ -1048,8 +1023,8 @@ end subroutine setup_initialize_system module subroutine setup_pl(self, n, param) implicit none - class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object - integer(I4B), intent(in) :: n !! Number of particles to allocate space for + class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object + integer(I4B), intent(in) :: n !! Number of particles to allocate space for class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine setup_pl diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 75d8e0e73..b929b5040 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -1,12 +1,13 @@ module symba_classes !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! - !! Definition of classes and methods specific to the Democratic SyMBAcentric Method - !! Adapted from David E. Kaufmann's Swifter routine: helio.f90 + !! Definition of classes and methods specific to the SyMBA integrator + !! Adapted from David E. Kaufmann's Swifter routine: module_symba.f90 use swiftest_globals use swiftest_classes, only : swiftest_parameters, swiftest_base, swiftest_encounter, swiftest_particle_info, netcdf_parameters use helio_classes, only : helio_cb, helio_pl, helio_tp, helio_nbody_system use rmvs_classes, only : rmvs_chk_ind + use fraggle_classes, only : fraggle_colliders, fraggle_fragments implicit none public @@ -66,7 +67,7 @@ module symba_classes real(DP), dimension(:), allocatable :: atp !! semimajor axis following perihelion passage type(symba_kinship), dimension(:), allocatable :: kin !! Array of merger relationship structures that can account for multiple pairwise mergers in a single step contains - procedure :: make_family => symba_collision_make_family_pl !! When a single body is involved in more than one collision in a single step, it becomes part of a family + procedure :: make_colliders => symba_collision_make_colliders_pl !! When a single body is involved in more than one collision in a single step, it becomes part of a family procedure :: index => symba_util_index_eucl_plpl !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix procedure :: discard => symba_discard_pl !! Process massive body discards procedure :: drift => symba_drift_pl !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level @@ -190,11 +191,11 @@ module subroutine symba_collision_encounter_extract_collisions(self, system, par class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine - module subroutine symba_collision_make_family_pl(self,idx) + module subroutine symba_collision_make_colliders_pl(self,idx) implicit none class(symba_pl), intent(inout) :: self !! SyMBA massive body object integer(I4B), dimension(2), intent(in) :: idx !! Array holding the indices of the two bodies involved in the collision - end subroutine symba_collision_make_family_pl + end subroutine symba_collision_make_colliders_pl module subroutine symba_collision_resolve_fragmentations(self, system, param) implicit none @@ -291,50 +292,44 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc logical :: lany_encounter !! Returns true if there is at least one close encounter end function symba_encounter_check_tp - module function symba_collision_casedisruption(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) + module function symba_collision_casedisruption(system, param, colliders, frag) result(status) + use fraggle_classes, only : fraggle_colliders, fraggle_fragments implicit none - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision - real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision - real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision - real(DP), dimension(:), intent(inout) :: mass_res !! The distribution of fragment mass obtained by the regime calculation - real(DP), intent(inout) :: Qloss !! Energy lost during collisionn - integer(I4B) :: status !! Status flag assigned to this outcome + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object + integer(I4B) :: status !! Status flag assigned to this outcome end function symba_collision_casedisruption - module function symba_collision_casehitandrun(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) - implicit none - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision - real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision - real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision - real(DP), dimension(:), intent(inout) :: mass_res !! The distribution of fragment mass obtained by the regime calculation - real(DP), intent(inout) :: Qloss !! Energy lost during collision - integer(I4B) :: status !! Status flag assigned to this outcome + module function symba_collision_casehitandrun(system, param, colliders, frag) result(status) + use fraggle_classes, only : fraggle_colliders, fraggle_fragments + implicit none + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object + integer(I4B) :: status !! Status flag assigned to this outcome end function symba_collision_casehitandrun - module function symba_collision_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) result(status) + module function symba_collision_casemerge(system, param, colliders, frag) result(status) + use fraggle_classes, only : fraggle_colliders, fraggle_fragments implicit none - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision - real(DP), dimension(:,:), intent(in) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision - real(DP), dimension(:), intent(in) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collisio - integer(I4B) :: status !! Status flag assigned to this outcome + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object + integer(I4B) :: status !! Status flag assigned to this outcome end function symba_collision_casemerge - module function symba_collision_casesupercatastrophic(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) + module function symba_collision_casesupercatastrophic(system, param, colliders, frag) result(status) + use fraggle_classes, only : fraggle_colliders, fraggle_fragments implicit none - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision - real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision - real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision - real(DP), dimension(:), intent(inout) :: mass_res !! The distribution of fragment mass obtained by the regime calculation - real(DP), intent(inout) :: Qloss !! Energy lost during collision - integer(I4B) :: status !! Status flag assigned to this outcome + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object + integer(I4B) :: status !! Status flag assigned to this outcome end function symba_collision_casesupercatastrophic module subroutine symba_util_index_eucl_plpl(self, param) diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90 index 1b8d13b45..9675f289d 100644 --- a/src/modules/whm_classes.f90 +++ b/src/modules/whm_classes.f90 @@ -34,7 +34,6 @@ module whm_classes procedure :: j2h => whm_coord_j2h_pl !! Convert position and velcoity vectors from Jacobi to helliocentric coordinates procedure :: vh2vj => whm_coord_vh2vj_pl !! Convert velocity vectors from heliocentric to Jacobi coordinates procedure :: drift => whm_drift_pl !! Loop through massive bodies and call Danby drift routine to jacobi coordinates - procedure :: index => whm_util_index_eucl_plpl !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix procedure :: accel_gr => whm_gr_kick_getacch_pl !! Acceleration term arising from the post-Newtonian correction procedure :: gr_pos_kick => whm_gr_p4_pl !! Position kick due to p**4 term in the post-Newtonian correction procedure :: accel => whm_kick_getacch_pl !! Compute heliocentric accelerations of massive bodies diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 504faca35..41c0d9fba 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -81,7 +81,6 @@ module subroutine setup_encounter(self, n) class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter structure integer(I4B), intent(in) :: n !! Number of encounters to allocate space for - self%nenc = n if (n < 0) return if (allocated(self%lvdotr)) deallocate(self%lvdotr) @@ -97,6 +96,7 @@ module subroutine setup_encounter(self, n) if (allocated(self%v2)) deallocate(self%v2) if (allocated(self%t)) deallocate(self%t) + self%nenc = n if (n == 0) return allocate(self%lvdotr(n)) @@ -206,7 +206,6 @@ module subroutine setup_body(self, n, param) ! Internals integer(I4B) :: i - self%nbody = n if (n < 0) return self%lfirst = .true. @@ -234,6 +233,7 @@ module subroutine setup_body(self, n, param) if (allocated(self%omega)) deallocate(self%omega) if (allocated(self%capm)) deallocate(self%capm) + self%nbody = n if (n == 0) return allocate(self%info(n)) diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 384102ffc..a3bf27e99 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -5,81 +5,37 @@ integer(I4B), parameter :: NFRAG_SUPERCAT = 20 contains - module function symba_collision_casedisruption(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) + module function symba_collision_casedisruption(system, param, colliders, frag) result(status) !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Create the fragments resulting from a non-catastrophic disruption collision !! implicit none ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision - real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision - real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision - real(DP), dimension(:), intent(inout) :: mass_res !! The distribution of fragment mass obtained by the regime calculation - real(DP), intent(inout) :: Qloss !! Energy lost during collision + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object ! Result - integer(I4B) :: status !! Status flag assigned to this outcome + integer(I4B) :: status !! Status flag assigned to this outcome ! Internals - integer(I4B) :: i, istart, nfrag, nfamily, nstart, nend - real(DP) :: mtot, avg_dens - real(DP), dimension(NDIM) :: xcom, vcom, Ip_new - real(DP), dimension(2) :: vol - real(DP), dimension(:, :), allocatable :: vb_frag, xb_frag, rot_frag, Ip_frag - real(DP), dimension(:), allocatable :: m_frag, rad_frag - integer(I4B), dimension(:), allocatable :: id_frag - logical :: lfailure + integer(I4B) :: i, nfrag + logical :: lfailure character(len=STRMAX) :: collider_message collider_message = "Disruption between" - call symba_collision_collider_message(system%pl, family, collider_message) + call symba_collision_collider_message(system%pl, colliders%idx, collider_message) write(*,*) trim(adjustl(collider_message)) ! Collisional fragments will be uniformly distributed around the pre-impact barycenter nfrag = NFRAG_DISRUPT - allocate(m_frag(nfrag)) - allocate(rad_frag(nfrag)) - allocate(xb_frag(NDIM, nfrag)) - allocate(vb_frag(NDIM, nfrag)) - allocate(rot_frag(NDIM, nfrag)) - allocate(Ip_frag(NDIM, nfrag)) - allocate(id_frag(nfrag)) - - mtot = sum(mass(:)) - xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot - vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot - - ! Get mass weighted mean of Ip and average density - Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mtot - vol(:) = 4._DP / 3._DP * PI * radius(:)**3 - avg_dens = mtot / sum(vol(:)) - - ! Distribute the mass among fragments, with a branch to check for the size of the second largest fragment - m_frag(1) = mass_res(1) - if (mass_res(2) > mass_res(1) / 3._DP) then - m_frag(2) = mass_res(2) - istart = 3 - else - istart = 2 - end if - ! Distribute remaining mass among the remaining bodies - do i = istart, nfrag - m_frag(i) = (mtot - sum(m_frag(1:istart - 1))) / (nfrag - istart + 1) - end do + call frag%setup(nfrag, param) + call frag%set_mass_dist(colliders) + frag%id(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)] + param%maxid = frag%id(nfrag) - ! Distribute any residual mass if there is any and set the radius - m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) - rad_frag(1:nfrag) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) - id_frag(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)] - param%maxid = id_frag(nfrag) - - do i = 1, nfrag - Ip_frag(:, i) = Ip_new(:) - end do - - call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & - nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) + ! Generate the position and velocity distributions of the fragments + call frag%generate_fragments(colliders, system, param, lfailure) if (lfailure) then write(*,*) 'No fragment solution found, so treat as a pure hit-and-run' @@ -87,61 +43,45 @@ module function symba_collision_casedisruption(system, param, family, x, v, mass nfrag = 0 select type(pl => system%pl) class is (symba_pl) - pl%status(family(:)) = status - pl%ldiscard(family(:)) = .false. - pl%lcollision(family(:)) = .false. + pl%status(colliders%idx(:)) = status + pl%ldiscard(colliders%idx(:)) = .false. + pl%lcollision(colliders%idx(:)) = .false. end select else ! Populate the list of new bodies write(*,'("Generating ",I2.0," fragments")') nfrag status = DISRUPTION - call symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, status) + call symba_collision_mergeaddsub(system, param, colliders, frag, status) end if return end function symba_collision_casedisruption - module function symba_collision_casehitandrun(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) + module function symba_collision_casehitandrun(system, param, colliders, frag) result(status) !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Create the fragments resulting from a non-catastrophic hit-and-run collision !! implicit none ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision - real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision - real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision - real(DP), dimension(:), intent(inout) :: mass_res !! The distribution of fragment mass obtained by the regime calculation - real(DP), intent(inout) :: Qloss !! Energy lost during collision + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object ! Result - integer(I4B) :: status !! Status flag assigned to this outcome + integer(I4B) :: status !! Status flag assigned to this outcom ! Internals - integer(I4B) :: i, j, nfrag, jproj, jtarg, idstart, ibiggest, nfamily - real(DP) :: mtot, avg_dens - real(DP), dimension(NDIM) :: xcom, vcom - real(DP), dimension(2) :: vol - real(DP), dimension(:, :), allocatable :: vb_frag, xb_frag, rot_frag, Ip_frag - real(DP), dimension(:), allocatable :: m_frag, rad_frag - integer(I4B), dimension(:), allocatable :: id_frag - logical :: lpure - logical, dimension(system%pl%nbody) :: lmask + integer(I4B) :: i, ibiggest, nfrag, jtarg, jproj + logical :: lpure character(len=STRMAX) :: collider_message character(len=NAMELEN) :: idstr collider_message = "Hit and run between" - call symba_collision_collider_message(system%pl, family, collider_message) + call symba_collision_collider_message(system%pl, colliders%idx, collider_message) write(*,*) trim(adjustl(collider_message)) - mtot = sum(mass(:)) - xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot - vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot - lpure = .false. - - ! The largest body will stay untouched - if (mass(1) > mass(2)) then + if (colliders%mass(1) > colliders%mass(2)) then jtarg = 1 jproj = 2 else @@ -149,44 +89,23 @@ module function symba_collision_casehitandrun(system, param, family, x, v, mass, jproj = 1 end if - if (mass_res(2) > 0.9_DP * mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched + if (frag%mass_dist(2) > 0.9_DP * colliders%mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched write(*,*) 'Pure hit and run. No new fragments generated.' nfrag = 0 lpure = .true. else ! Imperfect hit and run, so we'll keep the largest body and destroy the other nfrag = NFRAG_DISRUPT - 1 lpure = .false. - allocate(m_frag(nfrag)) - allocate(id_frag(nfrag)) - allocate(rad_frag(nfrag)) - allocate(xb_frag(NDIM, nfrag)) - allocate(vb_frag(NDIM, nfrag)) - allocate(rot_frag(NDIM, nfrag)) - allocate(Ip_frag(NDIM, nfrag)) - m_frag(1) = mass(jtarg) - ibiggest = family(maxloc(system%pl%Gmass(family(:)), dim=1)) - id_frag(1) = system%pl%id(ibiggest) - rad_frag(1) = radius(jtarg) - xb_frag(:, 1) = x(:, jtarg) - vb_frag(:, 1) = v(:, jtarg) - Ip_frag(:,1) = Ip(:, jtarg) - - ! Get mass weighted mean of Ip and average density - vol(:) = 4._DP / 3._DP * pi * radius(:)**3 - avg_dens = mass(jproj) / vol(jproj) - m_frag(2:nfrag) = (mtot - m_frag(1)) / (nfrag - 1) - rad_frag(2:nfrag) = (3 * m_frag(2:nfrag) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) - m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) - id_frag(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] - param%maxid = id_frag(nfrag) - - do i = 1, nfrag - Ip_frag(:, i) = Ip(:, jproj) - end do + call frag%setup(nfrag, param) + call frag%set_mass_dist(colliders) + ibiggest = colliders%idx(maxloc(system%pl%Gmass(colliders%idx(:)), dim=1)) + frag%id(1) = system%pl%id(ibiggest) + frag%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] + param%maxid = frag%id(nfrag) + + ! Generate the position and velocity distributions of the fragments + call frag%generate_fragments(colliders, system, param, lpure) - ! Put the fragments on the circle surrounding the center of mass of the system - call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & - 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 @@ -198,20 +117,20 @@ module function symba_collision_casehitandrun(system, param, family, x, v, mass, status = HIT_AND_RUN_PURE select type(pl => system%pl) class is (symba_pl) - pl%status(family(:)) = ACTIVE - pl%ldiscard(family(:)) = .false. - pl%lcollision(family(:)) = .false. + pl%status(colliders%idx(:)) = ACTIVE + pl%ldiscard(colliders%idx(:)) = .false. + pl%lcollision(colliders%idx(:)) = .false. end select else status = HIT_AND_RUN_DISRUPT - call symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, status) + call symba_collision_mergeaddsub(system, param, colliders, frag, status) end if return end function symba_collision_casehitandrun - module function symba_collision_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) result(status) + module function symba_collision_casemerge(system, param, colliders, frag) result(status) !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Merge planets. @@ -221,76 +140,48 @@ module function symba_collision_casemerge(system, param, family, x, v, mass, rad !! Adapted from Hal Levison's Swift routines symba5_merge.f and discard_mass_merge.f implicit none ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision - real(DP), dimension(:,:), intent(in) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision - real(DP), dimension(:), intent(in) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object ! Result - integer(I4B) :: status !! Status flag assigned to this outcome + integer(I4B) :: status !! Status flag assigned to this outcome ! Internals - integer(I4B) :: i, j, k, ibiggest, nfamily - real(DP) :: volume_new, pe - real(DP), dimension(NDIM) :: xc, vc, xcrossv - real(DP), dimension(2) :: vol - real(DP), dimension(NDIM) :: L_orb_old, L_spin_old + integer(I4B) :: i, j, k, ibiggest + real(DP), dimension(2) :: volume, density + real(DP) :: pe real(DP), dimension(NDIM) :: L_spin_new - logical, dimension(system%pl%nbody) :: lmask - real(DP), dimension(NDIM, 1) :: vb_frag, xb_frag, rot_frag, Ip_frag - real(DP), dimension(1) :: m_frag, rad_frag - integer(I4B), dimension(1) :: id_frag character(len=STRMAX) :: collider_message character(len=NAMELEN) :: idstr collider_message = "Merging" - call symba_collision_collider_message(system%pl, family, collider_message) + call symba_collision_collider_message(system%pl, colliders%idx, collider_message) write(*,*) trim(adjustl(collider_message)) select type(pl => system%pl) class is (symba_pl) - ibiggest = family(maxloc(pl%Gmass(family(:)), dim=1)) - id_frag(1) = pl%id(ibiggest) - - m_frag(1) = sum(mass(:)) - - ! Merged body is created at the barycenter of the original bodies - xb_frag(:,1) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / m_frag(1) - vb_frag(:,1) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / m_frag(1) + call frag%setup(1, param) + call frag%set_mass_dist(colliders) + ibiggest = colliders%idx(maxloc(pl%Gmass(colliders%idx(:)), dim=1)) + frag%id(1) = pl%id(ibiggest) + frag%xb(:,1) = frag%xbcom(:) + frag%vb(:,1) = frag%vbcom(:) - ! Get mass weighted mean of Ip and - vol(:) = 4._DP / 3._DP * PI * radius(:)**3 - volume_new = sum(vol(:)) - rad_frag(1) = (3 * volume_new / (4 * PI))**(1._DP / 3._DP) - - L_orb_old(:) = 0.0_DP - - ! Compute orbital angular momentum of pre-impact system - do i = 1, 2 - xc(:) = x(:, i) - xb_frag(:,1) - vc(:) = v(:, i) - vb_frag(:,1) - xcrossv(:) = xc(:) .cross. vc(:) - L_orb_old(:) = L_orb_old(:) + mass(i) * xcrossv(:) - end do - if (param%lrotation) then - Ip_frag(:,1) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / m_frag(1) - L_spin_old(:) = L_spin(:,1) + L_spin(:,2) - ! Conserve angular momentum by putting pre-impact orbital momentum into spin of the new body - L_spin_new(:) = L_orb_old(:) + L_spin_old(:) + L_spin_new(:) = colliders%L_orbit(:,1) + colliders%L_orbit(:,2) + colliders%L_spin(:,1) + colliders%L_spin(:,2) ! Assume prinicpal axis rotation on 3rd Ip axis - rot_frag(:,1) = L_spin_new(:) / (Ip_frag(3,1) * m_frag(1) * rad_frag(1)**2) + frag%rot(:,1) = L_spin_new(:) / (frag%Ip(3,1) * frag%mass(1) * frag%radius(1)**2) else ! If spin is not enabled, we will consider the lost pre-collision angular momentum as "escaped" and add it to our bookkeeping variable - param%Lescape(:) = param%Lescape(:) + L_orb_old(:) + param%Lescape(:) = param%Lescape(:) + colliders%L_orbit(:,1) + colliders%L_orbit(:,2) end if - ! Keep track of the component of potential energy due to the pre-impact family for book-keeping - nfamily = size(family(:)) + ! Keep track of the component of potential energy due to the pre-impact colliders%idx for book-keeping pe = 0.0_DP - do j = 1, nfamily - do i = j + 1, nfamily + do j = 1, colliders%ncoll + do i = j + 1, colliders%ncoll pe = pe - pl%Gmass(i) * pl%mass(j) / norm2(pl%xb(:, i) - pl%xb(:, j)) end do end do @@ -299,8 +190,8 @@ module function symba_collision_casemerge(system, param, family, x, v, mass, rad ! Update any encounter lists that have the removed bodies in them so that they instead point to the new do k = 1, system%plplenc_list%nenc - do j = 1, nfamily - i = family(j) + do j = 1, colliders%ncoll + i = colliders%idx(j) if (i == ibiggest) cycle if (system%plplenc_list%id1(k) == pl%id(i)) then system%plplenc_list%id1(k) = pl%id(ibiggest) @@ -316,7 +207,7 @@ module function symba_collision_casemerge(system, param, family, x, v, mass, rad status = MERGED - call symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, status) + call symba_collision_mergeaddsub(system, param, colliders, frag, status) end select @@ -324,80 +215,37 @@ module function symba_collision_casemerge(system, param, family, x, v, mass, rad end function symba_collision_casemerge - module function symba_collision_casesupercatastrophic(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) + module function symba_collision_casesupercatastrophic(system, param, colliders, frag) result(status) !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Create the fragments resulting from a supercatastrophic collision !! implicit none ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision - real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision - real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision - real(DP), dimension(:), intent(inout) :: mass_res !! The distribution of fragment mass obtained by the regime calculation - real(DP), intent(inout) :: Qloss !! Energy lost during collision + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object ! Result - integer(I4B) :: status !! Status flag assigned to this outcome + integer(I4B) :: status !! Status flag assigned to this outcom ! Internals - integer(I4B) :: i, j, nfrag, ibiggest, nfamily, nstart, nend - real(DP) :: mtot, avg_dens, min_frag_mass - real(DP), dimension(NDIM) :: xcom, vcom - real(DP), dimension(2) :: vol - real(DP), dimension(NDIM) :: Ip_new - real(DP), dimension(:, :), allocatable :: vb_frag, xb_frag, rot_frag, Ip_frag - real(DP), dimension(:), allocatable :: m_frag, rad_frag - integer(I4B), dimension(:), allocatable :: id_frag + integer(I4B) :: i, j, nfrag logical :: lfailure - logical, dimension(system%pl%nbody) :: lmask character(len=STRMAX) :: collider_message character(len=NAMELEN) :: idstr collider_message = "Supercatastrophic disruption between" - call symba_collision_collider_message(system%pl, family, collider_message) + call symba_collision_collider_message(system%pl, colliders%idx, collider_message) write(*,*) trim(adjustl(collider_message)) - ! Collisional fragments will be uniformly distributed around the pre-impact barycenter nfrag = NFRAG_SUPERCAT - allocate(m_frag(nfrag)) - allocate(rad_frag(nfrag)) - allocate(id_frag(nfrag)) - allocate(xb_frag(NDIM, nfrag)) - allocate(vb_frag(NDIM, nfrag)) - allocate(rot_frag(NDIM, nfrag)) - allocate(Ip_frag(NDIM, nfrag)) - - mtot = sum(mass(:)) - xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot - vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot - - ! Get mass weighted mean of Ip and average density - Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mtot - vol(:) = 4._DP / 3._DP * pi * radius(:)**3 - avg_dens = mtot / sum(vol(:)) - - ! If we are adding the first and largest fragment (lr), check to see if its mass is SMALLER than an equal distribution of - ! mass between all fragments. If so, we will just distribute the mass equally between the fragments - min_frag_mass = mtot / nfrag - if (mass_res(1) < min_frag_mass) then - m_frag(:) = min_frag_mass - else - m_frag(1) = mass_res(1) - m_frag(2:nfrag) = (mtot - mass_res(1)) / (nfrag - 1) - end if - ! Distribute any residual mass if there is any and set the radius - m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) - rad_frag(1:nfrag) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) - id_frag(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)] - param%maxid = id_frag(nfrag) - - do i = 1, nfrag - Ip_frag(:, i) = Ip_new(:) - end do + call frag%setup(nfrag, param) + call frag%set_mass_dist(colliders) + frag%id(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)] + param%maxid = frag%id(nfrag) - call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & - nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) + ! Generate the position and velocity distributions of the fragments + call frag%generate_fragments(colliders, system, param, lfailure) if (lfailure) then write(*,*) 'No fragment solution found, so treat as a pure hit-and-run' @@ -405,22 +253,22 @@ module function symba_collision_casesupercatastrophic(system, param, family, x, nfrag = 0 select type(pl => system%pl) class is (symba_pl) - pl%status(family(:)) = status - pl%ldiscard(family(:)) = .false. - pl%lcollision(family(:)) = .false. + pl%status(colliders%idx(:)) = status + pl%ldiscard(colliders%idx(:)) = .false. + pl%lcollision(colliders%idx(:)) = .false. end select else ! Populate the list of new bodies write(*,'("Generating ",I2.0," fragments")') nfrag status = SUPERCATASTROPHIC - call symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, status) + call symba_collision_mergeaddsub(system, param, colliders, frag, status) end if return end function symba_collision_casesupercatastrophic - subroutine symba_collision_collider_message(pl, family, collider_message) + subroutine symba_collision_collider_message(pl, collidx, collider_message) !! author: David A. Minton !! !! Prints a nicely formatted message about which bodies collided, including their names and ids. @@ -428,19 +276,19 @@ subroutine symba_collision_collider_message(pl, family, collider_message) implicit none ! Arguments class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object - integer(I4B), dimension(:), intent(in) :: family !! Index of collisional family members + integer(I4B), dimension(:), intent(in) :: collidx !! Index of collisional colliders%idx members character(*), intent(inout) :: collider_message !! The message to print to the screen. ! Internals integer(I4B) :: i, n character(len=STRMAX) :: idstr - n = size(family) + n = size(collidx) if (n == 0) return do i = 1, n if (i > 1) collider_message = trim(adjustl(collider_message)) // " and " - collider_message = " " // trim(adjustl(collider_message)) // " " // trim(adjustl(pl%info(family(i))%name)) - write(idstr, '(I10)') pl%id(family(i)) + collider_message = " " // trim(adjustl(collider_message)) // " " // trim(adjustl(pl%info(collidx(i))%name)) + write(idstr, '(I10)') pl%id(collidx(i)) collider_message = trim(adjustl(collider_message)) // " (" // trim(adjustl(idstr)) // ") " end do @@ -534,8 +382,8 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec self%x2(:,k) = pl%xh(:,j) + system%cb%xb(:) self%v2(:,k) = pl%vb(:,j) if (lcollision(k)) then - ! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional family - if (pl%lcollision(i) .or. pl%lcollision(j)) call pl%make_family([i,j]) + ! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional colliders%idx + if (pl%lcollision(i) .or. pl%lcollision(j)) call pl%make_colliders([i,j]) ! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step pl%lcollision([i, j]) = .true. @@ -619,10 +467,10 @@ pure elemental function symba_collision_check_one(xr, yr, zr, vxr, vyr, vzr, Gmt end function symba_collision_check_one - function symba_collision_consolidate_familes(pl, cb, param, idx_parent, family, x, v, mass, radius, L_spin, Ip) result(lflag) + function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, colliders) result(lflag) !! author: David A. Minton !! - !! Loops through the pl-pl collision list and groups families together by index. Outputs the indices of all family members, + !! Loops through the pl-pl collision list and groups families together by index. Outputs the indices of all colliders%idx members, !! and pairs of quantities (x and v vectors, mass, radius, L_spin, and Ip) that can be used to resolve the collisional outcome. implicit none ! Arguments @@ -630,22 +478,21 @@ function symba_collision_consolidate_familes(pl, cb, param, idx_parent, family, class(symba_cb), intent(inout) :: cb !! SyMBA central body object class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(2), intent(inout) :: idx_parent !! Index of the two bodies considered the "parents" of the collision - integer(I4B), dimension(:), allocatable, intent(out) :: family !! List of indices of all bodies inovlved in the collision - real(DP), dimension(NDIM,2), intent(out) :: x, v, L_spin, Ip !! Output values that represent a 2-body equivalent of a possibly 2+ body collision - real(DP), dimension(2), intent(out) :: mass, radius !! Output values that represent a 2-body equivalent of a possibly 2+ body collision + class(fraggle_colliders), intent(out) :: colliders ! Result - logical :: lflag !! Logical flag indicating whether a family was successfully created or not + logical :: lflag !! Logical flag indicating whether a colliders%idx was successfully created or not ! Internals - type family_array + type collidx_array integer(I4B), dimension(:), allocatable :: id integer(I4B), dimension(:), allocatable :: idx - end type family_array - type(family_array), dimension(2) :: parent_child_index_array + end type collidx_array + type(collidx_array), dimension(2) :: parent_child_index_array integer(I4B), dimension(2) :: nchild - integer(I4B) :: i, j, fam_size, idx_child + integer(I4B) :: i, j, ncolliders, idx_child real(DP), dimension(2) :: volume, density real(DP) :: mchild, mtot, volchild real(DP), dimension(NDIM) :: xc, vc, xcom, vcom, xchild, vchild, xcrossv + real(DP), dimension(NDIM,2) :: mxc, vcc nchild(:) = pl%kin(idx_parent(:))%nchild ! If all of these bodies share a parent, but this is still a unique collision, move the last child @@ -663,9 +510,9 @@ function symba_collision_consolidate_familes(pl, cb, param, idx_parent, family, pl%kin(idx_parent(2))%parent = idx_parent(1) end if - mass(:) = pl%mass(idx_parent(:)) ! Note: This is meant to mass, not G*mass, as the collisional regime determination uses mass values that will be converted to Si - radius(:) = pl%radius(idx_parent(:)) - volume(:) = (4.0_DP / 3.0_DP) * PI * radius(:)**3 + colliders%mass(:) = pl%mass(idx_parent(:)) ! Note: This is meant to mass, not G*mass, as the collisional regime determination uses mass values that will be converted to Si + colliders%radius(:) = pl%radius(idx_parent(:)) + volume(:) = (4.0_DP / 3.0_DP) * PI * colliders%radius(:)**3 ! Group together the ids and indexes of each collisional parent and its children do j = 1, 2 @@ -682,23 +529,24 @@ function symba_collision_consolidate_familes(pl, cb, param, idx_parent, family, end associate end do - ! Consolidate the groups of collsional parents with any children they may have into a single "family" index array - fam_size = 2 + sum(nchild(:)) - allocate(family(fam_size)) - family = [parent_child_index_array(1)%idx(:),parent_child_index_array(2)%idx(:)] - fam_size = count(pl%lcollision(family(:))) - family = pack(family(:), pl%lcollision(family(:))) - L_spin(:,:) = 0.0_DP - Ip(:,:) = 0.0_DP + ! Consolidate the groups of collsional parents with any children they may have into a single "colliders%idx" index array + ncolliders = 2 + sum(nchild(:)) + allocate(colliders%idx(ncolliders)) + colliders%idx = [parent_child_index_array(1)%idx(:),parent_child_index_array(2)%idx(:)] + + colliders%ncoll = count(pl%lcollision(colliders%idx(:))) + colliders%idx = pack(colliders%idx(:), pl%lcollision(colliders%idx(:))) + colliders%L_spin(:,:) = 0.0_DP + colliders%Ip(:,:) = 0.0_DP ! Find the barycenter of each body along with its children, if it has any do j = 1, 2 - x(:, j) = pl%xh(:, idx_parent(j)) + cb%xb(:) - v(:, j) = pl%vb(:, idx_parent(j)) + colliders%xb(:, j) = pl%xh(:, idx_parent(j)) + cb%xb(:) + colliders%vb(:, j) = pl%vb(:, idx_parent(j)) ! Assume principal axis rotation about axis corresponding to highest moment of inertia (3rd Ip) if (param%lrotation) then - Ip(:, j) = mass(j) * pl%Ip(:, idx_parent(j)) - L_spin(:, j) = Ip(3, j) * radius(j)**2 * pl%rot(:, idx_parent(j)) + colliders%Ip(:, j) = colliders%mass(j) * pl%Ip(:, idx_parent(j)) + colliders%L_spin(:, j) = colliders%Ip(3, j) * colliders%radius(j)**2 * pl%rot(:, idx_parent(j)) end if if (nchild(j) > 0) then @@ -713,39 +561,47 @@ function symba_collision_consolidate_familes(pl, cb, param, idx_parent, family, ! Get angular momentum of the child-parent pair and add that to the spin ! Add the child's spin if (param%lrotation) then - xcom(:) = (mass(j) * x(:,j) + mchild * xchild(:)) / (mass(j) + mchild) - vcom(:) = (mass(j) * v(:,j) + mchild * vchild(:)) / (mass(j) + mchild) - xc(:) = x(:, j) - xcom(:) - vc(:) = v(:, j) - vcom(:) + xcom(:) = (colliders%mass(j) * colliders%xb(:,j) + mchild * xchild(:)) / (colliders%mass(j) + mchild) + vcom(:) = (colliders%mass(j) * colliders%vb(:,j) + mchild * vchild(:)) / (colliders%mass(j) + mchild) + xc(:) = colliders%xb(:, j) - xcom(:) + vc(:) = colliders%vb(:, j) - vcom(:) xcrossv(:) = xc(:) .cross. vc(:) - L_spin(:, j) = L_spin(:, j) + mass(j) * xcrossv(:) + colliders%L_spin(:, j) = colliders%L_spin(:, j) + colliders%mass(j) * xcrossv(:) xc(:) = xchild(:) - xcom(:) vc(:) = vchild(:) - vcom(:) xcrossv(:) = xc(:) .cross. vc(:) - L_spin(:, j) = L_spin(:, j) + mchild * xcrossv(:) + colliders%L_spin(:, j) = colliders%L_spin(:, j) + mchild * xcrossv(:) - L_spin(:, j) = L_spin(:, j) + mchild * pl%Ip(3, idx_child) * pl%radius(idx_child)**2 * pl%rot(:, idx_child) - Ip(:, j) = Ip(:, j) + mchild * pl%Ip(:, idx_child) + colliders%L_spin(:, j) = colliders%L_spin(:, j) + mchild * pl%Ip(3, idx_child) * pl%radius(idx_child)**2 * pl%rot(:, idx_child) + colliders%Ip(:, j) = colliders%Ip(:, j) + mchild * pl%Ip(:, idx_child) end if ! Merge the child and parent - mass(j) = mass(j) + mchild - x(:, j) = xcom(:) - v(:, j) = vcom(:) + colliders%mass(j) = colliders%mass(j) + mchild + colliders%xb(:, j) = xcom(:) + colliders%vb(:, j) = vcom(:) end do end if - density(j) = mass(j) / volume(j) - radius(j) = ((3 * mass(j)) / (density(j) * 4 * pi))**(1.0_DP / 3.0_DP) - if (param%lrotation) Ip(:, j) = Ip(:, j) / mass(j) + density(j) = colliders%mass(j) / volume(j) + colliders%radius(j) = (3 * volume(j) / (4 * PI))**(1.0_DP / 3.0_DP) + if (param%lrotation) colliders%Ip(:, j) = colliders%Ip(:, j) / colliders%mass(j) + + xcom(:) = (colliders%mass(1) * colliders%xb(:, 1) + colliders%mass(2) * colliders%xb(:, 2)) / sum(colliders%mass(:)) + vcom(:) = (colliders%mass(1) * colliders%vb(:, 1) + colliders%mass(2) * colliders%vb(:, 2)) / sum(colliders%mass(:)) + mxc(:, 1) = colliders%mass(1) * (colliders%xb(:, 1) - xcom(:)) + mxc(:, 2) = colliders%mass(2) * (colliders%xb(:, 2) - xcom(:)) + vcc(:, 1) = colliders%vb(:, 1) - vcom(:) + vcc(:, 2) = colliders%vb(:, 2) - vcom(:) + colliders%L_orbit(:,:) = mxc(:,:) .cross. vcc(:,:) end do lflag = .true. - ! Destroy the kinship relationships for all members of this family - call pl%reset_kinship(family(:)) + ! Destroy the kinship relationships for all members of this colliders%idx + call pl%reset_kinship(colliders%idx(:)) return - end function symba_collision_consolidate_familes + end function symba_collision_consolidate_colliders module subroutine symba_collision_encounter_extract_collisions(self, system, param) @@ -817,10 +673,10 @@ module subroutine symba_collision_encounter_extract_collisions(self, system, par end subroutine symba_collision_encounter_extract_collisions - module subroutine symba_collision_make_family_pl(self, idx) + module subroutine symba_collision_make_colliders_pl(self, idx) !! author: Jennifer L.L. Pouplin, Carlisle A. wishard, and David A. Minton !! - !! When a single body is involved in more than one collision in a single step, it becomes part of a family. + !! When a single body is involved in more than one collision in a single step, it becomes part of a colliders%idx. !! The largest body involved in a multi-body collision is the "parent" and all bodies that collide with it are its "children," !! including those that collide with the children. !! @@ -875,26 +731,23 @@ module subroutine symba_collision_make_family_pl(self, idx) end associate return - end subroutine symba_collision_make_family_pl + end subroutine symba_collision_make_colliders_pl - subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, status) + subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) !! author: David A. Minton !! !! Fills the pl_discards and pl_adds with removed and added bodies !! implicit none ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision - integer(I4B), dimension(:), intent(in) :: id_frag !! List of fragment ids - real(DP), dimension(:), intent(in) :: m_frag, rad_frag !! Distribution of fragment mass and radii - real(DP), dimension(:,:), intent(in) :: Ip_frag !! Fragment rotational inertia vectors - real(DP), dimension(:,:), intent(in) :: xb_frag, vb_frag, rot_frag !! Fragment barycentric position, barycentric velocity, and rotation vectors - integer(I4B), intent(in) :: status !! Status flag to assign to adds + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object + integer(I4B), intent(in) :: status !! Status flag to assign to adds ! Internals - integer(I4B) :: i, ibiggest, ismallest, iother, nstart, nend, nfamily, nfrag + integer(I4B) :: i, ibiggest, ismallest, iother, nstart, nend, ncolliders, nfrag logical, dimension(system%pl%nbody) :: lmask class(symba_pl), allocatable :: plnew, plsub character(*), parameter :: FRAGFMT = '("Newbody",I0.7)' @@ -905,87 +758,87 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, select type(pl_discards => system%pl_discards) class is (symba_merger) associate(info => pl%info, pl_adds => system%pl_adds, cb => system%cb, npl => pl%nbody) - ! Add the family bodies to the subtraction list - nfamily = size(family(:)) - nfrag = size(m_frag(:)) + ! Add the colliders%idx bodies to the subtraction list + ncolliders = colliders%ncoll + nfrag = size(frag%mass(:)) ! Setup new bodies allocate(plnew, mold=pl) call plnew%setup(nfrag, param) - ibiggest = family(maxloc(pl%Gmass(family(:)), dim=1)) - ismallest = family(minloc(pl%Gmass(family(:)), dim=1)) + ibiggest = colliders%idx(maxloc(pl%Gmass(colliders%idx(:)), dim=1)) + ismallest = colliders%idx(minloc(pl%Gmass(colliders%idx(:)), dim=1)) ! Copy over identification, information, and physical properties of the new bodies from the fragment list - plnew%id(1:nfrag) = id_frag(1:nfrag) - plnew%xb(:, 1:nfrag) = xb_frag(:, 1:nfrag) - plnew%vb(:, 1:nfrag) = vb_frag(:, 1:nfrag) + plnew%id(1:nfrag) = frag%id(1:nfrag) + plnew%xb(:, 1:nfrag) = frag%xb(:, 1:nfrag) + plnew%vb(:, 1:nfrag) = frag%vb(:, 1:nfrag) call pl%vb2vh(cb) call pl%xh2xb(cb) do i = 1, nfrag - plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) - plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) + plnew%xh(:,i) = frag%xb(:, i) - cb%xb(:) + plnew%vh(:,i) = frag%vb(:, i) - cb%vb(:) end do - plnew%mass(1:nfrag) = m_frag(1:nfrag) - plnew%Gmass(1:nfrag) = param%GU * m_frag(1:nfrag) - plnew%radius(1:nfrag) = rad_frag(1:nfrag) - plnew%density(1:nfrag) = m_frag(1:nfrag) / rad_frag(1:nfrag) + plnew%mass(1:nfrag) = frag%mass(1:nfrag) + plnew%Gmass(1:nfrag) = param%GU * frag%mass(1:nfrag) + plnew%radius(1:nfrag) = frag%radius(1:nfrag) + plnew%density(1:nfrag) = frag%mass(1:nfrag) / frag%radius(1:nfrag) call plnew%set_rhill(cb) select case(status) case(DISRUPTION) plnew%status(1:nfrag) = NEW_PARTICLE do i = 1, nfrag - write(newname, FRAGFMT) id_frag(i) + write(newname, FRAGFMT) frag%id(i) call plnew%info(i)%set_value(origin_type="Disruption", origin_time=param%t, name=newname, origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i)) end do - do i = 1, nfamily - if (family(i) == ibiggest) then + do i = 1, ncolliders + if (colliders%idx(i) == ibiggest) then iother = ismallest else iother = ibiggest end if - call pl%info(family(i))%set_value(status="Disruption", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) + call pl%info(colliders%idx(i))%set_value(status="Disruption", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) end do case(SUPERCATASTROPHIC) plnew%status(1:nfrag) = NEW_PARTICLE do i = 1, nfrag - write(newname, FRAGFMT) id_frag(i) + write(newname, FRAGFMT) frag%id(i) call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=param%t, name=newname, origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i)) end do - do i = 1, nfamily - if (family(i) == ibiggest) then + do i = 1, ncolliders + if (colliders%idx(i) == ibiggest) then iother = ismallest else iother = ibiggest end if - call pl%info(family(i))%set_value(status="Supercatastrophic", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) + call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) end do case(HIT_AND_RUN_DISRUPT) call plnew%info(1)%copy(pl%info(ibiggest)) plnew%status(1) = OLD_PARTICLE do i = 2, nfrag - write(newname, FRAGFMT) id_frag(i) + write(newname, FRAGFMT) frag%id(i) call plnew%info(i)%set_value(origin_type="Hit and run fragmention", origin_time=param%t, name=newname, origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i)) end do - do i = 1, nfamily - if (family(i) == ibiggest) cycle + do i = 1, ncolliders + if (colliders%idx(i) == ibiggest) cycle iother = ibiggest - call pl%info(family(i))%set_value(status="Hit and run fragmention", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) + call pl%info(colliders%idx(i))%set_value(status="Hit and run fragmention", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) end do case(MERGED) call plnew%info(1)%copy(pl%info(ibiggest)) plnew%status(1) = OLD_PARTICLE - do i = 1, nfamily - if (family(i) == ibiggest) cycle + do i = 1, ncolliders + if (colliders%idx(i) == ibiggest) cycle iother = ibiggest - call pl%info(family(i))%set_value(status="MERGED", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) + call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) end do end select if (param%lrotation) then - plnew%Ip(:, 1:nfrag) = Ip_frag(:, 1:nfrag) - plnew%rot(:, 1:nfrag) = rot_frag(:, 1:nfrag) + plnew%Ip(:, 1:nfrag) = frag%Ip(:, 1:nfrag) + plnew%rot(:, 1:nfrag) = frag%rot(:, 1:nfrag) end if if (param%ltides) then @@ -1008,11 +861,11 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, pl_adds%ncomp(nstart:nend) = plnew%nbody ! Add the discarded bodies to the discard list - pl%status(family(:)) = MERGED - pl%ldiscard(family(:)) = .true. - pl%lcollision(family(:)) = .true. + pl%status(colliders%idx(:)) = MERGED + pl%ldiscard(colliders%idx(:)) = .true. + pl%lcollision(colliders%idx(:)) = .true. lmask(:) = .false. - lmask(family(:)) = .true. + lmask(colliders%idx(:)) = .true. call plnew%setup(0, param) deallocate(plnew) @@ -1021,11 +874,11 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, call pl%spill(plsub, lmask, ldestructive=.false.) nstart = pl_discards%nbody + 1 - nend = pl_discards%nbody + nfamily - call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, nfamily)]) + nend = pl_discards%nbody + ncolliders + call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, ncolliders)]) ! Record how many bodies were subtracted in this event - pl_discards%ncomp(nstart:nend) = nfamily + pl_discards%ncomp(nstart:nend) = ncolliders call plsub%setup(0, param) deallocate(plsub) @@ -1049,18 +902,13 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions ! Internals ! Internals - integer(I4B), dimension(:), allocatable :: family !! List of indices of all bodies inovlved in the collision integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision real(DP), dimension(NDIM,2) :: x, v, L_spin, Ip !! Output values that represent a 2-body equivalent of a possibly 2+ body collision real(DP), dimension(2) :: mass, radius !! Output values that represent a 2-body equivalent of a possibly 2+ body collision logical :: lgoodcollision - integer(I4B) :: i, jtarg, jproj, regime - real(DP), dimension(2) :: radius_si, mass_si, density_si - real(DP) :: min_mfrag_si, Mcb_si - real(DP), dimension(NDIM) :: x1_si, v1_si, x2_si, v2_si - real(DP) :: mlr, mslr, mtot, dentot, msys, msys_new, Qloss, impact_parameter - integer(I4B), parameter :: NRES = 3 !! Number of collisional product results - real(DP), dimension(NRES) :: mass_res + integer(I4B) :: i + type(fraggle_colliders) :: colliders !! Fraggle colliders object + type(fraggle_fragments) :: frag !! Fraggle fragmentation system object associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2) select type(pl => system%pl) @@ -1070,52 +918,20 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) do i = 1, ncollisions idx_parent(1) = pl%kin(idx1(i))%parent idx_parent(2) = pl%kin(idx2(i))%parent - lgoodcollision = symba_collision_consolidate_familes(pl, cb, param, idx_parent, family, x, v, mass, radius, L_spin, Ip) + lgoodcollision = symba_collision_consolidate_colliders(pl, cb, param, idx_parent, colliders) if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLISION)) cycle + + call colliders%regime(frag, system, param) - ! Convert all quantities to SI units and determine which of the pair is the projectile vs. target before sending them - ! to symba_regime - if (mass(1) > mass(2)) then - jtarg = 1 - jproj = 2 - else - jtarg = 2 - jproj = 1 - end if - mass_si(:) = (mass(:)) * param%MU2KG !! The collective mass of the parent and its children - radius_si(:) = radius(:) * param%DU2M !! The collective radius of the parent and its children - x1_si(:) = plplcollision_list%x1(:,i) * param%DU2M !! The position of the parent from inside the step (at collision) - v1_si(:) = plplcollision_list%v1(:,i) * param%DU2M / param%TU2S !! The velocity of the parent from inside the step (at collision) - x2_si(:) = plplcollision_list%x2(:,i) * param%DU2M !! The position of the parent from inside the step (at collision) - v2_si(:) = plplcollision_list%v2(:,i) * param%DU2M / param%TU2S !! The velocity of the parent from inside the step (at collision) - density_si(:) = mass_si(:) / (4.0_DP / 3._DP * PI * radius_si(:)**3) !! The collective density of the parent and its children - Mcb_si = cb%mass * param%MU2KG - min_mfrag_si = (param%min_GMfrag / param%GU) * param%MU2KG - - mass_res(:) = 0.0_DP - - mtot = sum(mass_si(:)) - dentot = sum(mass_si(:) * density_si(:)) / mtot - - !! Use the positions and velocities of the parents from indside the step (at collision) to calculate the collisional regime - call fragmentation_regime(Mcb_si, mass_si(jtarg), mass_si(jproj), radius_si(jtarg), radius_si(jproj), x1_si(:), x2_si(:),& - v1_si(:), v2_si(:), density_si(jtarg), density_si(jproj), regime, mlr, mslr, min_mfrag_si, Qloss) - - mass_res(1) = min(max(mlr, 0.0_DP), mtot) - mass_res(2) = min(max(mslr, 0.0_DP), mtot) - mass_res(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot) - mass_res(:) = (mass_res(:) / param%MU2KG) - Qloss = Qloss * (param%TU2S / param%DU2M)**2 / param%MU2KG - - select case (regime) + select case (frag%regime) case (COLLRESOLVE_REGIME_DISRUPTION) - plplcollision_list%status(i) = symba_collision_casedisruption(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) + plplcollision_list%status(i) = symba_collision_casedisruption(system, param, colliders, frag) case (COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - plplcollision_list%status(i) = symba_collision_casesupercatastrophic(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) + plplcollision_list%status(i) = symba_collision_casesupercatastrophic(system, param, colliders, frag) case (COLLRESOLVE_REGIME_HIT_AND_RUN) - plplcollision_list%status(i) = symba_collision_casehitandrun(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) + plplcollision_list%status(i) = symba_collision_casehitandrun(system, param, colliders, frag) case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - plplcollision_list%status(i) = symba_collision_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) + plplcollision_list%status(i) = symba_collision_casemerge(system, param, colliders, frag) case default write(*,*) "Error in symba_collision, unrecognized collision regime" call util_exit(FAILURE) @@ -1140,12 +956,13 @@ module subroutine symba_collision_resolve_mergers(self, system, param) class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions ! Internals - integer(I4B), dimension(:), allocatable :: family !! List of indices of all bodies inovlved in the collision integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision real(DP), dimension(NDIM,2) :: x, v, L_spin, Ip !! Output values that represent a 2-body equivalent of a possibly 2+ body collision real(DP), dimension(2) :: mass, radius !! Output values that represent a 2-body equivalent of a possibly 2+ body collision logical :: lgoodcollision integer(I4B) :: i + type(fraggle_colliders) :: colliders !! Fraggle colliders object + type(fraggle_fragments) :: frag !! Fraggle fragmentation system object associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2) select type(pl => system%pl) @@ -1155,11 +972,12 @@ module subroutine symba_collision_resolve_mergers(self, system, param) do i = 1, ncollisions idx_parent(1) = pl%kin(idx1(i))%parent idx_parent(2) = pl%kin(idx2(i))%parent - lgoodcollision = symba_collision_consolidate_familes(pl, cb, param, idx_parent, family, x, v, mass, radius, L_spin, Ip) + lgoodcollision = symba_collision_consolidate_colliders(pl, cb, param, idx_parent, colliders) if (.not. lgoodcollision) cycle if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved - - plplcollision_list%status(i) = symba_collision_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) + + frag%regime = COLLRESOLVE_REGIME_MERGE + plplcollision_list%status(i) = symba_collision_casemerge(system, param, colliders, frag) end do end select end select diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 4a40100c0..b68635b0e 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -289,8 +289,6 @@ module subroutine symba_util_index_eucl_plpl(self, param) associate(pl => self) npl = int(self%nbody, kind=I8B) - call pl%sort("mass", ascending=.false.) - nplm = count(.not. pl%lmtiny(1:npl)) pl%nplm = int(nplm, kind=I4B) @@ -476,6 +474,7 @@ module subroutine symba_util_rearray_pl(self, system, param) deallocate(ldump_mask) ! Reindex the new list of bodies + call pl%sort("mass", ascending=.false.) call pl%index(param) ! Reset the kinship trackers diff --git a/src/whm/whm_setup.f90 b/src/whm/whm_setup.f90 index 8dee15982..3c5cdeec3 100644 --- a/src/whm/whm_setup.f90 +++ b/src/whm/whm_setup.f90 @@ -81,6 +81,8 @@ module subroutine whm_setup_initialize_system(self, param) call setup_initialize_system(self, param) ! First we need to make sure that the massive bodies are sorted by heliocentric distance before computing jacobies call util_set_ir3h(self%pl) + call self%pl%sort("ir3h", ascending=.false.) + call self%pl%index(param) ! Make sure that the discard list gets allocated initially call self%tp_discards%setup(0, param) diff --git a/src/whm/whm_util.f90 b/src/whm/whm_util.f90 index bedf359ae..e8ed91854 100644 --- a/src/whm/whm_util.f90 +++ b/src/whm/whm_util.f90 @@ -67,22 +67,6 @@ module subroutine whm_util_fill_pl(self, inserts, lfill_list) end subroutine whm_util_fill_pl - module subroutine whm_util_index_eucl_plpl(self, param) - !! author: David A. Minton - !! - !! Wrapper for the indexing method for WHM massive bodies. Sorts the massive bodies by heliocentric distance and then flattens the pl-pl upper triangular matrix - implicit none - ! Arguments - class(whm_pl), intent(inout) :: self !! WHM massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - - call self%sort("ir3h", ascending=.false.) - call util_index_eucl_plpl(self, param) - - return - end subroutine whm_util_index_eucl_plpl - - module subroutine whm_util_resize_pl(self, nnew) !! author: David A. Minton !!