diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 01fc3cf18..ad6ef5c87 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -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) @@ -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 @@ -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 diff --git a/src/fraggle/fraggle_io.f90 b/src/fraggle/fraggle_io.f90 new file mode 100644 index 000000000..dbe412167 --- /dev/null +++ b/src/fraggle/fraggle_io.f90 @@ -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 \ No newline at end of file diff --git a/src/fraggle/fraggle_regime.f90 b/src/fraggle/fraggle_regime.f90 index ab786917e..bdc11f7ef 100644 --- a/src/fraggle/fraggle_regime.f90 +++ b/src/fraggle/fraggle_regime.f90 @@ -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 @@ -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 diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index e7b327e93..7f35cb0d2 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -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 diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index d6ea0ac32..669d2a91c 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -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 !******************************************************************************************************************************* @@ -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 diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 2694fcc25..363ed6f66 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -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 @@ -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. @@ -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 diff --git a/src/obl/obl.f90 b/src/obl/obl.f90 index 035a54b18..977fc620c 100644 --- a/src/obl/obl.f90 +++ b/src/obl/obl.f90 @@ -106,7 +106,8 @@ module subroutine obl_acc_tp(self, system) end subroutine obl_acc_tp - module subroutine obl_pot(npl, Mcb, Mpl, j2rp2, j4rp4, xh, irh, oblpot) + + module subroutine obl_pot_system(self) !! author: David A. Minton !! !! Compute the contribution to the total gravitational potential due solely to the oblateness of the central body @@ -118,33 +119,59 @@ module subroutine obl_pot(npl, Mcb, Mpl, j2rp2, j4rp4, xh, irh, oblpot) !! Adapted from Hal Levison's Swift routine obl_pot.f implicit none ! Arguments - 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 + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + ! Internals + integer(I4B) :: i + real(DP), dimension(self%pl%nbody) :: oblpot_arr + + associate(system => self, pl => self%pl, npl => self%pl%nbody, cb => self%cb) + if (.not. any(pl%lmask(1:npl))) return + do concurrent (i = 1:npl, pl%lmask(i)) + oblpot_arr(i) = obl_pot_one(npl, cb%Gmass, pl%Gmass(i), cb%j2rp2, cb%j4rp4, pl%xh(3,i), 1.0_DP / norm2(pl%xh(:,i))) + end do + system%oblpot = sum(oblpot_arr, pl%lmask(1:npl)) + end associate + + return + end subroutine obl_pot_system + + + elemental function obl_pot_one(npl, GMcb, GMpl, j2rp2, j4rp4, zh, irh) result(oblpot) + !! author: David A. Minton + !! + !! Compute the contribution to the total gravitational potential due solely to the oblateness of the central body from a single massive body + !! Returned value does not include monopole term or terms higher than J4 + !! + !! Reference: MacMillan, W. D. 1958. The Theory of the Potential, (Dover Publications), 363. + !! + !! Adapted from David E. Kaufmann's Swifter routine: obl_pot.f90 + !! Adapted from Hal Levison's Swift routine obl_pot.f + implicit none + ! Arguments + integer(I4B), intent(in) :: npl !! Number of massive bodies + real(DP), intent(in) :: GMcb !! G*mass of the central body + real(DP), intent(in) :: GMpl !! G*mass of the massive body + real(DP), intent(in) :: j2rp2 !! J_2 / R**2 of the central body + real(DP), intent(in) :: j4rp4 !! J_2 / R**4 of the central body + real(DP), intent(in) :: zh !! z-component of the heliocentric distance vector of the massive body + real(DP), intent(in) :: irh !! Inverse of the heliocentric distance magnitude of the massive body + ! Result + real(DP) :: oblpot !! Gravitational potential ! Internals integer(I4B) :: i real(DP) :: rinv2, t0, t1, t2, t3, p2, p4, mu - oblpot = 0.0_DP - mu = Mcb - do i = 1, npl - rinv2 = irh(i)**2 - t0 = mu * Mpl(i) * rinv2 * irh(i) - t1 = j2rp2 - t2 = xh(3, i) * xh(3, i) * rinv2 - t3 = j4rp4 * rinv2 - p2 = 0.5_DP * (3 * t2 - 1.0_DP) - p4 = 0.125_DP * ((35 * t2 - 30.0_DP) * t2 + 3.0_DP) - oblpot = oblpot + t0 * (t1 * p2 + t3 * p4) - end do + rinv2 = irh**2 + t0 = GMcb * GMpl * rinv2 * irh + t1 = j2rp2 + t2 = zh**2 * rinv2 + t3 = j4rp4 * rinv2 + p2 = 0.5_DP * (3 * t2 - 1.0_DP) + p4 = 0.125_DP * ((35 * t2 - 30.0_DP) * t2 + 3.0_DP) + oblpot = t0 * (t1 * p2 + t3 * p4) return - end subroutine obl_pot - + end function obl_pot_one end submodule s_obl diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 71717aabe..d7741e38f 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -22,9 +22,9 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms integer(I4B) :: ilength, ifirst, ilast !! Variables used to parse input file character(STRMAX) :: line !! Line of the input file character (len=:), allocatable :: line_trim,param_name, param_value !! Strings used to parse the param file - integer(I4B) :: nseeds, nseeds_from_file, i - logical :: seed_set = .false. !! Is the random seed set in the input file? - character(len=*),parameter :: linefmt = '(A)' + integer(I4B) :: nseeds, nseeds_from_file, i + logical :: seed_set = .false. !! Is the random seed set in the input file? + character(len=*),parameter :: linefmt = '(A)' associate(param => self) call io_param_reader(param, unit, iotype, v_list, iostat, iomsg) @@ -81,7 +81,6 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms end do 1 continue - if (self%GMTINY < 0.0_DP) then write(iomsg,*) "GMTINY invalid or not set: ", self%GMTINY iostat = -1 @@ -100,6 +99,8 @@ 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 diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index 70fbb16df..803b7bac1 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -19,7 +19,6 @@ module subroutine util_get_energy_momentum_system(self, param) real(DP), dimension(self%pl%nbody) :: irh, kepl, kespinpl, pecb real(DP), dimension(self%pl%nbody) :: Lplorbitx, Lplorbity, Lplorbitz real(DP), dimension(self%pl%nbody) :: Lplspinx, Lplspiny, Lplspinz - logical, dimension(self%pl%nbody) :: lstatus real(DP), dimension(NDIM) :: Lcborbit, Lcbspin real(DP), dimension(:), allocatable :: pepl logical, dimension(:), allocatable :: lstatpl @@ -41,13 +40,13 @@ module subroutine util_get_energy_momentum_system(self, param) Lplspiny(:) = 0.0_DP Lplspinz(:) = 0.0_DP - lstatus(1:npl) = pl%status(1:npl) /= INACTIVE + pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE - system%GMtot = cb%Gmass + sum(pl%Gmass(1:npl), lstatus(1:npl)) + system%GMtot = cb%Gmass + sum(pl%Gmass(1:npl), pl%lmask(1:npl)) kecb = cb%mass * dot_product(cb%vb(:), cb%vb(:)) Lcborbit(:) = cb%mass * (cb%xb(:) .cross. cb%vb(:)) - do concurrent (i = 1:npl, lstatus(i)) + do concurrent (i = 1:npl, pl%lmask(i)) hx = pl%xb(2,i) * pl%vb(3,i) - pl%xb(3,i) * pl%vb(2,i) hy = pl%xb(3,i) * pl%vb(1,i) - pl%xb(1,i) * pl%vb(3,i) hz = pl%xb(1,i) * pl%vb(2,i) - pl%xb(2,i) * pl%vb(1,i) @@ -67,7 +66,7 @@ module subroutine util_get_energy_momentum_system(self, param) ! For simplicity, we always assume that the rotation pole is the 3rd principal axis Lcbspin(:) = cb%Ip(3) * cb%mass * cb%radius**2 * cb%rot(:) - do concurrent (i = 1:npl, lstatus(i)) + do concurrent (i = 1:npl, pl%lmask(i)) ! Currently we assume that the rotation pole is the 3rd principal axis ! Angular momentum from spin Lplspinx(i) = pl%mass(i) * pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(1,i) @@ -83,11 +82,11 @@ module subroutine util_get_energy_momentum_system(self, param) end if ! Do the central body potential energy component first - where(.not. lstatus(1:npl)) + where(.not. pl%lmask(1:npl)) pecb(1:npl) = 0.0_DP end where - do concurrent(i = 1:npl, lstatus(i)) + do concurrent(i = 1:npl, pl%lmask(i)) pecb(i) = -cb%Gmass * pl%mass(i) / norm2(pl%xb(:,i)) end do @@ -97,7 +96,7 @@ module subroutine util_get_energy_momentum_system(self, param) do concurrent (k = 1:nplpl) i = pl%k_plpl(1,k) j = pl%k_plpl(2,k) - lstatpl(k) = (lstatus(i) .and. lstatus(j)) + lstatpl(k) = (pl%lmask(i) .and. pl%lmask(j)) end do where(.not.lstatpl(1:nplpl)) @@ -110,29 +109,26 @@ module subroutine util_get_energy_momentum_system(self, param) pepl(k) = -(pl%Gmass(i) * pl%mass(j)) / norm2(pl%xb(:, i) - pl%xb(:, j)) end do - system%pe = sum(pepl(:), lstatpl(:)) + sum(pecb(1:npl), lstatus(1:npl)) + system%pe = sum(pepl(:), lstatpl(:)) + sum(pecb(1:npl), pl%lmask(1:npl)) deallocate(lstatpl, pepl) - system%ke_orbit = 0.5_DP * (kecb + sum(kepl(1:npl), lstatus(:))) - if (param%lrotation) system%ke_spin = 0.5_DP * (kespincb + sum(kespinpl(1:npl), lstatus(:))) + system%ke_orbit = 0.5_DP * (kecb + sum(kepl(1:npl), pl%lmask(1:npl))) + if (param%lrotation) system%ke_spin = 0.5_DP * (kespincb + sum(kespinpl(1:npl), pl%lmask(1:npl))) ! Potential energy from the oblateness term if (param%loblatecb) then - do concurrent(i = 1:npl, lstatus(i)) - irh(i) = 1.0_DP / norm2(pl%xh(:,i)) - end do - call obl_pot(npl, cb%Gmass, pl%mass, cb%j2rp2, cb%j4rp4, pl%xh, irh, oblpot) - system%pe = system%pe + oblpot + call system%obl_pot() + system%pe = system%pe + system%oblpot end if - system%Lorbit(1) = Lcborbit(1) + sum(Lplorbitx(1:npl), lstatus(1:npl)) - system%Lorbit(2) = Lcborbit(2) + sum(Lplorbity(1:npl), lstatus(1:npl)) - system%Lorbit(3) = Lcborbit(3) + sum(Lplorbitz(1:npl), lstatus(1:npl)) + system%Lorbit(1) = Lcborbit(1) + sum(Lplorbitx(1:npl), pl%lmask(1:npl)) + system%Lorbit(2) = Lcborbit(2) + sum(Lplorbity(1:npl), pl%lmask(1:npl)) + system%Lorbit(3) = Lcborbit(3) + sum(Lplorbitz(1:npl), pl%lmask(1:npl)) if (param%lrotation) then - system%Lspin(1) = Lcbspin(1) + sum(Lplspinx(1:npl), lstatus(1:npl)) - system%Lspin(2) = Lcbspin(2) + sum(Lplspiny(1:npl), lstatus(1:npl)) - system%Lspin(3) = Lcbspin(3) + sum(Lplspinz(1:npl), lstatus(1:npl)) + system%Lspin(1) = Lcbspin(1) + sum(Lplspinx(1:npl), pl%lmask(1:npl)) + system%Lspin(2) = Lcbspin(2) + sum(Lplspiny(1:npl), pl%lmask(1:npl)) + system%Lspin(3) = Lcbspin(3) + sum(Lplspinz(1:npl), pl%lmask(1:npl)) end if system%te = system%ke_orbit + system%ke_spin + system%pe