Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Rearranged the Fraggle generation code a bit, improved the stability …
Browse files Browse the repository at this point in the history 
…of the obl_pot call and made it a method, and added a logfile to Fraggle results.
  • Loading branch information
@daminton
daminton committed Sep 3, 2021
1 parent 0a30daf commit c3d0993
Show file tree
Hide file tree
Showing 9 changed files with 257 additions and 105 deletions.
112 changes: 66 additions & 46 deletions src/fraggle/fraggle_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,14 +21,17 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead?
! Internals
integer(I4B) :: i
integer(I4B) :: try
real(DP) :: r_max_start, f_spin, dEtot, dLmag
integer(I4B), parameter :: MAXTRY = 3000
character(len=*), parameter :: fmtlabel = "(A14,10(ES11.4,1X,:))"
logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes
logical :: lk_plpl

integer(I4B) :: i
integer(I4B) :: try
real(DP) :: r_max_start, f_spin, dEtot, dLmag
integer(I4B), parameter :: MAXTRY = 3000
character(len=*), parameter :: fmtlabel = "(A14,10(ES11.4,1X,:))"
logical :: lk_plpl
logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes
logical, dimension(size(IEEE_USUAL)) :: fpe_flag

! 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
call ieee_get_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily
fpe_quiet_modes(:) = .false.
call ieee_set_halting_mode(IEEE_ALL,fpe_quiet_modes)
Expand All @@ -54,17 +57,22 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa

! Start out the fragments close to the initial separation distance. This will be increased if there is any overlap or we fail to find a solution
r_max_start = 1 * norm2(colliders%xb(:,2) - colliders%xb(:,1))
try = 1
lfailure = .false.
try = 1
do while (try < MAXTRY)
if (lfailure) then
call frag%restructure(colliders, try, f_spin, r_max_start)
call frag%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)

! Initial velocity guess will be the barycentric velocity of the colliding system.
! This ensures that our budgets are based on the residual velocities needed in the
! collisional frame
! 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(:)
end do
Expand All @@ -73,54 +81,66 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa
call frag%set_budgets(colliders)

call fraggle_generate_spins(frag, colliders, f_spin, lfailure)
if (.not.lfailure) then
call fraggle_generate_tan_vel(frag, colliders, lfailure)
if (.not. lfailure) then
call fraggle_generate_rad_vel(frag, colliders, lfailure)
! if (lfailure) write(*,*) 'Failed to find radial velocities'
if (.not.lfailure) then
call frag%get_energy_and_momentum(colliders, system, param, lbefore=.false.)
dEtot = frag%Etot_after - frag%Etot_before
dLmag = .mag. (frag%Ltot_after(:) - frag%Ltot_before(:))

