diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index ad6ef5c87..f66904599 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -25,10 +25,10 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa 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 :: lk_plpl logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes logical, dimension(size(IEEE_USUAL)) :: fpe_flag + character(len=STRMAX) :: message ! 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 @@ -39,7 +39,8 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa 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." + write(message,*) "Fraggle needs at least ",NFRAG_MIN," fragments, but only ",nfrag," were given." + call fraggle_io_log_one_message(message) lfailure = .true. return end if @@ -60,6 +61,8 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa lfailure = .false. try = 1 do while (try < MAXTRY) + write(message,*) try + call fraggle_io_log_one_message("Fraggle try " // trim(adjustl(message))) if (lfailure) then call frag%restructure(colliders, try, f_spin, r_max_start) call frag%reset() @@ -82,19 +85,19 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa call fraggle_generate_spins(frag, colliders, f_spin, lfailure) if (lfailure) then - write(*,*) "Fraggle failed to find spins" + call fraggle_io_log_one_message("Fraggle failed to find spins") cycle end if call fraggle_generate_tan_vel(frag, colliders, lfailure) if (lfailure) then - write(*,*) "Fraggle failed to find tangential velocities" + call fraggle_io_log_one_message("Fraggle failed to find tangential velocities") cycle end if call fraggle_generate_rad_vel(frag, colliders, lfailure) if (lfailure) then - write(*,*) "Fraggle failed to find radial velocities" + call fraggle_io_log_one_message("Fraggle failed to find radial velocities") cycle end if @@ -104,13 +107,15 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa lfailure = ((abs(dEtot + frag%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) if (lfailure) then - write(*,*) "Fraggle failed due to high energy error: ",dEtot, abs(dEtot + frag%Qloss) / FRAGGLE_ETOL + write(message, *) dEtot, abs(dEtot + frag%Qloss) / FRAGGLE_ETOL + call fraggle_io_log_one_message("Fraggle failed due to high energy error: " // trim(adjustl(message))) cycle end if lfailure = ((abs(dLmag) / (.mag.frag%Ltot_before)) > FRAGGLE_LTOL) if (lfailure) then - write(*,*) "Fraggle failed due to high angular momentum error: ", dLmag / (.mag.frag%Ltot_before(:)) + write(message,*) dLmag / (.mag.frag%Ltot_before(:)) + call fraggle_io_log_one_message("Fraggle failed due to high angular momentum error: " // trim(adjustl(message))) cycle end if @@ -118,24 +123,18 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa call ieee_get_flag(ieee_usual, fpe_flag) lfailure = any(fpe_flag) if (.not.lfailure) exit - write(*,*) "Fraggle failed due to a floating point exception: ", fpe_flag + write(message,*) "Fraggle failed due to a floating point exception: ", fpe_flag + call fraggle_io_log_one_message(message) end do - write(*, "(' -------------------------------------------------------------------------------------')") - write(*, "(' Fraggle final diagnostic')") - write(*, "(' -------------------------------------------------------------------------------------')") + write(message,*) try if (lfailure) then - write(*,*) "Fraggle fragment generation failed after: ",try," tries" + call fraggle_io_log_one_message("Fraggle fragment generation failed after " // trim(adjustl(message)) // " 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) + call fraggle_io_log_one_message("Fraggle fragment generation succeeded after " // trim(adjustl(message)) // " tries") + call fraggle_io_log_generate(frag) end if - write(*, "(' -------------------------------------------------------------------------------------')") + call frag%set_original_scale(colliders) ! Restore the big array @@ -226,7 +225,8 @@ subroutine fraggle_generate_spins(frag, colliders, f_spin, lfailure) 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 + integer(I4B) :: i + character(len=STRMAX) :: message associate(nfrag => frag%nbody) lfailure = .false. @@ -252,6 +252,20 @@ subroutine fraggle_generate_spins(frag, colliders, f_spin, lfailure) frag%ke_spin = 0.5_DP * frag%ke_spin lfailure = ((frag%ke_budget - frag%ke_spin - frag%ke_orbit) < 0.0_DP) + + if (lfailure) then + call fraggle_io_log_one_message(" ") + call fraggle_io_log_one_message("Spin failure diagnostics") + write(message, *) frag%ke_budget + call fraggle_io_log_one_message("ke_budget : " // trim(adjustl(message))) + write(message, *) frag%ke_spin + call fraggle_io_log_one_message("ke_spin : " // trim(adjustl(message))) + write(message, *) frag%ke_orbit + call fraggle_io_log_one_message("ke_orbit : " // trim(adjustl(message))) + write(message, *) frag%ke_budget - frag%ke_spin - frag%ke_orbit + call fraggle_io_log_one_message("ke_remainder : " // trim(adjustl(message))) + end if + end associate return @@ -286,27 +300,12 @@ subroutine fraggle_generate_tan_vel(frag, colliders, lfailure) real(DP), dimension(frag%nbody) :: kefrag type(lambda_obj) :: spinfunc type(lambda_obj_err) :: objective_function - real(DP), dimension(NDIM) :: Li, L_remainder + real(DP), dimension(NDIM) :: Li, L_remainder, L_frag_tot + character(len=STRMAX) :: message 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 @@ -351,14 +350,22 @@ subroutine fraggle_generate_tan_vel(frag, colliders, lfailure) ! 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 + if (lfailure) then + call fraggle_io_log_one_message(" ") + call fraggle_io_log_one_message("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 fraggle_io_log_one_message("|L_remainder| : " // trim(adjustl(message))) + write(message, *) frag%ke_budget + call fraggle_io_log_one_message("ke_budget : " // trim(adjustl(message))) + write(message, *) frag%ke_spin + call fraggle_io_log_one_message("ke_spin : " // trim(adjustl(message))) + write(message, *) frag%ke_orbit + call fraggle_io_log_one_message("ke_tangential : " // trim(adjustl(message))) + write(message, *) frag%ke_budget - frag%ke_spin - frag%ke_orbit + call fraggle_io_log_one_message("ke_radial : " // trim(adjustl(message))) + end if end associate return @@ -465,6 +472,7 @@ subroutine fraggle_generate_rad_vel(frag, colliders, lfailure) real(DP), dimension(:), allocatable :: v_r_initial, v_r_sigma real(DP), dimension(:,:), allocatable :: v_r type(lambda_obj) :: objective_function + character(len=STRMAX) :: message associate(nfrag => frag%nbody) ! Set the "target" ke for the radial component @@ -497,13 +505,19 @@ subroutine fraggle_generate_rad_vel(frag, colliders, lfailure) 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. + lfailure = abs((frag%ke_budget - (frag%ke_orbit + frag%ke_spin)) / frag%ke_budget) > FRAGGLE_ETOL + if (lfailure) then + call fraggle_io_log_one_message(" ") + call fraggle_io_log_one_message("Radial velocity failure diagnostics") + write(message, *) frag%ke_budget + call fraggle_io_log_one_message("ke_budget : " // trim(adjustl(message))) + write(message, *) frag%ke_spin + call fraggle_io_log_one_message("ke_spin : " // trim(adjustl(message))) + write(message, *) frag%ke_orbit + call fraggle_io_log_one_message("ke_orbit : " // trim(adjustl(message))) + write(message, *) frag%ke_budget - (frag%ke_orbit + frag%ke_spin) + call fraggle_io_log_one_message("ke_remainder : " // trim(adjustl(message))) + end if end associate return diff --git a/src/fraggle/fraggle_io.f90 b/src/fraggle/fraggle_io.f90 index dbe412167..72d4adfa1 100644 --- a/src/fraggle/fraggle_io.f90 +++ b/src/fraggle/fraggle_io.f90 @@ -3,20 +3,178 @@ contains - module subroutine fraggle_io_log_regime(param, colliders, frag) + module subroutine fraggle_io_log_generate(frag) !! author: David A. Minton !! - !! Writes a log of the results of the collisional regime determination + !! Writes a log of the results of the fragment generation implicit none ! Arguments - class(swiftest_parameters), intent(in) :: param - class(fraggle_colliders), intent(in) :: colliders class(fraggle_fragments), intent(in) :: frag ! Internals + integer(I4B) :: i + character(STRMAX) :: errmsg + character(len=*), parameter :: fmtlabel = "(A14,10(ES11.4,1X,:))" + + open(unit=FRAGGLE_LOG_UNIT, file=FRAGGLE_LOG_OUT, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) + write(FRAGGLE_LOG_UNIT, *, err = 667, iomsg = errmsg) + write(FRAGGLE_LOG_UNIT, *) "--------------------------------------------------------------------" + write(FRAGGLE_LOG_UNIT, *) " Fraggle fragment generation results" + write(FRAGGLE_LOG_UNIT, *) "--------------------------------------------------------------------" + write(FRAGGLE_LOG_UNIT, "(' dL_tot should be very small' )") + write(FRAGGLE_LOG_UNIT,fmtlabel) ' dL_tot |', (.mag.(frag%Ltot_after(:) - frag%Ltot_before(:))) / (.mag.frag%Ltot_before(:)) + write(FRAGGLE_LOG_UNIT, "(' dE_tot should be negative and equal to Qloss' )") + write(FRAGGLE_LOG_UNIT,fmtlabel) ' dE_tot |', (frag%Etot_after - frag%Etot_before) / abs(frag%Etot_before) + write(FRAGGLE_LOG_UNIT,fmtlabel) ' Qloss |', -frag%Qloss / abs(frag%Etot_before) + write(FRAGGLE_LOG_UNIT,fmtlabel) ' dE - Qloss |', (frag%Etot_after - frag%Etot_before + frag%Qloss) / abs(frag%Etot_before) + write(FRAGGLE_LOG_UNIT, "(' -------------------------------------------------------------------------------------')") + write(FRAGGLE_LOG_UNIT, *) "Individual fragment values (collisional system natural units)" + write(FRAGGLE_LOG_UNIT, *) "mass" + do i = 1, frag%nbody + write(FRAGGLE_LOG_UNIT, *) i, frag%mass(i) + end do + write(FRAGGLE_LOG_UNIT, *) "x_coll" + do i = 1, frag%nbody + write(FRAGGLE_LOG_UNIT, *) i, frag%x_coll(:,i) + end do + write(FRAGGLE_LOG_UNIT, *) "v_coll" + do i = 1, frag%nbody + write(FRAGGLE_LOG_UNIT, *) i, frag%v_coll(:,i) + end do + write(FRAGGLE_LOG_UNIT, *) "xb" + do i = 1, frag%nbody + write(FRAGGLE_LOG_UNIT, *) i, frag%xb(:,i) + end do + write(FRAGGLE_LOG_UNIT, *) "vb" + do i = 1, frag%nbody + write(FRAGGLE_LOG_UNIT, *) i, frag%vb(:,i) + end do + write(FRAGGLE_LOG_UNIT, *) "rot" + do i = 1, frag%nbody + write(FRAGGLE_LOG_UNIT, *) i, frag%rot(:,i) + end do + + close(FRAGGLE_LOG_UNIT) + + return + 667 continue + write(*,*) "Error writing Fraggle message to log file: " // trim(adjustl(errmsg)) + end subroutine fraggle_io_log_generate + + + module subroutine fraggle_io_log_one_message(message) + !! author: David A. Minton + !! + !! Writes a single message to the fraggle log file + implicit none + ! Arguments + character(len=*), intent(in) :: message + ! Internals + character(STRMAX) :: errmsg + + open(unit=FRAGGLE_LOG_UNIT, file=FRAGGLE_LOG_OUT, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) + write(FRAGGLE_LOG_UNIT, *) trim(adjustl(message)) + close(FRAGGLE_LOG_UNIT) + + return + 667 continue + write(*,*) "Error writing Fraggle message to log file: " // trim(adjustl(errmsg)) + end subroutine fraggle_io_log_one_message + + + module subroutine fraggle_io_log_pl(pl, param) + !! author: David A. Minton + !! + !! Writes a single message to the fraggle log file + implicit none + ! Arguments + class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object (only the new bodies generated in a collision) + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + ! Internals + integer(I4B) :: i + character(STRMAX) :: errmsg + + open(unit=FRAGGLE_LOG_UNIT, file=FRAGGLE_LOG_OUT, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) + write(FRAGGLE_LOG_UNIT, *, err = 667, iomsg = errmsg) + + write(FRAGGLE_LOG_UNIT, *) "--------------------------------------------------------------------" + write(FRAGGLE_LOG_UNIT, *) " Fraggle fragment final body properties" + write(FRAGGLE_LOG_UNIT, *) "--------------------------------------------------------------------" + write(FRAGGLE_LOG_UNIT, *) "id, name" + do i = 1, pl%nbody + write(FRAGGLE_LOG_UNIT, *) i, pl%id(i), pl%info(i)%name + end do + write(FRAGGLE_LOG_UNIT, *) "mass, Gmass" + do i = 1, pl%nbody + write(FRAGGLE_LOG_UNIT, *) i, pl%mass(i), pl%Gmass(i) + end do + write(FRAGGLE_LOG_UNIT, *) "radius" + do i = 1, pl%nbody + write(FRAGGLE_LOG_UNIT, *) i, pl%radius(i) + end do + write(FRAGGLE_LOG_UNIT, *) "xb" + do i = 1, pl%nbody + write(FRAGGLE_LOG_UNIT, *) i, pl%xb(:,i) + end do + write(FRAGGLE_LOG_UNIT, *) "vb" + do i = 1, pl%nbody + write(FRAGGLE_LOG_UNIT, *) i, pl%vb(:,i) + end do + write(FRAGGLE_LOG_UNIT, *) "xh" + do i = 1, pl%nbody + write(FRAGGLE_LOG_UNIT, *) i, pl%xh(:,i) + end do + write(FRAGGLE_LOG_UNIT, *) "vh" + do i = 1, pl%nbody + write(FRAGGLE_LOG_UNIT, *) i, pl%vh(:,i) + end do + + if (param%lrotation) then + write(FRAGGLE_LOG_UNIT, *) "rot" + do i = 1, pl%nbody + write(FRAGGLE_LOG_UNIT, *) i, pl%rot(:,i) + end do + write(FRAGGLE_LOG_UNIT, *) "Ip" + do i = 1, pl%nbody + write(FRAGGLE_LOG_UNIT, *) i, pl%Ip(:,i) + end do + end if + + if (param%ltides) then + write(FRAGGLE_LOG_UNIT, *) "Q" + do i = 1, pl%nbody + write(FRAGGLE_LOG_UNIT, *) i, pl%Q(i) + end do + write(FRAGGLE_LOG_UNIT, *) "k2" + do i = 1, pl%nbody + write(FRAGGLE_LOG_UNIT, *) i, pl%k2(i) + end do + write(FRAGGLE_LOG_UNIT, *) "tlag" + do i = 1, pl%nbody + write(FRAGGLE_LOG_UNIT, *) i, pl%tlag(i) + end do + end if + + close(FRAGGLE_LOG_UNIT) + + return + 667 continue + write(*,*) "Error writing Fraggle message to log file: " // trim(adjustl(errmsg)) + end subroutine fraggle_io_log_pl + + + module subroutine fraggle_io_log_regime(colliders, frag) + !! 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 + ! Internals character(STRMAX) :: errmsg open(unit=FRAGGLE_LOG_UNIT, file=FRAGGLE_LOG_OUT, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) - write(FRAGGLE_LOG_UNIT, *) + write(FRAGGLE_LOG_UNIT, *, err = 667, iomsg = errmsg) write(FRAGGLE_LOG_UNIT, *) "--------------------------------------------------------------------" write(FRAGGLE_LOG_UNIT, *) " Fraggle collisional regime determination results" write(FRAGGLE_LOG_UNIT, *) "--------------------------------------------------------------------" @@ -61,8 +219,9 @@ module subroutine fraggle_io_log_regime(param, colliders, frag) write(FRAGGLE_LOG_UNIT, *) "Remaining fragment mass : ", frag%mass_dist(3) write(FRAGGLE_LOG_UNIT, *) "Center of mass position : ", frag%xbcom(:) write(FRAGGLE_LOG_UNIT, *) "Center of mass velocity : ", frag%vbcom(:) - write(FRAGGLE_LOG_UNIT, *) "Energy loss : ",frag%Qloss + write(FRAGGLE_LOG_UNIT, *) "Energy loss : ", frag%Qloss write(FRAGGLE_LOG_UNIT, *) "--------------------------------------------------------------------" + close(FRAGGLE_LOG_UNIT) return 667 continue @@ -86,6 +245,7 @@ module subroutine fraggle_io_log_start(param) open(unit=FRAGGLE_LOG_UNIT, file=FRAGGLE_LOG_OUT, status="REPLACE", err = 667, iomsg = errmsg) write(FRAGGLE_LOG_UNIT, *, err = 667, iomsg = errmsg) "Fraggle logfile" end if + close(FRAGGLE_LOG_UNIT) return diff --git a/src/fraggle/fraggle_regime.f90 b/src/fraggle/fraggle_regime.f90 index bdc11f7ef..eb614d951 100644 --- a/src/fraggle/fraggle_regime.f90 +++ b/src/fraggle/fraggle_regime.f90 @@ -64,7 +64,7 @@ module subroutine fraggle_regime_colliders(self, frag, system, param) 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 - call fraggle_io_log_regime(param, colliders, frag) + call fraggle_io_log_regime(colliders, frag) end associate return diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index 669d2a91c..8dc81c303 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -9,7 +9,8 @@ module fraggle_classes integer(I4B), parameter :: FRAGGLE_NMASS_DIST = 3 !! Number of mass bins returned by the regime calculation (largest fragment, second largest, and remainder) character(len=*), parameter :: FRAGGLE_LOG_OUT = "fraggle.log" !! Name of log file for Fraggle diagnostic information - integer(I4B), parameter :: FRAGGLE_LOG_UNIT = 88 !! Unit number for Fraggle log file + integer(I4B), parameter :: FRAGGLE_LOG_UNIT = 76 !! Unit number for Fraggle log file + !******************************************************************************************************************************** ! fraggle_colliders class definitions and method interfaces !******************************************************************************************************************************* @@ -107,9 +108,25 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? end subroutine fraggle_generate_fragments - module subroutine fraggle_io_log_regime(param, colliders, frag) + module subroutine fraggle_io_log_generate(frag) + implicit none + class(fraggle_fragments), intent(in) :: frag + end subroutine fraggle_io_log_generate + + module subroutine fraggle_io_log_one_message(message) + implicit none + character(len=*), intent(in) :: message + character(STRMAX) :: errmsg + end subroutine fraggle_io_log_one_message + + module subroutine fraggle_io_log_pl(pl, param) + implicit none + class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object (only the new bodies generated in a collision) + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + end subroutine fraggle_io_log_pl + + module subroutine fraggle_io_log_regime(colliders, frag) implicit none - class(swiftest_parameters), intent(in) :: param class(fraggle_colliders), intent(in) :: colliders class(fraggle_fragments), intent(in) :: frag end subroutine fraggle_io_log_regime diff --git a/src/modules/swiftest_globals.f90 b/src/modules/swiftest_globals.f90 index 72d204a49..1f9c6028c 100644 --- a/src/modules/swiftest_globals.f90 +++ b/src/modules/swiftest_globals.f90 @@ -24,6 +24,8 @@ module swiftest_globals real(DP), parameter :: THIRD = 0.333333333333333333333333333333333333333_DP !! Definition of 1 / 3 real(DP), parameter :: DEG2RAD = PI / 180.0_DP !! Definition of conversion factor from degrees to radians real(DP), parameter :: RAD2DEG = 180.0_DP / PI !! Definition of conversion factor from degrees to radians + real(DP), parameter :: GC = 6.6743E-11_DP !! Universal gravitational constant in SI units + real(DP), parameter :: einsteinC = 299792458.0_DP !! Speed of light in SI units integer(I4B), parameter :: LOWERCASE_BEGIN = iachar('a') !! ASCII character set parameter for lower to upper conversion - start of lowercase integer(I4B), parameter :: LOWERCASE_END = iachar('z') !! ASCII character set parameter for lower to upper conversion - end of lowercase @@ -69,7 +71,6 @@ module swiftest_globals integer(I4B), parameter :: USAGE = -2 !! Symbolic name for function return/flag code for printing the usage message integer(I4B), parameter :: HELP = -3 !! Symbolic name for function return/flag code for printing the usage message - character(*), parameter :: SUCCESS_MSG = '(/, "Normal termination of Swiftest (version ", f3.1, ")")' character(*), parameter :: FAIL_MSG = '(/, "Terminating Swiftest (version ", f3.1, ") due to error!!")' character(*), parameter :: USAGE_MSG = '("Usage: swiftest [bs|helio|ra15|rmvs|symba|tu4|whm] ")' @@ -126,15 +127,11 @@ module swiftest_globals character(*), parameter :: PARTICLE_OUTFILE = 'particle.dat' integer(I4B), parameter :: PARTICLEUNIT = 44 !! File unit number for the binary particle info output file - !> Miscellaneous constants: integer(I4B), parameter :: NDIM = 3 !! Number of dimensions in our reality integer(I4B), parameter :: NDIM2 = 2 * NDIM !! 2x the number of dimensions real(DP), parameter :: VSMALL = 2 * epsilon(1._DP) !! Very small number used to prevent floating underflow - real(DP), parameter :: GC = 6.6743E-11_DP !! Universal gravitational constant in SI units - real(DP), parameter :: einsteinC = 299792458.0_DP !! Speed of light in SI units - !> NetCDF variable names and constants character(*), parameter :: NETCDF_OUTFILE = 'bin.nc' !! Default output file name character(*), parameter :: TIME_DIMNAME = "time" !! NetCDF name of the time dimension diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index b929b5040..21986b920 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -322,16 +322,6 @@ module function symba_collision_casemerge(system, param, colliders, frag) result integer(I4B) :: status !! Status flag assigned to this outcome end function symba_collision_casemerge - 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 - 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) use swiftest_classes, only : swiftest_parameters implicit none diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index a3bf27e99..5327fefd9 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -21,14 +21,20 @@ module function symba_collision_casedisruption(system, param, colliders, frag) ! Internals integer(I4B) :: i, nfrag logical :: lfailure - character(len=STRMAX) :: collider_message - - collider_message = "Disruption between" - call symba_collision_collider_message(system%pl, colliders%idx, collider_message) - write(*,*) trim(adjustl(collider_message)) + character(len=STRMAX) :: message + + select case(frag%regime) + case(COLLRESOLVE_REGIME_DISRUPTION) + message = "Disruption between" + nfrag = NFRAG_DISRUPT + case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + message = "Supercatastrophic disruption between" + nfrag = NFRAG_SUPERCAT + end select + call symba_collision_collider_message(system%pl, colliders%idx, message) + call fraggle_io_log_one_message(message) ! Collisional fragments will be uniformly distributed around the pre-impact barycenter - nfrag = NFRAG_DISRUPT call frag%setup(nfrag, param) call frag%set_mass_dist(colliders) frag%id(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)] @@ -38,7 +44,7 @@ module function symba_collision_casedisruption(system, param, colliders, frag) call frag%generate_fragments(colliders, system, param, lfailure) if (lfailure) then - write(*,*) 'No fragment solution found, so treat as a pure hit-and-run' + call fraggle_io_log_one_message("No fragment solution found, so treat as a pure hit-and-run") status = ACTIVE nfrag = 0 select type(pl => system%pl) @@ -49,8 +55,14 @@ module function symba_collision_casedisruption(system, param, colliders, frag) end select else ! Populate the list of new bodies - write(*,'("Generating ",I2.0," fragments")') nfrag - status = DISRUPTION + write(message, *) nfrag + call fraggle_io_log_one_message("Generating " // trim(adjustl(message)) // " fragments") + select case(frag%regime) + case(COLLRESOLVE_REGIME_DISRUPTION) + status = DISRUPTION + case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + status = SUPERCATASTROPHIC + end select call symba_collision_mergeaddsub(system, param, colliders, frag, status) end if @@ -74,12 +86,12 @@ module function symba_collision_casehitandrun(system, param, colliders, frag) r ! Internals integer(I4B) :: i, ibiggest, nfrag, jtarg, jproj logical :: lpure - character(len=STRMAX) :: collider_message + character(len=STRMAX) :: message character(len=NAMELEN) :: idstr - collider_message = "Hit and run between" - call symba_collision_collider_message(system%pl, colliders%idx, collider_message) - write(*,*) trim(adjustl(collider_message)) + message = "Hit and run between" + call symba_collision_collider_message(system%pl, colliders%idx, message) + call fraggle_io_log_one_message(trim(adjustl(message))) if (colliders%mass(1) > colliders%mass(2)) then jtarg = 1 @@ -90,7 +102,7 @@ module function symba_collision_casehitandrun(system, param, colliders, frag) r end if 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.' + call fraggle_io_log_one_message("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 @@ -107,10 +119,11 @@ module function symba_collision_casehitandrun(system, param, colliders, frag) r call frag%generate_fragments(colliders, system, param, lpure) if (lpure) then - write(*,*) 'Should have been a pure hit and run instead' + call fraggle_io_log_one_message("Should have been a pure hit and run instead") nfrag = 0 else - write(*,'("Generating ",I2.0," fragments")') nfrag + write(message, *) nfrag + call fraggle_io_log_one_message("Generating " // trim(adjustl(message)) // " fragments") end if end if if (lpure) then ! Reset these bodies back to being active so that nothing further is done to them @@ -133,7 +146,7 @@ end function symba_collision_casehitandrun 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. + !! Merge massive bodies. !! !! Adapted from David E. Kaufmann's Swifter routines symba_merge_pl.f90 and symba_discard_merge_pl.f90 !! @@ -151,12 +164,12 @@ module function symba_collision_casemerge(system, param, colliders, frag) resul real(DP), dimension(2) :: volume, density real(DP) :: pe real(DP), dimension(NDIM) :: L_spin_new - character(len=STRMAX) :: collider_message character(len=NAMELEN) :: idstr + character(len=STRMAX) :: message - collider_message = "Merging" - call symba_collision_collider_message(system%pl, colliders%idx, collider_message) - write(*,*) trim(adjustl(collider_message)) + message = "Merging" + call symba_collision_collider_message(system%pl, colliders%idx, message) + call fraggle_io_log_one_message(message) select type(pl => system%pl) class is (symba_pl) @@ -215,59 +228,6 @@ module function symba_collision_casemerge(system, param, colliders, frag) resul end function symba_collision_casemerge - 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 - 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 outcom - ! Internals - integer(I4B) :: i, j, nfrag - logical :: lfailure - character(len=STRMAX) :: collider_message - character(len=NAMELEN) :: idstr - - collider_message = "Supercatastrophic disruption between" - call symba_collision_collider_message(system%pl, colliders%idx, collider_message) - write(*,*) trim(adjustl(collider_message)) - - nfrag = NFRAG_SUPERCAT - 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) - - ! 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' - status = ACTIVE - nfrag = 0 - select type(pl => system%pl) - class is (symba_pl) - 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, colliders, frag, status) - end if - - return - end function symba_collision_casesupercatastrophic - - subroutine symba_collision_collider_message(pl, collidx, collider_message) !! author: David A. Minton !! @@ -320,7 +280,8 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec integer(I4B) :: i, j, k, nenc real(DP) :: rlim, Gmtot logical :: isplpl - character(len=STRMAX) :: timestr, idstri, idstrj + character(len=STRMAX) :: timestr, idstri, idstrj, message + lany_collision = .false. if (self%nenc == 0) return @@ -401,9 +362,10 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec write(idstrj, *) tp%id(j) write(timestr, *) t call tp%info(j)%set_value(status="DISCARDED_PLR", discard_time=t, discard_xh=tp%xh(:,j), discard_vh=tp%vh(:,j)) - write(*, *) "Particle " // trim(adjustl(tp%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & + write(message, *) "Particle " // trim(adjustl(tp%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & // " collided with massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & // " at t = " // trim(adjustl(timestr)) + call fraggle_io_log_one_message(message) end if end if end do @@ -752,6 +714,7 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) class(symba_pl), allocatable :: plnew, plsub character(*), parameter :: FRAGFMT = '("Newbody",I0.7)' character(len=NAMELEN) :: newname + character(len=STRMAX) :: message select type(pl => system%pl) class is (symba_pl) @@ -852,6 +815,9 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) plnew%ldiscard(1:nfrag) = .false. plnew%levelg(1:nfrag) = pl%levelg(ibiggest) plnew%levelm(1:nfrag) = pl%levelm(ibiggest) + + ! Log the properties of the new bodies + call fraggle_io_log_pl(plnew, param) ! Append the new merged body to the list nstart = pl_adds%nbody + 1 @@ -924,10 +890,8 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) call colliders%regime(frag, system, param) select case (frag%regime) - case (COLLRESOLVE_REGIME_DISRUPTION) + case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) plplcollision_list%status(i) = symba_collision_casedisruption(system, param, colliders, frag) - case (COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - 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, colliders, frag) case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) @@ -1024,7 +988,10 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, ir do write(timestr,*) t - write(*, *) "Collision between massive bodies detected at time t = " // trim(adjustl(timestr)) + call fraggle_io_log_one_message("") + call fraggle_io_log_one_message("--------------------------------------------------------------------") + call fraggle_io_log_one_message("Collision between massive bodies detected at time t = " // trim(adjustl(timestr))) + call fraggle_io_log_one_message("--------------------------------------------------------------------") allocate(tmp_param, source=param) tmp_param%t = t if (param%lfragmentation) then diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 24f40b1bd..704d11f85 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -20,7 +20,7 @@ subroutine symba_discard_cb_pl(pl, system, param) ! Internals integer(I4B) :: i, j real(DP) :: energy, vb2, rb2, rh2, rmin2, rmax2, rmaxu2 - character(len=STRMAX) :: idstr, timestr + character(len=STRMAX) :: idstr, timestr, message associate(npl => pl%nbody, cb => system%cb) call system%set_msys() @@ -36,7 +36,12 @@ subroutine symba_discard_cb_pl(pl, system, param) pl%status(i) = DISCARDED_RMAX write(idstr, *) pl%id(i) write(timestr, *) param%t - write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too far from the central body at t = " // trim(adjustl(timestr)) + write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too far from the central body at t = " // trim(adjustl(timestr)) + call fraggle_io_log_one_message("") + call fraggle_io_log_one_message("***********************************************************************************************************************") + call fraggle_io_log_one_message(message) + call fraggle_io_log_one_message("***********************************************************************************************************************") + call fraggle_io_log_one_message("") call pl%info(i)%set_value(status="DISCARDED_RMAX", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i)) else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then pl%ldiscard(i) = .true. @@ -44,7 +49,12 @@ subroutine symba_discard_cb_pl(pl, system, param) pl%status(i) = DISCARDED_RMIN write(idstr, *) pl%id(i) write(timestr, *) param%t - write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too close to the central body at t = " // trim(adjustl(timestr)) + write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too close to the central body at t = " // trim(adjustl(timestr)) + call fraggle_io_log_one_message("") + call fraggle_io_log_one_message("***********************************************************************************************************************") + call fraggle_io_log_one_message(message) + call fraggle_io_log_one_message("***********************************************************************************************************************") + call fraggle_io_log_one_message("") call pl%info(i)%set_value(status="DISCARDED_RMIN", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=cb%id) else if (param%rmaxu >= 0.0_DP) then rb2 = dot_product(pl%xb(:,i), pl%xb(:,i)) @@ -56,7 +66,12 @@ subroutine symba_discard_cb_pl(pl, system, param) pl%status(i) = DISCARDED_RMAXU write(idstr, *) pl%id(i) write(timestr, *) param%t - write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) + write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) + call fraggle_io_log_one_message("") + call fraggle_io_log_one_message("***********************************************************************************************************************") + call fraggle_io_log_one_message(message) + call fraggle_io_log_one_message("***********************************************************************************************************************") + call fraggle_io_log_one_message("") call pl%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i)) end if end if diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index d7741e38f..42f36ea64 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -99,8 +99,6 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms write(*,*) "SEED: N,VAL = ",size(param%seed), param%seed(:) if (param%min_GMfrag < 0.0_DP) param%min_GMfrag = param%GMTINY write(*,*) "MIN_GMFRAG = ", self%min_GMfrag - ! For the Fraggle log file, delete it if this is a new run and the file exists - call fraggle_io_log_start(param) end if if (.not.self%lclose) then @@ -108,6 +106,8 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms iostat = -1 return end if + ! All reporting of collision information in SyMBA (including mergers) is now recorded in the Fraggle logfile + call fraggle_io_log_start(param) end associate iostat = 0