diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 54b431287..bd1df3227 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -38,6 +38,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes logical, dimension(size(IEEE_USUAL)) :: fpe_flag character(len=STRMAX) :: message + real(DP), dimension(NDIM) :: runit, vunit ! The minimization and linear solvers can sometimes lead to floating point exceptions. Rather than halting the code entirely if this occurs, we ! can simply fail the attempt and try again. So we need to turn off any floating point exception halting modes temporarily @@ -45,7 +46,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa fpe_quiet_modes(:) = .false. call ieee_set_halting_mode(IEEE_ALL,fpe_quiet_modes) - associate(frag => self, nfrag => self%nbody, pl => system%pl) + associate(fragments => self, nfrag => self%nbody, pl => system%pl) write(message,*) nfrag call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle generating " // trim(adjustl(message)) // " fragments.") @@ -55,7 +56,16 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa lfailure = .true. return end if - f_spin = F_SPIN_FIRST + + ! Get the unit vectors for the relative position and velocity vectors. These are used to shift the fragment cloud depending on the + runit(:) = colliders%rb(:,2) - colliders%rb(:,1) + runit(:) = runit(:) / (.mag. runit(:)) + + vunit(:) = colliders%vb(:,2) - colliders%vb(:,1) + vunit(:) = vunit(:) / (.mag. vunit(:)) + + ! This is a factor that will "distort" the shape of the frgment cloud in the direction of the impact velocity + f_spin= .mag. (runit(:) .cross. vunit(:)) if (param%lflatten_interactions) then lk_plpl = allocated(pl%k_plpl) @@ -64,73 +74,74 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa lk_plpl = .false. end if - call frag%set_natural_scale(colliders) + call fragments%set_natural_scale(colliders) - call frag%reset() + call fragments%reset() ! Calculate the initial energy of the system without the collisional family - call frag%get_energy_and_momentum(colliders, system, param, lbefore=.true.) + call fragments%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 = 2 * norm2(colliders%rb(:,2) - colliders%rb(:,1)) + r_max_start = 1.2_DP * (norm2(colliders%rb(:,2) - colliders%rb(:,1))) lfailure = .false. try = 1 do while (try < MAXTRY) write(message,*) try call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle try " // trim(adjustl(message))) if (lfailure) then - call frag%restructure(colliders, try, f_spin, r_max_start) - call frag%reset() + call fragments%restructure(colliders, try, f_spin, r_max_start) + call fragments%reset() try = try + 1 end if lfailure = .false. call ieee_set_flag(ieee_all, .false.) ! Set all fpe flags to quiet - call fraggle_generate_pos_vec(frag, colliders, r_max_start) - call frag%set_coordinate_system(colliders) + call fraggle_generate_pos_vec(fragments, colliders, r_max_start) + call fragments%set_coordinate_system(colliders) ! Initial velocity guess will be the barycentric velocity of the colliding system so that the budgets are based on the much smaller collisional-frame velocities do concurrent (i = 1:nfrag) - frag%vb(:, i) = frag%vbcom(:) + fragments%vb(:, i) = fragments%vbcom(:) end do - call frag%get_energy_and_momentum(colliders, system, param, lbefore=.false.) - call frag%set_budgets() + call fragments%get_energy_and_momentum(colliders, system, param, lbefore=.false.) + call fragments%set_budgets() - call fraggle_generate_spins(frag, f_spin, lfailure) + call fraggle_generate_spins(fragments, colliders, f_spin, lfailure) if (lfailure) then call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed to find spins") cycle end if - call fraggle_generate_tan_vel(frag, lfailure) + call fraggle_generate_tan_vel(fragments, colliders, lfailure) if (lfailure) then call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed to find tangential velocities") cycle end if - call fraggle_generate_rad_vel(frag, lfailure) + call fraggle_generate_rad_vel(fragments, colliders, lfailure) if (lfailure) then call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed to find radial velocities") cycle end if - 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(:)) + call fragments%get_energy_and_momentum(colliders, system, param, lbefore=.false.) + dEtot = fragments%Etot_after - fragments%Etot_before + dLmag = .mag. (fragments%Ltot_after(:) - fragments%Ltot_before(:)) + exit - lfailure = ((abs(dEtot + frag%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) + lfailure = ((abs(dEtot + fragments%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) if (lfailure) then - write(message, *) dEtot, abs(dEtot + frag%Qloss) / FRAGGLE_ETOL + write(message, *) dEtot, abs(dEtot + fragments%Qloss) / FRAGGLE_ETOL call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed due to high energy error: " // & trim(adjustl(message))) cycle end if - lfailure = ((abs(dLmag) / (.mag.frag%Ltot_before)) > FRAGGLE_LTOL) + lfailure = ((abs(dLmag) / (.mag.fragments%Ltot_before)) > FRAGGLE_LTOL) if (lfailure) then - write(message,*) dLmag / (.mag.frag%Ltot_before(:)) + write(message,*) dLmag / (.mag.fragments%Ltot_before(:)) call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed due to high angular momentum error: " // & trim(adjustl(message))) cycle @@ -142,8 +153,8 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa if (.not.lfailure) exit write(message,*) "Fraggle failed due to a floating point exception: ", fpe_flag call io_log_one_message(FRAGGLE_LOG_OUT, message) - end do + end do write(message,*) try if (lfailure) then call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle fragment generation failed after " // & @@ -153,7 +164,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa trim(adjustl(message)) // " tries") end if - call frag%set_original_scale(colliders) + call fragments%set_original_scale(colliders) ! Restore the big array if (lk_plpl) call pl%flatten(param) @@ -164,7 +175,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa end subroutine fraggle_generate_fragments - subroutine fraggle_generate_pos_vec(frag, colliders, r_max_start) + subroutine fraggle_generate_pos_vec(fragments, 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 @@ -172,72 +183,83 @@ subroutine fraggle_generate_pos_vec(frag, colliders, r_max_start) !! The initial positions do not conserve energy or momentum, so these need to be adjusted later. implicit none ! Arguments - class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragment system object + class(fraggle_fragments), intent(inout) :: fragments !! 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 - real(DP), dimension(NDIM) :: runit + real(DP) :: dis, rad, r_max, fdistort + real(DP), dimension(NDIM) :: runit, vunit logical, dimension(:), allocatable :: loverlap integer(I4B) :: i, j - logical :: lfixdir + logical :: lnoncat, lhitandrun - associate(nfrag => frag%nbody) + associate(nfrag => fragments%nbody) allocate(loverlap(nfrag)) + lnoncat = (fragments%regime /= COLLRESOLVE_REGIME_SUPERCATASTROPHIC) ! For non-catastrophic impacts, make the fragments act like ejecta and point away from the impact point + lhitandrun = (fragments%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) ! Disruptive hit and runs have their own fragment distribution + ! 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(:)) - lfixdir = (frag%regime /= COLLRESOLVE_REGIME_SUPERCATASTROPHIC) ! For this style of impact, make the fragments act like eject and point away from the impact point - runit(:) = colliders%rb(:,2) - colliders%rb(:,1) + ! Get the unit vectors for the relative position and velocity vectors. These are used to shift the fragment cloud depending on the + runit(:) = colliders%rb(:,2) - colliders%rb(:,1) runit(:) = runit(:) / (.mag. runit(:)) + vunit(:) = colliders%vb(:,2) - colliders%vb(:,1) + vunit(:) = vunit(:) / (.mag. vunit(:)) + + ! This is a factor that will "distort" the shape of the frgment cloud in the direction of the impact velocity + fdistort = .mag. (runit(:) .cross. vunit(:)) + ! 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)) + call random_number(fragments%r_coll(:,3:nfrag)) loverlap(:) = .true. do while (any(loverlap(3:nfrag))) - frag%x_coll(:, 1) = colliders%rb(:, 1) - frag%rbcom(:) - frag%x_coll(:, 2) = colliders%rb(:, 2) - frag%rbcom(:) + fragments%r_coll(:, 1) = colliders%rb(:, 1) - fragments%rbcom(:) + fragments%r_coll(:, 2) = colliders%rb(:, 2) - fragments%rbcom(:) 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 - frag%x_coll(:, i) = frag%x_coll(:, i) + (frag%rbimp(:) - frag%rbcom(:)) ! Shift the center of the fragment cloud to the impact point rather than the CoM - if (lfixdir .and. dot_product(frag%x_coll(:,i), runit(:)) < 0.0_DP) frag%x_coll(:, i) = -frag%x_coll(:, i) ! Make sure the fragment cloud points away from the impact point + call random_number(fragments%r_coll(:,i)) + fragments%r_coll(:,i) = 2 * (fragments%r_coll(:, i) - 0.5_DP) + fragments%r_coll(:, i) = fragments%r_coll(:,i) + fdistort * vunit(:) + fragments%r_coll(:, i) = r_max * fragments%r_coll(:, i) + fragments%r_coll(:, i) = fragments%r_coll(:, i) + (fragments%rbimp(:) - fragments%rbcom(:)) ! Shift the center of the fragment cloud to the impact point rather than the CoM + !if (lnoncat .and. dot_product(fragments%r_coll(:,i), runit(:)) < 0.0_DP) fragments%r_coll(:, i) = -fragments%r_coll(:, i) ! Make sure the fragment cloud points away from the impact point 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))) + dis = norm2(fragments%r_coll(:,j) - fragments%r_coll(:,i)) + loverlap(i) = loverlap(i) .or. (dis <= (fragments%radius(i) + fragments%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) + call fraggle_util_shift_vector_to_origin(fragments%mass, fragments%r_coll) + call fragments%set_coordinate_system(colliders) do concurrent(i = 1:nfrag) - frag%rb(:,i) = frag%x_coll(:,i) + frag%rbcom(:) + fragments%rb(:,i) = fragments%r_coll(:,i) + fragments%rbcom(:) end do - frag%rbcom(:) = 0.0_DP + fragments%rbcom(:) = 0.0_DP do concurrent(i = 1:nfrag) - frag%rbcom(:) = frag%rbcom(:) + frag%mass(i) * frag%rb(:,i) + fragments%rbcom(:) = fragments%rbcom(:) + fragments%mass(i) * fragments%rb(:,i) end do - frag%rbcom(:) = frag%rbcom(:) / frag%mtot + fragments%rbcom(:) = fragments%rbcom(:) / fragments%mtot end associate return end subroutine fraggle_generate_pos_vec - subroutine fraggle_generate_spins(frag, f_spin, lfailure) + subroutine fraggle_generate_spins(fragments, 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. @@ -245,51 +267,77 @@ subroutine fraggle_generate_spins(frag, f_spin, lfailure) !! A failure will trigger a restructuring of the fragments so we will try new values of the radial position distribution. implicit none ! Arguments - class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragment system object + class(fraggle_fragments), intent(inout) :: fragments !! Fraggle fragment system object + class(fraggle_colliders), intent(inout) :: 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 + real(DP), dimension(NDIM) :: L_remainder, rot_L, rot_ke, L + real(DP), dimension(NDIM,fragments%nbody) :: frot_rand ! The random rotation factor applied to fragments + real(DP), parameter :: frot_rand_mag = 1.50_DP ! The magnitude of the rotation variation to apply to the fragments integer(I4B) :: i character(len=STRMAX) :: message + real(DP) :: ke_remainder, ke - associate(nfrag => frag%nbody) + associate(nfrag => fragments%nbody) lfailure = .false. + L_remainder(:) = fragments%L_budget(:) + ke_remainder = fragments%ke_budget - ! Start the first two bodies with the same rotation as the original two impactors, then distribute the remaining angular momentum among the rest - L_remainder(:) = frag%L_budget(:) - frag%rot(:,:) = 0.0_DP + ! Add a fraction of the orbit angular momentum of the second body to the spin of the first body to kick things off + L(:) = colliders%L_spin(:,1) + f_spin * (colliders%L_orbit(:,2) + colliders%L_spin(:,2)) + fragments%rot(:,1) = L(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1)) + L_remainder(:) = L_remainder(:) - L(:) - frag%ke_spin = 0.0_DP + ! Partition the spin momentum of the second body into the spin of the fragments, with some random variation + L(:) = colliders%L_spin(:,2) / (nfrag - 1) + + call random_number(frot_rand(:,2:nfrag)) + frot_rand(:,2:nfrag) = 2 * (frot_rand(:,2:nfrag) - 0.5_DP) * frot_rand_mag + + do i = 2, nfrag + rot_L(:) = L(:) + frot_rand(:,i) * .mag.L(:) + fragments%rot(:,i) = rot_L(:) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,1)) + L_remainder(:) = L_remainder(:) - rot_L(:) + end do + + ! Make sure we didn't blow our kinetic energy budget + do i = 1, nfrag + ke_remainder = ke_remainder - 0.5_DP * fragments%mass(i) * fragments%Ip(3, i) * fragments%radius(i)**2 * norm2(fragments%rot(:, i)) + end do + + ! Distributed most of the remaining angular momentum amongst all the particles + fragments%ke_spin = 0.0_DP if (norm2(L_remainder(:)) > FRAGGLE_LTOL) then - do i = 1, nfrag + do i = nfrag, 1, -1 ! 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)) + rot_ke(:) = sqrt(2 * f_spin * ke_remainder / (i * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3, i))) * L_remainder(:) / norm2(L_remainder(:)) + rot_L(:) = f_spin * L_remainder(:) / (i * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3, i)) if (norm2(rot_ke) < norm2(rot_L)) then - frag%rot(:,i) = rot_ke(:) + fragments%rot(:,i) = fragments%rot(:,i) + rot_ke(:) else - frag%rot(:, i) = rot_L(:) + fragments%rot(:, i) = fragments%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)) + ke = 0.5_DP * fragments%mass(i) * fragments%Ip(3, i) * fragments%radius(i)**2 * norm2(fragments%rot(:, i)) + L(:) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3, i) * fragments%rot(:, i) + ke_remainder = ke_remainder - ke + L_remainder(:) = L_remainder(:) - L(:) + fragments%ke_spin = fragments%ke_spin + ke end do end if - frag%ke_spin = 0.5_DP * frag%ke_spin - lfailure = ((frag%ke_budget - frag%ke_spin - frag%ke_orbit) < 0.0_DP) + lfailure = ((fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit) < 0.0_DP) if (lfailure) then call io_log_one_message(FRAGGLE_LOG_OUT, " ") call io_log_one_message(FRAGGLE_LOG_OUT, "Spin failure diagnostics") - write(message, *) frag%ke_budget + write(message, *) fragments%ke_budget call io_log_one_message(FRAGGLE_LOG_OUT, "ke_budget : " // trim(adjustl(message))) - write(message, *) frag%ke_spin + write(message, *) fragments%ke_spin call io_log_one_message(FRAGGLE_LOG_OUT, "ke_spin : " // trim(adjustl(message))) - write(message, *) frag%ke_orbit + write(message, *) fragments%ke_orbit call io_log_one_message(FRAGGLE_LOG_OUT, "ke_orbit : " // trim(adjustl(message))) - write(message, *) frag%ke_budget - frag%ke_spin - frag%ke_orbit + write(message, *) fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit call io_log_one_message(FRAGGLE_LOG_OUT, "ke_remainder : " // trim(adjustl(message))) end if @@ -299,12 +347,12 @@ subroutine fraggle_generate_spins(frag, f_spin, lfailure) end subroutine fraggle_generate_spins - subroutine fraggle_generate_tan_vel(frag, lfailure) + subroutine fraggle_generate_tan_vel(fragments, 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. + !! our fragment kinetic energy (fragments%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. @@ -314,86 +362,86 @@ subroutine fraggle_generate_tan_vel(frag, lfailure) !! A failure will trigger a restructuring of the fragments so we will try new values of the radial position distribution. implicit none ! Arguments - class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragment system object + class(fraggle_fragments), intent(inout) :: fragments !! Fraggle fragment system object + class(fraggle_colliders), intent(inout) :: 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 - real(DP), parameter :: VNOISE_MAG = 1e-3_DP !! Magnitude of the noise to apply to initial conditions to help minimizer find a solution in case of failure + integer(I4B) :: i, try + real(DP), parameter :: TOL_MIN = 1e0_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-12_DP + real(DP), parameter :: VNOISE_MAG = 1e-1_DP !! Magnitude of the noise to apply to initial conditions to help minimizer find a solution in case of failure integer(I4B), parameter :: MAXLOOP = 10 + integer(I4B), parameter :: MAXTRY = 100 real(DP) :: tol - real(DP), dimension(:), allocatable :: v_t_initial - real(DP), dimension(frag%nbody) :: kefrag, vnoise + real(DP), dimension(:), allocatable :: v_t_initial, v_t_output + real(DP), dimension(fragments%nbody) :: kefrag, vnoise type(lambda_obj_err) :: objective_function - real(DP), dimension(NDIM) :: Li, L_remainder, L_frag_tot + real(DP), dimension(NDIM) :: L_frag_tot character(len=STRMAX) :: message + real(DP) :: ke_diff - associate(nfrag => frag%nbody) + associate(nfrag => fragments%nbody) lfailure = .false. - allocate(v_t_initial, mold=frag%v_t_mag) - v_t_initial(:) = 0.0_DP - frag%v_coll(:,:) = 0.0_DP + allocate(v_t_initial, mold=fragments%v_t_mag) + do try = 1, MAXTRY + v_t_initial(1) = dot_product(colliders%vb(:,1),fragments%v_t_unit(:,1)) + do i = 2, nfrag + v_t_initial(i) = dot_product(colliders%vb(:,2), fragments%v_t_unit(:,i)) + end do + fragments%v_t_mag(:) = v_t_initial + + ! 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) + call util_minimize_bfgs(objective_function, nfrag-6, v_t_initial(7:nfrag), tol, MAXLOOP, lfailure, v_t_output) + fragments%v_t_mag(7:nfrag) = v_t_output(:) + ! 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) = fragments%v_t_mag(7:nfrag) + if (.not.lfailure) exit + tol = tol * 2_DP ! Keep increasing the tolerance until we converge on a solution + call random_number(vnoise(1:nfrag)) ! Adding a bit of noise to the initial conditions helps it find a solution more often + vnoise(:) = 1.0_DP + VNOISE_MAG * (2 * vnoise(:) - 1._DP) + v_t_initial(:) = v_t_initial(:) * vnoise(:) + end do + fragments%v_t_mag(1:nfrag) = solve_fragment_tan_vel(v_t_mag_input=v_t_initial(7:nfrag), lfailure=lfailure) - ! 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 + ! Perform one final shift of the radial velocity vectors to align with the center of mass of the collisional system (the origin) + fragments%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%v_r_unit(:,1:nfrag), fragments%v_t_mag(1:nfrag), & + fragments%v_t_unit(:,1:nfrag), fragments%mass(1:nfrag), fragments%vbcom(:)) + do concurrent (i = 1:nfrag) + fragments%v_coll(:,i) = fragments%vb(:,i) - fragments%vbcom(:) + 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) + ! 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) = fragments%mass(i) * dot_product(fragments%vb(:, i), fragments%vb(:, i)) + end do + fragments%ke_orbit = 0.5_DP * sum(kefrag(:)) - 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 we are over the energy budget, flag this as a failure so we can try again + ke_diff = fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit + lfailure = ke_diff < 0.0_DP if (.not.lfailure) exit - tol = tol * 2_DP ! Keep increasing the tolerance until we converge on a solution - call random_number(vnoise(1:nfrag)) ! Adding a bit of noise to the initial conditions helps it find a solution more often - vnoise(:) = 1.0_DP + VNOISE_MAG * (2 * vnoise(:) - 1._DP) - v_t_initial(:) = v_t_initial(:) * vnoise(:) - 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(:) + fragments%r_coll(:,:) = fragments%r_coll(:,:) * 1.1_DP 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) if (lfailure) then call io_log_one_message(FRAGGLE_LOG_OUT, " ") call io_log_one_message(FRAGGLE_LOG_OUT, "Tangential velocity failure diagnostics") - call frag%get_ang_mtm() - L_frag_tot = frag%L_spin(:) + frag%L_orbit(:) - write(message, *) .mag.(frag%L_budget(:) - L_frag_tot(:)) / (.mag.frag%Ltot_before(:)) + call fragments%get_ang_mtm() + L_frag_tot = fragments%L_spin(:) + fragments%L_orbit(:) + write(message, *) .mag.(fragments%L_budget(:) - L_frag_tot(:)) / (.mag.fragments%Ltot_before(:)) call io_log_one_message(FRAGGLE_LOG_OUT, "|L_remainder| : " // trim(adjustl(message))) - write(message, *) frag%ke_budget + write(message, *) fragments%ke_budget call io_log_one_message(FRAGGLE_LOG_OUT, "ke_budget : " // trim(adjustl(message))) - write(message, *) frag%ke_spin + write(message, *) fragments%ke_spin call io_log_one_message(FRAGGLE_LOG_OUT, "ke_spin : " // trim(adjustl(message))) - write(message, *) frag%ke_orbit + write(message, *) fragments%ke_orbit call io_log_one_message(FRAGGLE_LOG_OUT, "ke_tangential : " // trim(adjustl(message))) - write(message, *) frag%ke_budget - frag%ke_spin - frag%ke_orbit + write(message, *) fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit call io_log_one_message(FRAGGLE_LOG_OUT, "ke_radial : " // trim(adjustl(message))) end if end associate @@ -418,7 +466,7 @@ function solve_fragment_tan_vel(lfailure, v_t_mag_input) result(v_t_mag_output) 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) + associate(nfrag => fragments%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 @@ -427,17 +475,17 @@ function solve_fragment_tan_vel(lfailure, v_t_mag_input) result(v_t_mag_output) 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)) + A(1:3, i) = fragments%mass(i) * fragments%v_t_unit(:, i) + A(4:6, i) = fragments%mass(i) * fragments%rmag(i) * (fragments%v_r_unit(:, i) .cross. fragments%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(:)) + vtmp(:) = v_t_mag_input(i - 6) * fragments%v_t_unit(:, i) + L_lin_others(:) = L_lin_others(:) + fragments%mass(i) * vtmp(:) + L(:) = fragments%mass(i) * (fragments%r_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(:) + b(4:6) = fragments%L_budget(:) - fragments%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(:) @@ -458,19 +506,19 @@ function tangential_objective_function(v_t_mag_input, lfailure) result(fval) 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), dimension(NDIM,fragments%nbody) :: v_shift + real(DP), dimension(fragments%nbody) :: v_t_new, kearr real(DP) :: keo - associate(nfrag => frag%nbody) + associate(nfrag => fragments%nbody) lfailure = .false. v_t_new(:) = solve_fragment_tan_vel(v_t_mag_input=v_t_mag_input(:), lfailure=lfailure) - v_shift(:,:) = fraggle_util_vmag_to_vb(frag%v_r_mag, frag%v_r_unit, v_t_new, frag%v_t_unit, frag%mass, frag%vbcom) + v_shift(:,:) = fraggle_util_vmag_to_vb(fragments%v_r_mag, fragments%v_r_unit, v_t_new, fragments%v_t_unit, fragments%mass, fragments%vbcom) kearr = 0.0_DP do concurrent(i = 1:nfrag) - kearr(i) = frag%mass(i) * dot_product(v_shift(:, i), v_shift(:, i)) + kearr(i) = fragments%mass(i) * dot_product(v_shift(:, i), v_shift(:, i)) end do keo = 0.5_DP * sum(kearr(:)) fval = keo @@ -483,14 +531,15 @@ end function tangential_objective_function end subroutine fraggle_generate_tan_vel - subroutine fraggle_generate_rad_vel(frag, lfailure) + subroutine fraggle_generate_rad_vel(fragments, 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_fragments), intent(inout) :: fragments !! Fraggle fragment system object + class(fraggle_colliders), intent(inout) :: 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 @@ -500,55 +549,61 @@ subroutine fraggle_generate_rad_vel(frag, lfailure) real(DP) :: ke_radial, tol integer(I4B) :: i real(DP), dimension(:), allocatable :: v_r_initial - real(DP), dimension(frag%nbody) :: vnoise + real(DP), dimension(fragments%nbody) :: vnoise type(lambda_obj) :: objective_function character(len=STRMAX) :: message - associate(nfrag => frag%nbody) + associate(nfrag => fragments%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 related to equipartition of kinetic energy with some noise + allocate(v_r_initial, source=fragments%v_r_mag) + ! Initialize radial velocity magnitudes with a random value that related to equipartition of kinetic energy with some noise and scaled with respect to the initial distance + v_r_initial(1) = dot_product(colliders%vb(:,1),fragments%v_r_unit(:,1)) + fragments%ke_orbit = 0.5_DP * fragments%mass(1) * v_r_initial(1)**2 + + ke_radial = fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit call random_number(vnoise(1:nfrag)) vnoise(:) = 1.0_DP + VNOISE_MAG * (2 * vnoise(:) - 1.0_DP) - v_r_initial(1:nfrag) = sqrt(abs(2 * ke_radial) / (frag%mass(1:nfrag) * nfrag)) * vnoise(1:nfrag) + v_r_initial(2:nfrag) = -sqrt(abs(2 * ke_radial) / (fragments%mass(1:nfrag) * nfrag)) * vnoise(1: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) + call util_minimize_bfgs(objective_function, nfrag, v_r_initial, tol, MAXLOOP, lfailure, fragments%v_r_mag) if (.not.lfailure) exit tol = tol * 2 ! Keep increasing the tolerance until we converge on a solution - v_r_initial(:) = frag%v_r_mag(:) + v_r_initial(:) = fragments%v_r_mag(:) call random_number(vnoise(1:nfrag)) ! Adding a bit of noise to the initial conditions helps it find a solution more often vnoise(:) = 1.0_DP + VNOISE_MAG * (2 * vnoise(:) - 1._DP) v_r_initial(:) = v_r_initial(:) * vnoise(:) 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(:)) + fragments%ke_orbit = 0.0_DP + fragments%ke_spin = 0.0_DP + fragments%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%v_r_unit(:,1:nfrag), & + fragments%v_t_mag(1:nfrag), fragments%v_t_unit(:,1:nfrag), fragments%mass(1:nfrag), fragments%vbcom(:)) 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)) + fragments%v_coll(:, i) = fragments%vb(:, i) - fragments%vbcom(:) + fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * norm2(fragments%vb(:, i)) + fragments%ke_spin = fragments%ke_spin + fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * norm2(fragments%rot(:,i)) end do - frag%ke_orbit = 0.5_DP * frag%ke_orbit + fragments%ke_orbit = 0.5_DP * fragments%ke_orbit + fragments%ke_spin = 0.5_DP * fragments%ke_spin - lfailure = abs((frag%ke_budget - (frag%ke_orbit + frag%ke_spin)) / frag%ke_budget) > FRAGGLE_ETOL + lfailure = abs((fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin)) / fragments%ke_budget) > FRAGGLE_ETOL if (lfailure) then call io_log_one_message(FRAGGLE_LOG_OUT, " ") call io_log_one_message(FRAGGLE_LOG_OUT, "Radial velocity failure diagnostics") - write(message, *) frag%ke_budget + write(message, *) fragments%ke_budget call io_log_one_message(FRAGGLE_LOG_OUT, "ke_budget : " // trim(adjustl(message))) - write(message, *) frag%ke_spin + write(message, *) fragments%ke_spin call io_log_one_message(FRAGGLE_LOG_OUT, "ke_spin : " // trim(adjustl(message))) - write(message, *) frag%ke_orbit + write(message, *) fragments%ke_orbit call io_log_one_message(FRAGGLE_LOG_OUT, "ke_orbit : " // trim(adjustl(message))) - write(message, *) frag%ke_budget - (frag%ke_orbit + frag%ke_spin) + write(message, *) fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin) call io_log_one_message(FRAGGLE_LOG_OUT, "ke_remainder : " // trim(adjustl(message))) end if @@ -569,20 +624,20 @@ function radial_objective_function(v_r_mag_input) result(fval) ! Internals integer(I4B) :: i real(DP), dimension(:,:), allocatable :: v_shift - real(DP), dimension(frag%nbody) :: kearr + real(DP), dimension(fragments%nbody) :: kearr real(DP) :: keo, ke_radial, rotmag2, vmag2 - 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) - !$omp do simd firstprivate(frag) - do i = 1,frag%nbody - rotmag2 = frag%rot(1,i)**2 + frag%rot(2,i)**2 + frag%rot(3,i)**2 + allocate(v_shift, mold=fragments%vb) + v_shift(:,:) = fraggle_util_vmag_to_vb(v_r_mag_input, fragments%v_r_unit, fragments%v_t_mag, fragments%v_t_unit, fragments%mass, fragments%vbcom) + !$omp do simd firstprivate(fragments) + do i = 1,fragments%nbody + rotmag2 = fragments%rot(1,i)**2 + fragments%rot(2,i)**2 + fragments%rot(3,i)**2 vmag2 = v_shift(1,i)**2 + v_shift(2,i)**2 + v_shift(3,i)**2 - kearr(i) = frag%mass(i) * (frag%Ip(3, i) * frag%radius(i)**2 * rotmag2 + vmag2) + kearr(i) = fragments%mass(i) * (fragments%Ip(3, i) * fragments%radius(i)**2 * rotmag2 + vmag2) end do !$omp end do simd - keo = 2 * frag%ke_budget - sum(kearr(:)) - ke_radial = frag%ke_budget - frag%ke_orbit - frag%ke_spin + keo = 2 * fragments%ke_budget - sum(kearr(:)) + ke_radial = fragments%ke_budget - fragments%ke_orbit - fragments%ke_spin ! The following ensures that fval = 0 is a local minimum, which is what the BFGS method is searching for fval = (keo / (2 * ke_radial))**2 diff --git a/src/fraggle/fraggle_io.f90 b/src/fraggle/fraggle_io.f90 index 61650e700..3253bcdb2 100644 --- a/src/fraggle/fraggle_io.f90 +++ b/src/fraggle/fraggle_io.f90 @@ -224,14 +224,14 @@ module subroutine fraggle_io_write_frame(self, nc, param) end subroutine fraggle_io_write_frame - module subroutine fraggle_io_log_regime(colliders, frag) + module subroutine fraggle_io_log_regime(colliders, fragments) !! author: David A. Minton !! !! Writes a log of the results of the collisional regime determination implicit none ! Arguments class(fraggle_colliders), intent(in) :: colliders !! Fraggle collider system object - class(fraggle_fragments), intent(in) :: frag !! Fraggle fragment object + class(fraggle_fragments), intent(in) :: fragments !! Fraggle fragment object ! Internals character(STRMAX) :: errmsg @@ -242,7 +242,7 @@ module subroutine fraggle_io_log_regime(colliders, frag) write(LUN, *) "--------------------------------------------------------------------" write(LUN, *) "True number of colliders : ",colliders%ncoll write(LUN, *) "Index list of true colliders : ",colliders%idx(1:colliders%ncoll) - select case(frag%regime) + select case(fragments%regime) case(COLLRESOLVE_REGIME_MERGE) write(LUN, *) "Regime: Merge" case(COLLRESOLVE_REGIME_DISRUPTION) @@ -254,7 +254,7 @@ module subroutine fraggle_io_log_regime(colliders, frag) case(COLLRESOLVE_REGIME_HIT_AND_RUN) write(LUN, *) "Regime: Hit and run" end select - write(LUN, *) "Energy loss : ", frag%Qloss + write(LUN, *) "Energy loss : ", fragments%Qloss write(LUN, *) "--------------------------------------------------------------------" close(LUN) diff --git a/src/fraggle/fraggle_regime.f90 b/src/fraggle/fraggle_regime.f90 index 7962e6c25..dc5b07f58 100644 --- a/src/fraggle/fraggle_regime.f90 +++ b/src/fraggle/fraggle_regime.f90 @@ -12,7 +12,7 @@ contains - module subroutine fraggle_regime_colliders(self, frag, system, param) + module subroutine fraggle_regime_colliders(self, fragments, 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. @@ -20,7 +20,7 @@ module subroutine fraggle_regime_colliders(self, frag, system, param) implicit none ! Arguments class(fraggle_colliders), intent(inout) :: self !! Fraggle colliders object - class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragment system object + class(fraggle_fragments), intent(inout) :: fragments !! 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 @@ -60,27 +60,27 @@ module subroutine fraggle_regime_colliders(self, frag, system, param) !! 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) + min_mfrag_si, fragments%regime, mlr, mslr, fragments%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) + fragments%mass_dist(1) = min(max(mlr, 0.0_DP), mtot) + fragments%mass_dist(2) = min(max(mslr, 0.0_DP), mtot) + fragments%mass_dist(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot) ! Find the center of mass of the collisional system - frag%mtot = sum(colliders%mass(:)) - frag%rbcom(:) = (colliders%mass(1) * colliders%rb(:,1) + colliders%mass(2) * colliders%rb(:,2)) / frag%mtot - frag%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / frag%mtot + fragments%mtot = sum(colliders%mass(:)) + fragments%rbcom(:) = (colliders%mass(1) * colliders%rb(:,1) + colliders%mass(2) * colliders%rb(:,2)) / fragments%mtot + fragments%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / fragments%mtot ! Find the point of impact between the two bodies runit(:) = colliders%rb(:,2) - colliders%rb(:,1) runit(:) = runit(:) / (.mag. runit(:)) - frag%rbimp(:) = colliders%rb(:,1) + colliders%radius(1) * runit(:) + fragments%rbimp(:) = colliders%rb(:,1) + colliders%radius(1) * runit(:) ! 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 + fragments%mass_dist(:) = (fragments%mass_dist(:) / param%MU2KG) + fragments%Qloss = fragments%Qloss * (param%TU2S / param%DU2M)**2 / param%MU2KG - call fraggle_io_log_regime(colliders, frag) + call fraggle_io_log_regime(colliders, fragments) end associate return diff --git a/src/fraggle/fraggle_set.f90 b/src/fraggle/fraggle_set.f90 index 45cf41a92..3fbe50d2f 100644 --- a/src/fraggle/fraggle_set.f90 +++ b/src/fraggle/fraggle_set.f90 @@ -22,13 +22,13 @@ module subroutine fraggle_set_budgets_fragments(self) real(DP) :: dEtot real(DP), dimension(NDIM) :: dL - associate(frag => self) + associate(fragments => self) - dEtot = frag%Etot_after - frag%Etot_before - dL(:) = frag%Ltot_after(:) - frag%Ltot_before(:) + dEtot = fragments%Etot_after - fragments%Etot_before + dL(:) = fragments%Ltot_after(:) - fragments%Ltot_before(:) - frag%L_budget(:) = -dL(:) - frag%ke_budget = -(dEtot - 0.5_DP * frag%mtot * dot_product(frag%vbcom(:), frag%vbcom(:))) - frag%Qloss + fragments%L_budget(:) = -dL(:) + fragments%ke_budget = -(dEtot - 0.5_DP * fragments%mtot * dot_product(fragments%vbcom(:), fragments%vbcom(:))) - fragments%Qloss end associate return @@ -59,10 +59,10 @@ module subroutine fraggle_set_mass_dist_fragments(self, colliders, param) integer(I4B), parameter :: iMslr = 2 integer(I4B), parameter :: iMrem = 3 - associate(frag => self) + associate(fragments => self) ! Get mass weighted mean of Ip and density volume(1:2) = 4._DP / 3._DP * PI * colliders%radius(1:2)**3 - Ip_avg(:) = (colliders%mass(1) * colliders%Ip(:,1) + colliders%mass(2) * colliders%Ip(:,2)) / frag%mtot + Ip_avg(:) = (colliders%mass(1) * colliders%Ip(:,1) + colliders%mass(2) * colliders%Ip(:,2)) / fragments%mtot if (colliders%mass(1) > colliders%mass(2)) then jtarg = 1 jproj = 2 @@ -71,7 +71,7 @@ module subroutine fraggle_set_mass_dist_fragments(self, colliders, param) jproj = 1 end if - select case(frag%regime) + select case(fragments%regime) case(COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC, COLLRESOLVE_REGIME_HIT_AND_RUN) ! The first two bins of the mass_dist are the largest and second-largest fragments that came out of fraggle_regime. ! The remainder from the third bin will be distributed among nfrag-2 bodies. The following code will determine nfrag based on @@ -84,7 +84,7 @@ module subroutine fraggle_set_mass_dist_fragments(self, colliders, param) ! The number of fragments we generate is bracked by the minimum required by fraggle_generate (7) and the ! maximum set by the NFRAG_SIZE_MULTIPLIER which limits the total number of fragments to prevent the nbody ! code from getting an overwhelmingly large number of fragments - nfrag = ceiling(NFRAG_SIZE_MULTIPLIER * log(frag%mtot / min_mfrag)) + nfrag = ceiling(NFRAG_SIZE_MULTIPLIER * log(fragments%mtot / min_mfrag)) nfrag = max(min(nfrag, NFRAGMAX), NFRAGMIN) class default min_mfrag = 0.0_DP @@ -92,36 +92,36 @@ module subroutine fraggle_set_mass_dist_fragments(self, colliders, param) end select i = iMrem - mremaining = frag%mass_dist(iMrem) + mremaining = fragments%mass_dist(iMrem) do while (i <= nfrag) - mfrag = (1 + i - iMslr)**(-3._DP / BETA) * frag%mass_dist(iMslr) + mfrag = (1 + i - iMslr)**(-3._DP / BETA) * fragments%mass_dist(iMslr) if (mremaining - mfrag < 0.0_DP) exit mremaining = mremaining - mfrag i = i + 1 end do if (i < nfrag) nfrag = max(i, NFRAGMIN) ! The sfd would actually give us fewer fragments than our maximum - call frag%setup(nfrag, param) + call fragments%setup(nfrag, param) case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - call frag%setup(1, param) - frag%mass(1) = frag%mass_dist(1) - frag%radius(1) = colliders%radius(jtarg) - frag%density(1) = frag%mass_dist(1) / volume(jtarg) - if (param%lrotation) frag%Ip(:, 1) = colliders%Ip(:,1) + call fragments%setup(1, param) + fragments%mass(1) = fragments%mass_dist(1) + fragments%radius(1) = colliders%radius(jtarg) + fragments%density(1) = fragments%mass_dist(1) / volume(jtarg) + if (param%lrotation) fragments%Ip(:, 1) = colliders%Ip(:,1) return case default - write(*,*) "fraggle_set_mass_dist_fragments error: Unrecognized regime code",frag%regime + write(*,*) "fraggle_set_mass_dist_fragments error: Unrecognized regime code",fragments%regime end select ! Make the first two bins the same as the Mlr and Mslr values that came from fraggle_regime - frag%mass(1) = frag%mass_dist(iMlr) - frag%mass(2) = frag%mass_dist(iMslr) + fragments%mass(1) = fragments%mass_dist(iMlr) + fragments%mass(2) = fragments%mass_dist(iMslr) ! Distribute the remaining mass the 3:nfrag bodies following the model SFD given by slope BETA - mremaining = frag%mass_dist(iMrem) + mremaining = fragments%mass_dist(iMrem) do i = iMrem, nfrag - mfrag = (1 + i - iMslr)**(-3._DP / BETA) * frag%mass_dist(iMslr) - frag%mass(i) = mfrag + mfrag = (1 + i - iMslr)**(-3._DP / BETA) * fragments%mass_dist(iMslr) + fragments%mass(i) = mfrag mremaining = mremaining - mfrag end do @@ -131,27 +131,27 @@ module subroutine fraggle_set_mass_dist_fragments(self, colliders, param) else ! If the remainder is postiive, this means that the number of fragments required by the SFD is larger than our upper limit set by computational expediency. istart = iMslr ! We will increase the mass of the 2:nfrag bodies to compensate, which ensures that the second largest fragment remains the second largest end if - mfrag = 1._DP + mremaining / sum(frag%mass(istart:nfrag)) - frag%mass(istart:nfrag) = frag%mass(istart:nfrag) * mfrag + mfrag = 1._DP + mremaining / sum(fragments%mass(istart:nfrag)) + fragments%mass(istart:nfrag) = fragments%mass(istart:nfrag) * mfrag ! There may still be some small residual due to round-off error. If so, simply add it to the last bin of the mass distribution. - mremaining = frag%mtot - sum(frag%mass(1:nfrag)) - frag%mass(nfrag) = frag%mass(nfrag) + mremaining + mremaining = fragments%mtot - sum(fragments%mass(1:nfrag)) + fragments%mass(nfrag) = fragments%mass(nfrag) + mremaining ! Compute physical properties of the new fragments - select case(frag%regime) + select case(fragments%regime) case(COLLRESOLVE_REGIME_HIT_AND_RUN) ! The hit and run case always preserves the largest body intact, so there is no need to recompute the physical properties of the first fragment - frag%radius(1) = colliders%radius(jtarg) - frag%density(1) = frag%mass_dist(iMlr) / volume(jtarg) - frag%Ip(:, 1) = colliders%Ip(:,1) + fragments%radius(1) = colliders%radius(jtarg) + fragments%density(1) = fragments%mass_dist(iMlr) / volume(jtarg) + fragments%Ip(:, 1) = colliders%Ip(:,1) istart = 2 case default istart = 1 end select - frag%density(istart:nfrag) = frag%mtot / sum(volume(:)) - frag%radius(istart:nfrag) = (3 * frag%mass(istart:nfrag) / (4 * PI * frag%density(istart:nfrag)))**(1.0_DP / 3.0_DP) + fragments%density(istart:nfrag) = fragments%mtot / sum(volume(:)) + fragments%radius(istart:nfrag) = (3 * fragments%mass(istart:nfrag) / (4 * PI * fragments%density(istart:nfrag)))**(1.0_DP / 3.0_DP) do i = istart, nfrag - frag%Ip(:, i) = Ip_avg(:) + fragments%Ip(:, i) = Ip_avg(:) end do end associate @@ -174,7 +174,7 @@ module subroutine fraggle_set_coordinate_system(self, colliders) real(DP) :: r_col_norm, v_col_norm, L_mag real(DP), dimension(NDIM, self%nbody) :: L_sigma - associate(frag => self, nfrag => self%nbody) + associate(fragments => self, nfrag => self%nbody) delta_v(:) = colliders%vb(:, 2) - colliders%vb(:, 1) v_col_norm = .mag. delta_v(:) delta_r(:) = colliders%rb(:, 2) - colliders%rb(:, 1) @@ -183,29 +183,29 @@ module subroutine fraggle_set_coordinate_system(self, colliders) ! 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 + fragments%y_coll_unit(:) = delta_r(:) / r_col_norm L_mag = .mag.Ltot(:) if (L_mag > sqrt(tiny(L_mag))) then - frag%z_coll_unit(:) = Ltot(:) / L_mag + fragments%z_coll_unit(:) = Ltot(:) / L_mag else - call random_number(frag%z_coll_unit(:)) - frag%z_coll_unit(:) = frag%z_coll_unit(:) / (.mag.frag%z_coll_unit(:)) + call random_number(fragments%z_coll_unit(:)) + fragments%z_coll_unit(:) = fragments%z_coll_unit(:) / (.mag.fragments%z_coll_unit(:)) end if ! 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(:) + fragments%x_coll_unit(:) = fragments%y_coll_unit(:) .cross. fragments%z_coll_unit(:) - if (.not.any(frag%x_coll(:,:) > 0.0_DP)) return - frag%rmag(:) = .mag. frag%x_coll(:,:) + if (.not.any(fragments%r_coll(:,:) > 0.0_DP)) return + fragments%rmag(:) = .mag. fragments%r_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)) + do concurrent(i = 1:nfrag, fragments%rmag(i) > 0.0_DP) + fragments%v_r_unit(:, i) = fragments%r_coll(:, i) / fragments%rmag(i) + fragments%v_n_unit(:, i) = fragments%z_coll_unit(:) + 2e-1_DP * (L_sigma(:,i) - 0.5_DP) + fragments%v_n_unit(:, i) = fragments%v_n_unit(:, i) / (.mag. fragments%v_n_unit(:, i)) + fragments%v_t_unit(:, i) = fragments%v_n_unit(:, i) .cross. fragments%v_r_unit(:, i) + fragments%v_t_unit(:, i) = fragments%v_t_unit(:, i) / (.mag. fragments%v_t_unit(:, i)) end do end associate @@ -226,35 +226,35 @@ module subroutine fraggle_set_natural_scale_factors(self, colliders) ! Internals integer(I4B) :: i - associate(frag => self) + associate(fragments => self) ! Set scale factors - frag%Escale = 0.5_DP * (colliders%mass(1) * dot_product(colliders%vb(:,1), colliders%vb(:,1)) & + fragments%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 + fragments%dscale = sum(colliders%radius(:)) + fragments%mscale = fragments%mtot + fragments%vscale = sqrt(fragments%Escale / fragments%mscale) + fragments%tscale = fragments%dscale / fragments%vscale + fragments%Lscale = fragments%mscale * fragments%dscale * fragments%vscale ! Scale all dimensioned quantities of colliders and fragments - frag%rbcom(:) = frag%rbcom(:) / frag%dscale - frag%vbcom(:) = frag%vbcom(:) / frag%vscale - frag%rbimp(:) = frag%rbimp(:) / frag%dscale - colliders%rb(:,:) = colliders%rb(:,:) / 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 - colliders%L_orbit(:,:) = colliders%L_orbit(:,:) / frag%Lscale + fragments%rbcom(:) = fragments%rbcom(:) / fragments%dscale + fragments%vbcom(:) = fragments%vbcom(:) / fragments%vscale + fragments%rbimp(:) = fragments%rbimp(:) / fragments%dscale + colliders%rb(:,:) = colliders%rb(:,:) / fragments%dscale + colliders%vb(:,:) = colliders%vb(:,:) / fragments%vscale + colliders%mass(:) = colliders%mass(:) / fragments%mscale + colliders%radius(:) = colliders%radius(:) / fragments%dscale + colliders%L_spin(:,:) = colliders%L_spin(:,:) / fragments%Lscale + colliders%L_orbit(:,:) = colliders%L_orbit(:,:) / fragments%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 + fragments%mtot = fragments%mtot / fragments%mscale + fragments%mass = fragments%mass / fragments%mscale + fragments%radius = fragments%radius / fragments%dscale + fragments%Qloss = fragments%Qloss / fragments%Escale end associate return @@ -277,58 +277,58 @@ module subroutine fraggle_set_original_scale_factors(self, colliders) 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) + associate(fragments => self) ! Restore scale factors - frag%rbcom(:) = frag%rbcom(:) * frag%dscale - frag%vbcom(:) = frag%vbcom(:) * frag%vscale - frag%rbimp(:) = frag%rbimp(:) * frag%dscale + fragments%rbcom(:) = fragments%rbcom(:) * fragments%dscale + fragments%vbcom(:) = fragments%vbcom(:) * fragments%vscale + fragments%rbimp(:) = fragments%rbimp(:) * fragments%dscale - colliders%mass = colliders%mass * frag%mscale - colliders%radius = colliders%radius * frag%dscale - colliders%rb = colliders%rb * frag%dscale - colliders%vb = colliders%vb * frag%vscale - colliders%L_spin = colliders%L_spin * frag%Lscale + colliders%mass = colliders%mass * fragments%mscale + colliders%radius = colliders%radius * fragments%dscale + colliders%rb = colliders%rb * fragments%dscale + colliders%vb = colliders%vb * fragments%vscale + colliders%L_spin = colliders%L_spin * fragments%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%rot = frag%rot / frag%tscale - frag%x_coll = frag%x_coll * frag%dscale - frag%v_coll = frag%v_coll * frag%vscale + fragments%mtot = fragments%mtot * fragments%mscale + fragments%mass = fragments%mass * fragments%mscale + fragments%radius = fragments%radius * fragments%dscale + fragments%rot = fragments%rot / fragments%tscale + fragments%r_coll = fragments%r_coll * fragments%dscale + fragments%v_coll = fragments%v_coll * fragments%vscale - do i = 1, frag%nbody - frag%rb(:, i) = frag%x_coll(:, i) + frag%rbcom(:) - frag%vb(:, i) = frag%v_coll(:, i) + frag%vbcom(:) + do i = 1, fragments%nbody + fragments%rb(:, i) = fragments%r_coll(:, i) + fragments%rbcom(:) + fragments%vb(:, i) = fragments%v_coll(:, i) + fragments%vbcom(:) end do - frag%Qloss = frag%Qloss * frag%Escale + fragments%Qloss = fragments%Qloss * fragments%Escale - frag%Lorbit_before(:) = frag%Lorbit_before * frag%Lscale - frag%Lspin_before(:) = frag%Lspin_before * frag%Lscale - frag%Ltot_before(:) = frag%Ltot_before * frag%Lscale - frag%ke_orbit_before = frag%ke_orbit_before * frag%Escale - frag%ke_spin_before = frag%ke_spin_before * frag%Escale - frag%pe_before = frag%pe_before * frag%Escale - frag%Etot_before = frag%Etot_before * frag%Escale + fragments%Lorbit_before(:) = fragments%Lorbit_before * fragments%Lscale + fragments%Lspin_before(:) = fragments%Lspin_before * fragments%Lscale + fragments%Ltot_before(:) = fragments%Ltot_before * fragments%Lscale + fragments%ke_orbit_before = fragments%ke_orbit_before * fragments%Escale + fragments%ke_spin_before = fragments%ke_spin_before * fragments%Escale + fragments%pe_before = fragments%pe_before * fragments%Escale + fragments%Etot_before = fragments%Etot_before * fragments%Escale - frag%Lorbit_after(:) = frag%Lorbit_after * frag%Lscale - frag%Lspin_after(:) = frag%Lspin_after * frag%Lscale - frag%Ltot_after(:) = frag%Ltot_after * frag%Lscale - frag%ke_orbit_after = frag%ke_orbit_after * frag%Escale - frag%ke_spin_after = frag%ke_spin_after * frag%Escale - frag%pe_after = frag%pe_after * frag%Escale - frag%Etot_after = frag%Etot_after * frag%Escale + fragments%Lorbit_after(:) = fragments%Lorbit_after * fragments%Lscale + fragments%Lspin_after(:) = fragments%Lspin_after * fragments%Lscale + fragments%Ltot_after(:) = fragments%Ltot_after * fragments%Lscale + fragments%ke_orbit_after = fragments%ke_orbit_after * fragments%Escale + fragments%ke_spin_after = fragments%ke_spin_after * fragments%Escale + fragments%pe_after = fragments%pe_after * fragments%Escale + fragments%Etot_after = fragments%Etot_after * fragments%Escale - 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 + fragments%mscale = 1.0_DP + fragments%dscale = 1.0_DP + fragments%vscale = 1.0_DP + fragments%tscale = 1.0_DP + fragments%Lscale = 1.0_DP + fragments%Escale = 1.0_DP end associate call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) diff --git a/src/fraggle/fraggle_setup.f90 b/src/fraggle/fraggle_setup.f90 index ab31af995..bd660dc98 100644 --- a/src/fraggle/fraggle_setup.f90 +++ b/src/fraggle/fraggle_setup.f90 @@ -22,7 +22,7 @@ module subroutine fraggle_setup_reset_fragments(self) self%rb(:,:) = 0.0_DP self%vb(:,:) = 0.0_DP self%rot(:,:) = 0.0_DP - self%x_coll(:,:) = 0.0_DP + self%r_coll(:,:) = 0.0_DP self%v_coll(:,:) = 0.0_DP self%v_r_unit(:,:) = 0.0_DP self%v_t_unit(:,:) = 0.0_DP @@ -55,7 +55,7 @@ module subroutine fraggle_setup_fragments(self, n, param) call setup_pl(self, n, param) if (n < 0) return - if (allocated(self%x_coll)) deallocate(self%x_coll) + if (allocated(self%r_coll)) deallocate(self%r_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) @@ -67,7 +67,7 @@ module subroutine fraggle_setup_fragments(self, n, param) if (n == 0) return - allocate(self%x_coll(NDIM,n)) + allocate(self%r_coll(NDIM,n)) allocate(self%v_coll(NDIM,n)) allocate(self%v_r_unit(NDIM,n)) allocate(self%v_t_unit(NDIM,n)) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 3ed18f32d..c4fffe4a6 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -11,13 +11,13 @@ use swiftest contains - module subroutine fraggle_util_add_fragments_to_system(frag, colliders, system, param) + module subroutine fraggle_util_add_fragments_to_system(fragments, 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_fragments), intent(in) :: fragments !! 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 @@ -25,24 +25,24 @@ module subroutine fraggle_util_add_fragments_to_system(frag, colliders, system, integer(I4B) :: i, npl_before, npl_after logical, dimension(:), allocatable :: lexclude - associate(nfrag => frag%nbody, pl => system%pl, cb => system%cb) + associate(nfrag => fragments%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) + pl%mass(npl_before+1:npl_after) = fragments%mass(1:nfrag) + pl%Gmass(npl_before+1:npl_after) = fragments%mass(1:nfrag) * param%GU + pl%radius(npl_before+1:npl_after) = fragments%radius(1:nfrag) do concurrent (i = 1:nfrag) - pl%rb(:,npl_before+i) = frag%rb(:,i) - pl%vb(:,npl_before+i) = frag%vb(:,i) - pl%rh(:,npl_before+i) = frag%rb(:,i) - cb%rb(:) - pl%vh(:,npl_before+i) = frag%vb(:,i) - cb%vb(:) + pl%rb(:,npl_before+i) = fragments%rb(:,i) + pl%vb(:,npl_before+i) = fragments%vb(:,i) + pl%rh(:,npl_before+i) = fragments%rb(:,i) - cb%rb(:) + pl%vh(:,npl_before+i) = fragments%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) + pl%Ip(:,npl_before+1:npl_after) = fragments%Ip(:,1:nfrag) + pl%rot(:,npl_before+1:npl_after) = fragments%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. @@ -69,13 +69,13 @@ module subroutine fraggle_util_ang_mtm(self) ! Internals integer(I4B) :: i - associate(frag => self, nfrag => self%nbody) - frag%L_orbit(:) = 0.0_DP - frag%L_spin(:) = 0.0_DP + associate(fragments => self, nfrag => self%nbody) + fragments%L_orbit(:) = 0.0_DP + fragments%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) + fragments%L_orbit(:) = fragments%L_orbit(:) + fragments%mass(i) * (fragments%r_coll(:, i) .cross. fragments%v_coll(:, i)) + fragments%L_spin(:) = fragments%L_spin(:) + fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:, i) * fragments%rot(:, i) end do end associate @@ -83,13 +83,13 @@ module subroutine fraggle_util_ang_mtm(self) end subroutine fraggle_util_ang_mtm - module subroutine fraggle_util_construct_temporary_system(frag, system, param, tmpsys, tmpparam) + module subroutine fraggle_util_construct_temporary_system(fragments, system, param, tmpsys, tmpparam) !! Author: David A. Minton !! !! Constructs a temporary internal system consisting of active bodies and additional fragments. This internal temporary system is used to calculate system energy with and without fragments implicit none ! Arguments - class(fraggle_fragments), intent(in) :: frag !! Fraggle fragment system object + class(fraggle_fragments), intent(in) :: fragments !! Fraggle fragment system object class(swiftest_nbody_system), intent(in) :: system !! Original swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters class(swiftest_nbody_system), allocatable, intent(out) :: tmpsys !! Output temporary swiftest nbody system object @@ -98,7 +98,7 @@ module subroutine fraggle_util_construct_temporary_system(frag, system, param, t logical, dimension(:), allocatable :: linclude integer(I4B) :: npl_tot - associate(nfrag => frag%nbody, pl => system%pl, npl => system%pl%nbody, cb => system%cb) + associate(nfrag => fragments%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) @@ -123,7 +123,7 @@ module subroutine fraggle_util_construct_temporary_system(frag, system, param, t 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) + call tmpsys%rescale(tmpparam, fragments%mscale, fragments%dscale, fragments%tscale) end associate @@ -154,7 +154,7 @@ module subroutine fraggle_util_final_fragments(self) ! Arguments type(fraggle_fragments), intent(inout) :: self !! Fraggle encountar storage object - if (allocated(self%x_coll)) deallocate(self%x_coll) + if (allocated(self%r_coll)) deallocate(self%r_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) @@ -201,7 +201,7 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para class(swiftest_parameters), allocatable, save :: tmpparam integer(I4B) :: npl_before, npl_after - associate(frag => self, nfrag => self%nbody, pl => system%pl, cb => system%cb) + associate(fragments => 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 @@ -212,7 +212,7 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para npl_after = npl_before + nfrag if (lbefore) then - call fraggle_util_construct_temporary_system(frag, system, param, tmpsys, tmpparam) + call fraggle_util_construct_temporary_system(fragments, system, param, tmpsys, tmpparam) ! Build the exluded body logical mask for the *before* case: Only the original bodies are used to compute energy and momentum tmpsys%pl%status(colliders%idx(1:colliders%ncoll)) = ACTIVE tmpsys%pl%status(npl_before+1:npl_after) = INACTIVE @@ -223,7 +223,7 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para call util_exit(FAILURE) end if ! Build the exluded body logical mask for the *after* case: Only the new bodies are used to compute energy and momentum - call fraggle_util_add_fragments_to_system(frag, colliders, tmpsys, tmpparam) + call fraggle_util_add_fragments_to_system(fragments, colliders, tmpsys, tmpparam) tmpsys%pl%status(colliders%idx(1:colliders%ncoll)) = INACTIVE tmpsys%pl%status(npl_before+1:npl_after) = ACTIVE end if @@ -235,21 +235,21 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para ! 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 + fragments%Lorbit_before(:) = tmpsys%Lorbit(:) + fragments%Lspin_before(:) = tmpsys%Lspin(:) + fragments%Ltot_before(:) = tmpsys%Ltot(:) + fragments%ke_orbit_before = tmpsys%ke_orbit + fragments%ke_spin_before = tmpsys%ke_spin + fragments%pe_before = tmpsys%pe + fragments%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 - (frag%pe_after - frag%pe_before) ! Gotta be careful with PE when number of bodies changes. + fragments%Lorbit_after(:) = tmpsys%Lorbit(:) + fragments%Lspin_after(:) = tmpsys%Lspin(:) + fragments%Ltot_after(:) = tmpsys%Ltot(:) + fragments%ke_orbit_after = tmpsys%ke_orbit + fragments%ke_spin_after = tmpsys%ke_spin + fragments%pe_after = tmpsys%pe + fragments%Etot_after = tmpsys%te - (fragments%pe_after - fragments%pe_before) ! Gotta be careful with PE when number of bodies changes. end if end associate @@ -308,16 +308,16 @@ module subroutine fraggle_util_restructure(self, colliders, try, f_spin, r_max_s 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) + associate(fragments => 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_tot_deficit = - (fragments%ke_budget - fragments%ke_orbit - fragments%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_tot_deficit = ke_tot_deficit - (fragments%ke_budget - fragments%ke_orbit - fragments%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) diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index ebd5212cf..235451379 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -10,7 +10,7 @@ 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 Fraggle: *Frag*ment *g*eneration that conserves angular momentum (*L*) and energy (*E*) + !! Definition of classes and methods specific to Fraggle: *Fragment* *g*eneration that conserves angular momentum (*L*) and energy (*E*) use swiftest_globals use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system, swiftest_cb, swiftest_pl, swiftest_storage, netcdf_parameters use encounter_classes, only : encounter_snapshot, encounter_io_parameters, encounter_storage @@ -58,7 +58,7 @@ module fraggle_classes 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 :: r_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 @@ -114,9 +114,9 @@ module fraggle_classes !! NetCDF dimension and variable names for the enounter save object type, extends(encounter_io_parameters) :: fraggle_io_parameters - integer(I4B) :: stage_dimid !! ID for the stage dimension - integer(I4B) :: stage_varid !! ID for the stage variable - character(NAMELEN) :: stage_dimname = "stage" !! name of the stage dimension (before/after) + integer(I4B) :: stage_dimid !! ID for the stage dimension + integer(I4B) :: stage_varid !! ID for the stage variable + character(NAMELEN) :: stage_dimname = "stage" !! name of the stage dimension (before/after) character(len=6), dimension(2) :: stage_coords = ["before", "after"] !! The stage coordinate labels character(NAMELEN) :: event_dimname = "collision" !! Name of collision event dimension @@ -156,28 +156,28 @@ end subroutine fraggle_generate_fragments module subroutine fraggle_io_initialize_output(self, param) implicit none - class(fraggle_io_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset - class(swiftest_parameters), intent(in) :: param + class(fraggle_io_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine fraggle_io_initialize_output module subroutine fraggle_io_write_frame(self, nc, param) implicit none - class(fraggle_snapshot), intent(in) :: self !! Swiftest encounter structure + class(fraggle_snapshot), intent(in) :: self !! Swiftest encounter structure class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine fraggle_io_write_frame - module subroutine fraggle_io_log_regime(colliders, frag) + module subroutine fraggle_io_log_regime(colliders, fragments) implicit none class(fraggle_colliders), intent(in) :: colliders - class(fraggle_fragments), intent(in) :: frag + class(fraggle_fragments), intent(in) :: fragments end subroutine fraggle_io_log_regime !> 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(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time @@ -205,17 +205,17 @@ module subroutine fraggle_placeholder_step(self, system, param, t, dt) real(DP), intent(in) :: dt !! Stepsiz end subroutine fraggle_placeholder_step - module subroutine fraggle_regime_colliders(self, frag, system, param) + module subroutine fraggle_regime_colliders(self, fragments, system, param) implicit none class(fraggle_colliders), intent(inout) :: self !! Fraggle colliders object - class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragment system object + class(fraggle_fragments), intent(inout) :: fragments !! 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 end subroutine fraggle_regime_colliders module subroutine fraggle_set_budgets_fragments(self) implicit none - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object end subroutine fraggle_set_budgets_fragments module subroutine fraggle_set_coordinate_system(self, colliders) @@ -255,10 +255,10 @@ module subroutine fraggle_setup_reset_fragments(self) class(fraggle_fragments), intent(inout) :: self end subroutine fraggle_setup_reset_fragments - module subroutine fraggle_util_add_fragments_to_system(frag, colliders, system, param) + module subroutine fraggle_util_add_fragments_to_system(fragments, 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_fragments), intent(in) :: fragments !! 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 @@ -269,10 +269,10 @@ module subroutine fraggle_util_ang_mtm(self) 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) + module subroutine fraggle_util_construct_temporary_system(fragments, 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(fraggle_fragments), intent(in) :: fragments !! Fraggle fragment system object class(swiftest_nbody_system), intent(in) :: system !! Original swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters class(swiftest_nbody_system), allocatable, intent(out) :: tmpsys !! Output temporary swiftest nbody system object diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 2e32f8c1d..a3bf66ad7 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -1525,7 +1525,7 @@ module subroutine util_index_map_storage(self) class(swiftest_storage(*)), intent(inout) :: self !! Swiftest storage object end subroutine util_index_map_storage - module function util_minimize_bfgs(f, N, x0, eps, maxloop, lerr) result(x1) + module subroutine util_minimize_bfgs(f, N, x0, eps, maxloop, lerr, x1) use lambda_function implicit none integer(I4B), intent(in) :: N @@ -1534,8 +1534,8 @@ module function util_minimize_bfgs(f, N, x0, eps, maxloop, lerr) result(x1) real(DP), intent(in) :: eps logical, intent(out) :: lerr integer(I4B), intent(in) :: maxloop - real(DP), dimension(:), allocatable :: x1 - end function util_minimize_bfgs + real(DP), dimension(:), allocatable, intent(out) :: x1 + end subroutine util_minimize_bfgs module subroutine util_peri_tp(self, system, param) implicit none diff --git a/src/util/util_minimize_bfgs.f90 b/src/util/util_minimize_bfgs.f90 index cd9e8f8bc..970a0ae45 100644 --- a/src/util/util_minimize_bfgs.f90 +++ b/src/util/util_minimize_bfgs.f90 @@ -10,7 +10,7 @@ submodule (swiftest_classes) s_util_minimize_bfgs use swiftest contains - module function util_minimize_bfgs(f, N, x0, eps, maxloop, lerr) result(x1) + module subroutine util_minimize_bfgs(f, N, x0, eps, maxloop, lerr, x1) !! author: David A. Minton !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !! This function implements the Broyden-Fletcher-Goldfarb-Shanno method to determine the minimum of a function of N variables. @@ -36,7 +36,7 @@ module function util_minimize_bfgs(f, N, x0, eps, maxloop, lerr) result(x1) integer(I4B), intent(in) :: maxloop logical, intent(out) :: lerr ! Result - real(DP), dimension(:), allocatable :: x1 + real(DP), dimension(:), intent(out), allocatable :: x1 ! Internals integer(I4B) :: i, j, k, l, conv real(DP), parameter :: graddelta = 1e-4_DP !! Delta x for gradient calculations @@ -227,6 +227,7 @@ function minimize1D(f, x0, S, N, eps, lerr) result(astar) call bracket(f, x0, S, N, gam, step, alo, ahi, lerr) if (lerr) then !write(*,*) "BFGS bracketing step failed!" + !write(*,*) "alo: ",alo, "ahi: ", ahi return end if if (abs(alo - ahi) < eps) then @@ -588,5 +589,5 @@ subroutine quadfit(f, x0, S, N, eps, lo, hi, lerr) return end subroutine quadfit - end function util_minimize_bfgs + end subroutine util_minimize_bfgs end submodule s_util_minimize_bfgs \ No newline at end of file