if ((abs(dEtot + frag%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) then
! write(*,*) 'Failed due to high energy error: ',dEtot, abs(dEtot + frag%Qloss) / FRAGGLE_ETOL
lfailure = .true.
else if ((abs(dLmag) / (.mag.frag%Ltot_before)) > FRAGGLE_LTOL) then
! write(*,*) 'Failed due to high angular momentum error: ', dLmag / (.mag.frag%Ltot_before(:))
lfailure = .true.
end if
end if
end if
if (lfailure) then
write(*,*) "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"
cycle
end if

call fraggle_generate_rad_vel(frag, colliders, lfailure)
if (lfailure) then
write(*,*) "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(:))

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
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(:))
cycle
end if

! Check if any of the usual floating point exceptions happened, and fail the try if so
call ieee_get_flag(ieee_usual, fpe_flag)
lfailure = any(fpe_flag)
if (.not.lfailure) exit
call frag%restructure(colliders, try, f_spin, r_max_start)
call frag%reset()
try = try + 1
write(*,*) "Fraggle failed due to a floating point exception: ", fpe_flag
end do

! write(*, "(' -------------------------------------------------------------------------------------')")
! write(*, "(' Final diagnostic')")
! write(*, "(' -------------------------------------------------------------------------------------')")
write(*, "(' -------------------------------------------------------------------------------------')")
write(*, "(' Fraggle final diagnostic')")
write(*, "(' -------------------------------------------------------------------------------------')")
if (lfailure) then
write(*,*) "Fraggle fragment generation failed after: ",try," tries"
else
write(*,*) "Fraggle fragment generation succeeded after: ",try," tries"
! write(*, "(' dL_tot should be very small' )")
! write(*,fmtlabel) ' dL_tot |', (.mag.(frag%Ltot_after(:) - frag%Ltot_before(:))) / (.mag.frag%Ltot_before(:))
! write(*, "(' dE_tot should be negative and equal to Qloss' )")
! write(*,fmtlabel) ' dE_tot |', (frag%Etot_after - frag%Etot_before) / abs(frag%Etot_before)
! write(*,fmtlabel) ' Qloss |', -frag%Qloss / abs(frag%Etot_before)
! write(*,fmtlabel) ' dE - Qloss |', (frag%Etot_after - frag%Etot_before + frag%Qloss) / abs(frag%Etot_before)
write(*, "(' dL_tot should be very small' )")
write(*,fmtlabel) ' dL_tot |', (.mag.(frag%Ltot_after(:) - frag%Ltot_before(:))) / (.mag.frag%Ltot_before(:))
write(*, "(' dE_tot should be negative and equal to Qloss' )")
write(*,fmtlabel) ' dE_tot |', (frag%Etot_after - frag%Etot_before) / abs(frag%Etot_before)
write(*,fmtlabel) ' Qloss |', -frag%Qloss / abs(frag%Etot_before)
write(*,fmtlabel) ' dE - Qloss |', (frag%Etot_after - frag%Etot_before + frag%Qloss) / abs(frag%Etot_before)
end if
! write(*, "(' -------------------------------------------------------------------------------------')")
write(*, "(' -------------------------------------------------------------------------------------')")
call frag%set_original_scale(colliders)

! Restore the big array
if (lk_plpl) call pl%index(param)
end associate

call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily

return
Expand Down
96 changes: 96 additions & 0 deletions src/fraggle/fraggle_io.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,96 @@
submodule(fraggle_classes) s_fraggle_io
use swiftest

contains

module subroutine fraggle_io_log_regime(param, colliders, frag)
!! author: David A. Minton
!!
!! Writes a log of the results of the collisional regime determination
implicit none
! Arguments
class(swiftest_parameters), intent(in) :: param
class(fraggle_colliders), intent(in) :: colliders
class(fraggle_fragments), intent(in) :: frag
! 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, *) "--------------------------------------------------------------------"
write(FRAGGLE_LOG_UNIT, *) " Fraggle collisional regime determination results"
write(FRAGGLE_LOG_UNIT, *) "--------------------------------------------------------------------"
write(FRAGGLE_LOG_UNIT, *) "----------------------- Collider information -----------------------"
write(FRAGGLE_LOG_UNIT, *) "True number of colliders : ",colliders%ncoll
write(FRAGGLE_LOG_UNIT, *) "Index list of true colliders : ",colliders%idx(1:colliders%ncoll)
write(FRAGGLE_LOG_UNIT, *) "-------------------- Two-body equialent values ---------------------"
write(FRAGGLE_LOG_UNIT, *) "mass1 : ",colliders%mass(1)
write(FRAGGLE_LOG_UNIT, *) "radius1 : ",colliders%radius(1)
write(FRAGGLE_LOG_UNIT, *) "xb1 : ",colliders%xb(:,1)
write(FRAGGLE_LOG_UNIT, *) "vb1 : ",colliders%vb(:,1)
write(FRAGGLE_LOG_UNIT, *) "rot1 : ",colliders%rot(:,1)
write(FRAGGLE_LOG_UNIT, *) "Ip1 : ",colliders%Ip(:,1)
write(FRAGGLE_LOG_UNIT, *) "L_spin1 : ",colliders%L_spin(:,1)
write(FRAGGLE_LOG_UNIT, *) "L_orbit1 : ",colliders%L_orbit(:,1)
write(FRAGGLE_LOG_UNIT, *) "mass2 : ",colliders%mass(2)
write(FRAGGLE_LOG_UNIT, *) "radius2 : ",colliders%radius(2)
write(FRAGGLE_LOG_UNIT, *) "xb2 : ",colliders%xb(:,2)
write(FRAGGLE_LOG_UNIT, *) "vb2 : ",colliders%vb(:,2)
write(FRAGGLE_LOG_UNIT, *) "rot2 : ",colliders%rot(:,2)
write(FRAGGLE_LOG_UNIT, *) "Ip2 : ",colliders%Ip(:,2)
write(FRAGGLE_LOG_UNIT, *) "L_spin2 : ",colliders%L_spin(:,2)
write(FRAGGLE_LOG_UNIT, *) "L_orbit2 : ",colliders%L_orbit(:,2)
write(FRAGGLE_LOG_UNIT, *) "------------------------------ Regime -----------------------------"
select case(frag%regime)
case(COLLRESOLVE_REGIME_MERGE)
write(FRAGGLE_LOG_UNIT, *) "Merge"
case(COLLRESOLVE_REGIME_DISRUPTION)
write(FRAGGLE_LOG_UNIT, *) "Disruption"
case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC)
write(FRAGGLE_LOG_UNIT, *) "Supercatastrophic disruption"
case(COLLRESOLVE_REGIME_GRAZE_AND_MERGE)
write(FRAGGLE_LOG_UNIT, *) "Graze and merge"
case(COLLRESOLVE_REGIME_HIT_AND_RUN)
write(FRAGGLE_LOG_UNIT, *) "Hit and run"
end select
write(FRAGGLE_LOG_UNIT, *) "----------------------- Fragment information ----------------------"
write(FRAGGLE_LOG_UNIT, *) "Total mass of fragments : ", frag%mtot
write(FRAGGLE_LOG_UNIT, *) "Largest fragment mass : ", frag%mass_dist(1)
write(FRAGGLE_LOG_UNIT, *) "Second-largest fragment mass : ", frag%mass_dist(2)
write(FRAGGLE_LOG_UNIT, *) "Remaining fragment mass : ", frag%mass_dist(3)
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, *) "--------------------------------------------------------------------"

