From 25b27a9d7f44873231775984003e4290e6e63577 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 23 Aug 2021 23:22:50 -0400 Subject: [PATCH 1/2] Reduced energy error tolerance for faster solutions. --- src/fragmentation/fragmentation.f90 | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index 754a65b10..4787b025c 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -42,7 +42,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, integer(I4B), parameter :: NFRAG_MIN = 7 !! The minimum allowable number of fragments (set to 6 because that's how many unknowns are needed in the tangential velocity calculation) real(DP) :: r_max_start, r_max_start_old, r_max, f_spin real(DP), parameter :: Ltol = 10 * epsilon(1.0_DP) - real(DP), parameter :: Etol = 1e-10_DP + real(DP), parameter :: Etol = 1e-8_DP integer(I4B), parameter :: MAXTRY = 3000 logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes class(swiftest_nbody_system), allocatable :: tmpsys @@ -126,17 +126,17 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, if (.not.lfailure) call set_fragment_tan_vel(lfailure) if (lfailure) then - !write(*,*) 'Failed to find tangential velocities' + ! write(*,*) 'Failed to find tangential velocities' else call set_fragment_radial_velocities(lfailure) - !if (lfailure) write(*,*) 'Failed to find radial velocities' + ! if (lfailure) write(*,*) 'Failed to find radial velocities' if (.not.lfailure) then call calculate_system_energy(linclude_fragments=.true.) if ((abs(dEtot + Qloss) > Etol) .or. (dEtot > 0.0_DP)) then - !write(*,*) 'Failed due to high energy error: ',dEtot, abs(dEtot + Qloss) / Etol + ! write(*,*) 'Failed due to high energy error: ',dEtot, abs(dEtot + Qloss) / Etol lfailure = .true. else if (abs(dLmag) / Lmag_before > Ltol) then - !write(*,*) 'Failed due to high angular momentum error: ', dLmag / Lmag_before + ! write(*,*) 'Failed due to high angular momentum error: ', dLmag / Lmag_before lfailure = .true. end if end if @@ -153,11 +153,11 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, ! write(*, "(' Final diagnostic')") ! write(*, "(' -------------------------------------------------------------------------------------')") ! call calculate_system_energy(linclude_fragments=.true.) - ! if (lfailure) then - ! write(*,*) "symba_frag_pos failed after: ",try," tries" - ! do ii = 1, nfrag - ! vb_frag(:, ii) = vcom(:) - ! end do + if (lfailure) then + ! write(*,*) "symba_frag_pos failed after: ",try," tries" + do ii = 1, nfrag + vb_frag(:, ii) = vcom(:) + end do ! else ! write(*,*) "symba_frag_pos succeeded after: ",try," tries" ! write(*, "(' dL_tot should be very small' )") @@ -166,7 +166,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, ! write(*,fmtlabel) ' dE_tot |', dEtot / abs(Etot_before) ! write(*,fmtlabel) ' Qloss |', -Qloss / abs(Etot_before) ! write(*,fmtlabel) ' dE - Qloss |', (Etot_after - Etot_before + Qloss) / abs(Etot_before) - ! end if + end if ! write(*, "(' -------------------------------------------------------------------------------------')") call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily From e66ca96b77a84a08ad0a119d2cf74dfcf2cc56ac Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 24 Aug 2021 00:40:32 -0400 Subject: [PATCH 2/2] Fixed memory allocation issue in the orbital element conversion method --- src/io/io.f90 | 2 -- src/orbel/orbel.f90 | 21 +++++++++++---------- 2 files changed, 11 insertions(+), 12 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index 58e233dd1..3123710db 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -1660,8 +1660,6 @@ module subroutine io_write_frame_system(self, iu, param) call pl%write_frame(iu, param) call tp%write_frame(iu, param) - deallocate(cb, pl, tp) - close(iu, err = 667, iomsg = errmsg) return diff --git a/src/orbel/orbel.f90 b/src/orbel/orbel.f90 index c07200b14..0a93ec599 100644 --- a/src/orbel/orbel.f90 +++ b/src/orbel/orbel.f90 @@ -16,7 +16,7 @@ module subroutine orbel_el2xv_vec(self, cb) if (self%nbody == 0) return call self%set_mu(cb) - do i = 1, self%nbody + do concurrent (i = 1:self%nbody) call orbel_el2xv(self%mu(i), self%a(i), self%e(i), self%inc(i), self%capom(i), & self%omega(i), self%capm(i), self%xh(:, i), self%vh(:, i)) end do @@ -288,7 +288,7 @@ real(DP) pure function orbel_flon(e,icapn) end if return - end function orbel_flon + end function orbel_flon !********************************************************************** @@ -358,7 +358,7 @@ real(DP) pure function orbel_fget(e,capn) !write(*,*) 'fget : returning without complete convergence' return - end function orbel_fget + end function orbel_fget !********************************************************************** @@ -878,18 +878,19 @@ module subroutine orbel_xv2el_vec(self, cb) if (self%nbody == 0) return call self%set_mu(cb) - if (.not.allocated(self%a)) allocate(self%a(self%nbody)) - if (.not.allocated(self%e)) allocate(self%e(self%nbody)) - if (.not.allocated(self%inc)) allocate(self%inc(self%nbody)) - if (.not.allocated(self%capom)) allocate(self%capom(self%nbody)) - if (.not.allocated(self%omega)) allocate(self%omega(self%nbody)) - if (.not.allocated(self%capm)) allocate(self%capm(self%nbody)) - do i = 1, self%nbody + if (allocated(self%a)) deallocate(self%a); allocate(self%a(self%nbody)) + if (allocated(self%e)) deallocate(self%e); allocate(self%e(self%nbody)) + if (allocated(self%inc)) deallocate(self%inc); allocate(self%inc(self%nbody)) + if (allocated(self%capom)) deallocate(self%capom); allocate(self%capom(self%nbody)) + if (allocated(self%omega)) deallocate(self%omega); allocate(self%omega(self%nbody)) + if (allocated(self%capm)) deallocate(self%capm); allocate(self%capm(self%nbody)) + do concurrent (i = 1:self%nbody) call orbel_xv2el(self%mu(i), self%xh(:, i), self%vh(:, i), self%a(i), self%e(i), self%inc(i), & self%capom(i), self%omega(i), self%capm(i)) end do end subroutine orbel_xv2el_vec + module pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) !! author: David A. Minton !!