return
667 continue
write(*,*) "Error writing Fraggle regime information to log file: " // trim(adjustl(errmsg))
end subroutine fraggle_io_log_regime


module subroutine fraggle_io_log_start(param)
!! author: David A. Minton
!!
!! Checks to see if the Fraggle log file needs to be replaced if this is a new run, or appended if this is a restarted run
implicit none
! Arguments
class(swiftest_parameters), intent(in) :: param
! Internals
character(STRMAX) :: errmsg
logical :: fileExists

inquire(file=FRAGGLE_LOG_OUT, exist=fileExists)
if (.not.param%lrestart .or. .not.fileExists) then
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

return

667 continue
write(*,*) "Error writing Fraggle log file: " // trim(adjustl(errmsg))
end subroutine fraggle_io_log_start

end submodule s_fraggle_io
3 changes: 3 additions & 0 deletions src/fraggle/fraggle_regime.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ module subroutine fraggle_regime_colliders(self, frag, system, param)
real(DP) :: min_mfrag_si, Mcb_si
real(DP), dimension(NDIM) :: x1_si, v1_si, x2_si, v2_si
real(DP) :: mlr, mslr, mtot, dentot, msys, msys_new, Qloss, impact_parameter
logical :: fileExists

associate(colliders => self)
! Convert all quantities to SI units and determine which of the pair is the projectile vs. target before sending them to the regime determination subroutine
Expand Down Expand Up @@ -62,6 +63,8 @@ module subroutine fraggle_regime_colliders(self, frag, system, param)
frag%mtot = sum(colliders%mass(:))
frag%xbcom(:) = (colliders%mass(1) * colliders%xb(:,1) + colliders%mass(2) * colliders%xb(:,2)) / frag%mtot
frag%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / frag%mtot

call fraggle_io_log_regime(param, colliders, frag)
end associate

return
Expand Down
1 change: 0 additions & 1 deletion src/main/swiftest_driver.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,6 @@ program swiftest_driver
call nbody_system%dump(param)
end if


!> Define the maximum number of threads
nthreads = 1 ! In the *serial* case
!$ nthreads = omp_get_max_threads() ! In the *parallel* case
Expand Down
16 changes: 15 additions & 1 deletion src/modules/fraggle_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,9 @@ module fraggle_classes
implicit none
public

integer(I4B), parameter :: FRAGGLE_NMASS_DIST = 3 !! Number of mass bins returned by the regime calculation (largest fragment, second largest, and remainder)
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
!********************************************************************************************************************************
! fraggle_colliders class definitions and method interfaces
!*******************************************************************************************************************************
Expand Down Expand Up @@ -105,6 +107,18 @@ 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)
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

module subroutine fraggle_io_log_start(param)
implicit none
class(swiftest_parameters), intent(in) :: param
end subroutine fraggle_io_log_start

!> 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
Expand Down
14 changes: 5 additions & 9 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -396,6 +396,7 @@ module swiftest_classes
real(DP) :: ke_spin = 0.0_DP !! System spin kinetic energy
real(DP) :: pe = 0.0_DP !! System potential energy
real(DP) :: te = 0.0_DP !! System total energy
real(DP) :: oblpot = 0.0_DP !! System potential energy due to oblateness of the central body
real(DP), dimension(NDIM) :: Lorbit = 0.0_DP !! System orbital angular momentum vector
real(DP), dimension(NDIM) :: Lspin = 0.0_DP !! System spin angular momentum vector
real(DP), dimension(NDIM) :: Ltot = 0.0_DP !! System angular momentum vector
Expand All @@ -417,6 +418,7 @@ module swiftest_classes
procedure :: write_frame => io_write_frame_system !! Append a frame of output data to file
procedure :: write_hdr_bin => io_write_hdr_system !! Write a header for an output frame in Fortran binary format
procedure :: write_hdr_netcdf => netcdf_write_hdr_system !! Write a header for an output frame in NetCDF format
procedure :: obl_pot => obl_pot_system !! Compute the contribution to the total gravitational potential due solely to the oblateness of the central body
procedure :: initialize => setup_initialize_system !! Initialize the system from input files
procedure :: init_particle_info => setup_initialize_particle_info_system !! Initialize the system from input files
procedure :: step_spin => tides_step_spin_system !! Steps the spins of the massive & central bodies due to tides.
Expand Down Expand Up @@ -927,16 +929,10 @@ module subroutine obl_acc_tp(self, system)
class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
end subroutine obl_acc_tp

module subroutine obl_pot(npl, Mcb, Mpl, j2rp2, j4rp4, xh, irh, oblpot)
module subroutine obl_pot_system(self)
implicit none
integer(I4B), intent(in) :: npl
real(DP), intent(in) :: Mcb
real(DP), dimension(:), intent(in) :: Mpl
real(DP), intent(in) :: j2rp2, j4rp4
real(DP), dimension(:), intent(in) :: irh
real(DP), dimension(:, :), intent(in) :: xh
real(DP), intent(out) :: oblpot
end subroutine obl_pot
class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object
end subroutine obl_pot_system

module subroutine orbel_el2xv_vec(self, cb)
implicit none
Expand Down
Loading

0 comments on commit c3d0993

Please sign in to comment.