From e09ba73d2b33f92c8e854474c13c906c27768e76 Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 26 May 2021 22:53:30 -0400 Subject: [PATCH] Greatly improved robustness of the fragmentation code in symba_frag_pos. Reduced the number of failed initializations substantially --- src/modules/lambda_function.f90 | 55 +++- src/modules/module_interfaces.f90 | 12 +- src/symba/symba_collision.f90 | 34 ++- src/symba/symba_frag_pos.f90 | 358 ++++++++++++-------------- src/user/user_dump_param.f90 | 1 + src/user/user_udio_reader.f90 | 8 +- src/util/util_minimize_bfgs.f90 | 123 +++++---- src/util/util_solve_linear_system.f90 | 10 +- 8 files changed, 349 insertions(+), 252 deletions(-) diff --git a/src/modules/lambda_function.f90 b/src/modules/lambda_function.f90 index bdd8da74e..88a27342f 100644 --- a/src/modules/lambda_function.f90 +++ b/src/modules/lambda_function.f90 @@ -4,16 +4,29 @@ module lambda_function implicit none type, public :: lambda_obj + !! Base class for an lambda function object. This object takes no additional arguments other than the dependent variable x, an array of real numbers procedure(lambda0), pointer, nopass :: lambdaptr => null() - + real(DP) :: lastval + real(DP),dimension(:), allocatable :: lastarg contains generic :: init => lambda_init_0 procedure :: eval => lambda_eval_0 procedure, nopass :: lambda_init_0 final :: lambda_destroy end type + + type, public, extends(lambda_obj) :: lambda_obj_err + !! Extended class for an lambda function object. This object takes allows for the return of a logical error flag during evaluation of the function. + procedure(lambda0err), pointer, nopass :: lambdaptr_err => null() + logical :: lerr + contains + generic :: init => lambda_init_0_err + procedure :: eval => lambda_eval_0_err + procedure, nopass :: lambda_init_0_err + end type interface lambda_obj module procedure lambda_init_0 + module procedure lambda_init_0_err end interface abstract interface @@ -23,6 +36,13 @@ function lambda0(x) result(y) real(DP), dimension(:), intent(in) :: x real(DP) :: y end function + function lambda0err(x, lerr) result(y) + ! Template for a 0 argument function that returns an error value + import DP + real(DP), dimension(:), intent(in) :: x + logical, intent(out) :: lerr + real(DP) :: y + end function end interface contains @@ -34,22 +54,53 @@ type(lambda_obj) function lambda_init_0(lambda) lambda_init_0%lambdaptr => lambda return end function lambda_init_0 + + type(lambda_obj_err) function lambda_init_0_err(lambda, lerr) + implicit none + ! Arguments + procedure(lambda0err) :: lambda + logical, intent(in) :: lerr + lambda_init_0_err%lambdaptr_err => lambda + lambda_init_0_err%lerr = lerr + return + end function lambda_init_0_err function lambda_eval_0(self, x) result(y) implicit none ! Arguments - class(lambda_obj), intent(in) :: self + class(lambda_obj), intent(inout) :: self real(DP), dimension(:), intent(in) :: x ! Result real(DP) :: y if (associated(self%lambdaptr)) then y = self%lambdaptr(x) + self%lastval = y + if (allocated(self%lastarg)) deallocate(self%lastarg) + allocate(self%lastarg, source=x) else error stop "Lambda function was not initialized" end if end function lambda_eval_0 + function lambda_eval_0_err(self, x) result(y) + implicit none + ! Arguments + class(lambda_obj_err), intent(inout) :: self + real(DP), dimension(:), intent(in) :: x + ! Result + real(DP) :: y + + if (associated(self%lambdaptr_err)) then + y = self%lambdaptr_err(x, self%lerr) + self%lastval = y + if (allocated(self%lastarg)) deallocate(self%lastarg) + allocate(self%lastarg, source=x) + else + error stop "Lambda function was not initialized" + end if + end function lambda_eval_0_err + subroutine lambda_destroy(self) implicit none type(lambda_obj) :: self diff --git a/src/modules/module_interfaces.f90 b/src/modules/module_interfaces.f90 index 8257f0b57..f543a5d6c 100644 --- a/src/modules/module_interfaces.f90 +++ b/src/modules/module_interfaces.f90 @@ -1614,12 +1614,12 @@ function util_minimize_bfgs(f, N, x0_d, eps_d, lerr) result(x1_d) use swiftest_globals use lambda_function implicit none - integer(I4B), intent(in) :: N - class(lambda_obj), intent(in) :: f - real(DP), dimension(:), intent(in) :: x0_d - real(DP), intent(in) :: eps_d - logical, intent(out) :: lerr - real(DP), dimension(:), allocatable :: x1_d + integer(I4B), intent(in) :: N + class(lambda_obj), intent(inout) :: f + real(DP), dimension(:), intent(in) :: x0_d + real(DP), intent(in) :: eps_d + logical, intent(out) :: lerr + real(DP), dimension(:), allocatable :: x1_d end function util_minimize_bfgs end interface diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index e1432eb43..9f29859cc 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -34,7 +34,7 @@ subroutine symba_collision (t, symba_plA, nplplenc, plplenc_list, ldiscard, merg real(DP) :: mchild, mtot real(DP), dimension(NDIM) :: xc, vc, xcom, vcom, xchild, vchild, xcrossv real(DP) :: mtiny_si - real(DP) :: mlr, mslr, msys, msys_new, Qloss + real(DP) :: mlr, mslr, msys, msys_new, Qloss, impact_parameter integer(I4B), dimension(:), allocatable :: family, collision_idx, unique_parent_idx integer(I4B) :: fam_size, istart, ncollisions, nunique_parent type family_array @@ -153,6 +153,9 @@ subroutine symba_collision (t, symba_plA, nplplenc, plplenc_list, ldiscard, merg Ip(:, j) = mass(j) * symba_plA%helio%swiftest%Ip(:, idx_parent(j)) ! Assume principal axis rotation about axis corresponding to highest moment of inertia (3rd Ip) L_spin(:, j) = Ip(3, j) * radius(j)**2 * symba_plA%helio%swiftest%rot(:, idx_parent(j)) + + ! Use conservation of energy and angular momentum to adjust velocities and distances to be those equivalent to where they would be when the mutual radii are touching + !call vector_adjust(mass(j), x(:,j), v(:,j)) if (nchild(j) > 0) then do i = 1, nchild(j) ! Loop over all children and take the mass weighted mean of the properties idx_child = parent_child_index_array(j)%idx(i + 1) @@ -308,4 +311,33 @@ subroutine symba_collision (t, symba_plA, nplplenc, plplenc_list, ldiscard, merg return + contains + subroutine vector_adjust(mass, radius, x, v) + !! author: David A. Minton + !! + !! Uses conservation of angular momentum and energy to adjust the velocity vectors of two bodies to be what they would be if they were touching + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: mass, radius + real(DP), dimension(:), intent(inout) :: x,v + ! Internals + real(DP) :: mrat, vmag0, vmag, rmag0, rmag, B + real(DP), dimension(NDIM) :: h0 + + mrat = mass(1) * mass(2) / sum(mass(1:2)) + vmag0 = norm2(v(:)) + rmag0 = norm2(x(:)) + rmag = sum(radius(1:2)) + + ! Solve vis-viva to get the new magnitude of the velocity + vmag = 2 * sqrt((0.5_DP * vmag0**2 * mrat + mass(1) * mass(2) * (1.0_DP / rmag - 1.0_DP / rmag0)) / mrat) + + ! Use conservation of angular momentum to get the new impact parameter + call util_cross_product(x(:), v(:), h0(:)) + B = norm2(h0) / (rmag * vmag) + + + + end subroutine vector_adjust + end subroutine symba_collision diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index 23f1e7eb7..d2c8e6409 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -1,92 +1,3 @@ -module symba_frag_pos_lambda_implementation - use swiftest_globals - use lambda_function - implicit none - - ! Define the class and interface used to implement the lambda function - type, public, extends(lambda_obj) :: ke_constraint - procedure(abstract_objective_func), pointer, nopass :: ke_objective_func_ptr => null() - real(DP), dimension(:), allocatable :: v_t_mag - real(DP), dimension(:,:), allocatable :: v_r_unit, v_t_unit - real(DP), dimension(:), allocatable :: m_frag - real(DP), dimension(NDIM) :: vcom - real(DP) :: ke_target - contains - generic :: init => ke_objective_func_init - procedure :: eval => ke_objective_func_eval - procedure, nopass :: ke_objective_func_init - final :: ke_objective_func_destroy - end type ke_constraint - interface ke_constraint - module procedure ke_objective_func_init - end interface - - abstract interface - function abstract_objective_func(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom, ke_target) result(fnorm) - ! Template for the kinetic energy constraint function used for minimizing - import DP - real(DP), dimension(:), intent(in) :: v_r_mag !! Radial velocity mangitude - real(DP), dimension(:,:), intent(in) :: v_r_unit !! Radial velocity unit vector - real(DP), dimension(:), intent(in) :: v_t_mag !! Tangential velocity magnitude - real(DP), dimension(:,:), intent(in) :: v_t_unit !! Tangential velocity unit vector - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:), intent(in) :: vcom !! Barycentric velocity of collisional system center of mass - real(DP), intent(in) :: ke_target !! Target kinetic energ - real(DP) :: fnorm !! The objective function result: norm of the vector composed of the tangential momentum and energy - end function - end interface - - contains - type(ke_constraint) function ke_objective_func_init(lambda, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom, ke_target) - implicit none - ! Arguments - procedure(abstract_objective_func) :: lambda !! The lambda function - real(DP), dimension(:,:), intent(in) :: v_r_unit !! Radial velocity unit vector - real(DP), dimension(:), intent(in) :: v_t_mag !! Tangential velocity magnitude - real(DP), dimension(:,:), intent(in) :: v_t_unit !! Tangential velocity unit vector - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:), intent(in) :: vcom !! Barycentric velocity of collisional system center of mass - real(DP), intent(in) :: ke_target !! Target kinetic energ - ! Internals - associate(self => ke_objective_func_init) - self%ke_objective_func_ptr => lambda - allocate(self%v_r_unit, source=v_r_unit) - allocate(self%v_t_mag, source=v_t_mag) - allocate(self%v_t_unit, source=v_t_unit) - allocate(self%m_frag, source=m_frag) - self%vcom(:) = vcom(:) - self%ke_target = ke_target - end associate - return - end function ke_objective_func_init - - subroutine ke_objective_func_destroy(self) - implicit none - type(ke_constraint) :: self - if (allocated(self%v_r_unit)) deallocate(self%v_r_unit) - if (allocated(self%v_t_mag)) deallocate(self%v_t_mag) - if (allocated(self%v_t_unit)) deallocate(self%v_t_unit) - if (allocated(self%m_frag)) deallocate(self%m_frag) - if (associated(self%ke_objective_func_ptr)) nullify(self%ke_objective_func_ptr) - end subroutine ke_objective_func_destroy - - function ke_objective_func_eval(self, x) result(fval) - implicit none - ! Arguments - class(ke_constraint), intent(in) :: self - real(DP), dimension(:), intent(in) :: x - ! Result - real(DP) :: fval - - if (associated(self%ke_objective_func_ptr)) then - fval = self%ke_objective_func_ptr(x, self%v_r_unit, self%v_t_mag, self%v_t_unit, self%m_frag, self%vcom, self%ke_target) - else - error stop "KE Objective function was not initialized." - end if - end function ke_objective_func_eval - -end module symba_frag_pos_lambda_implementation - subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radius, & nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lmerge) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton @@ -98,6 +9,7 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi use module_helio use module_symba use module_swiftestalloc + use lambda_function use module_interfaces, except_this_one => symba_frag_pos implicit none ! Arguments @@ -114,7 +26,7 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi real(DP), intent(inout) :: Qloss ! Internals real(DP), parameter :: f_spin = 0.20_DP !! Fraction of pre-impact orbital angular momentum that is converted to fragment spin - real(DP) :: mscale = 1.0_DP, rscale = 1.0_DP, vscale = 1.0_DP, tscale = 1.0_DP, Lscale = 1.0_DP, Escale = 1.0_DP! Scale factors that reduce quantities to O(~1) in the collisional system + real(DP) :: mscale, rscale, vscale, tscale, Lscale, Escale ! Scale factors that reduce quantities to O(~1) in the collisional system real(DP) :: mtot real(DP), dimension(NDIM) :: xcom, vcom integer(I4B) :: ii, fam_size @@ -132,12 +44,20 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi character(len=*), parameter :: fmtlabel = "(A14,10(ES11.4,1X,:))" integer(I4B) :: try, ntry integer(I4B), parameter :: NFRAG_MIN = 6 !! 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 = 1.0_DP - real(DP), parameter :: Ltol = 2 * epsilon(1.0_DP) + real(DP) :: r_max_start + real(DP), parameter :: Ltol = 10 * epsilon(1.0_DP) real(DP), parameter :: Etol = 1e-8_DP integer(I4B), parameter :: MAXTRY = 1000 - integer(I4B) :: iflip = 1 - + integer(I4B) :: iflip = 1 + logical :: lreduce_fragment_number + + r_max_start = 1.0_DP + mscale = 1.0_DP + rscale = 1.0_DP + vscale = 1.0_DP + tscale = 1.0_DP + Lscale = 1.0_DP + Escale = 1.0_DP if (nfrag < NFRAG_MIN) then write(*,*) "symba_frag_pos needs at least ",NFRAG_MIN," fragments, but only ",nfrag," were given." lmerge = .true. @@ -172,61 +92,20 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi lmerge = .false. do while (nfrag >= NFRAG_MIN .and. try <= MAXTRY) lmerge = .false. - !write(*, "(' -------------------------------------------------------------------------------------')") - !write(*,*) "Try number: ",try, ' of ',ntry - !write(*, "(' -------------------------------------------------------------------------------------')") - !write(*, "(' First pass to get angular momentum ')") + lreduce_fragment_number = .false. call set_fragment_position_vectors() - call set_fragment_tangential_velocities() - call calculate_system_energy(linclude_fragments=.true.) - - !write(*, "(' -------------------------------------------------------------------------------------')") - !write(*, "(' | T_orb T_spin T pe Etot Ltot')") - !write(*, "(' -------------------------------------------------------------------------------------')") - !write(*,fmtlabel) ' change |',(ke_orbit_after - ke_orbit_before) / abs(Etot_before), & - ! (ke_spin_after - ke_spin_before)/ abs(Etot_before), & - ! (ke_orbit_after + ke_spin_after - ke_orbit_before - ke_spin_before)/ abs(Etot_before), & - ! (pe_after - pe_before) / abs(Etot_before), & - ! dEtot, dLmag - ! - ! write(*, "(' -------------------------------------------------------------------------------------')") - ! write(*, "(' If T_offset > 0, this will be a merge')") - ! write(*,fmtlabel) ' T_offset |',ke_offset / abs(Etot_before) - - call set_fragment_radial_velocities(lmerge) - - call calculate_system_energy(linclude_fragments=.true.) - + call set_fragment_tangential_velocities(lmerge) + if (.not.lmerge) call set_fragment_radial_velocities(lmerge) if (.not.lmerge) then + call calculate_system_energy(linclude_fragments=.true.) if (dEtot > 0.0_DP) then - write(*,*) 'Failed try ',try,': Positive energy' lmerge = .true. else if (abs((Etot_after - Etot_before + Qloss) / Etot_before) > Etol) then - write(*,*) 'Failed try ',try,': Energy error too big: ',dEtot lmerge = .true. else if (abs(dLmag) > Ltol) then - write(*,*) 'Failed try ',try,': Angular momentum error too big: ',dLmag lmerge = .true. end if - - !write(*, "(' -------------------------------------------------------------------------------------')") - !write(*, "(' Second pass to get energy ')") - !write(*, "(' -------------------------------------------------------------------------------------')") - !write(*,fmtlabel) ' T_frag targ |',ke_target / abs(Etot_before) - !write(*,fmtlabel) ' T_frag calc |',ke_frag / abs(Etot_before) - !write(*, "(' T_offset should now be small')") - !write(*,fmtlabel) ' T_offset |',ke_offset / abs(Etot_before) - !write(*, "(' -------------------------------------------------------------------------------------')") - !write(*, "(' | T_orb T_spin T pe Etot Ltot')") - !write(*, "(' -------------------------------------------------------------------------------------')") - !write(*,fmtlabel) ' change |',(ke_orbit_after - ke_orbit_before) / abs(Etot_before), & - ! (ke_spin_after - ke_spin_before)/ abs(Etot_before), & - ! (ke_orbit_after + ke_spin_after - ke_orbit_before - ke_spin_before)/ abs(Etot_before), & - ! (pe_after - pe_before) / abs(Etot_before), & - ! dEtot, dLmag - !else - !write(*,*) 'Failed try ',try,': Could not find radial velocities.' end if if (.not.lmerge) exit @@ -238,7 +117,6 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi write(*, "(' -------------------------------------------------------------------------------------')") if (lmerge) then write(*,*) "symba_frag_pos failed after: ",try," tries" - stop else write(*,*) "symba_frag_pos succeeded after: ",try," tries" write(*, "(' dL_tot should be very small' )") @@ -251,6 +129,7 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi write(*, "(' -------------------------------------------------------------------------------------')") call restore_scale_factors() + close(22) return @@ -348,6 +227,12 @@ subroutine define_coordinate_system() allocate(v_r_unit(NDIM,nfrag)) allocate(v_t_unit(NDIM,nfrag)) + rmag(:) = 0.0_DP + v_r_mag(:) = 0.0_DP + v_t_mag(:) = 0.0_DP + v_r_unit(:,:) = 0.0_DP + v_t_unit(:,:) = 0.0_DP + L_orb(:, :) = 0.0_DP ! Compute orbital angular momentum of pre-impact system do i = 1, 2 @@ -489,7 +374,6 @@ subroutine calculate_system_energy(linclude_fragments) ke_frag = ke_frag + 0.5_DP * m_frag(i) * dot_product(vb_frag(:, i), vb_frag(:, i)) end do ke_target = ke_family + (ke_spin_before - ke_spin) + (pe_before - pe) - Qloss - ke_offset = ke_frag - ke_target L_offset(:) = Ltot_before(:) - Ltot_after(:) else Ltot_before(:) = Lorbit(:) + Lspin(:) @@ -542,7 +426,6 @@ subroutine set_fragment_position_vectors() logical, dimension(:), allocatable :: loverlap integer(I4B) :: i, j - allocate(loverlap(nfrag)) ! Place the fragments into a region that is big enough that we should usually not have overlapping bodies @@ -588,26 +471,135 @@ subroutine set_fragment_position_vectors() return end subroutine set_fragment_position_vectors - subroutine set_fragment_tangential_velocities() + subroutine set_fragment_tangential_velocities(lerr) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Adjusts the positions, velocities, and spins of a collection of fragments such that they conserve angular momentum implicit none + ! Arguments + logical, intent(out) :: lerr ! Internals - integer(I4B) :: i - real(DP) :: L_orb_mag - real(DP), dimension(2 * NDIM, 2 * NDIM) :: A ! LHS of linear equation used to solve for momentum constraint in Gauss elimination code - 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 + integer(I4B) :: i + real(DP) :: L_orb_mag, fval + real(DP), parameter :: TOL = 1e-4_DP + real(DP), dimension(:), allocatable :: v_t_initial + type(lambda_obj_err) :: objective_function - v_frag(:,:) = 0.0_DP - ! Divide up the pre-impact spin angular momentum equally between the various bodies by mass do i = 1, nfrag rot_frag(:,i) = L_frag_spin(:) / nfrag / (Ip_frag(:, i) * m_frag(i) * rad_frag(i)**2) end do + vb_frag(:,:) = 0.0_DP + + call calculate_system_energy(linclude_fragments=.true.) + if (ke_target < 0.0_DP) then + write(*,*) 'Negative ke_target: ',ke_target + lerr = .true. + lreduce_fragment_number = .true. + return + end if + + lerr = .false. + if (nfrag > 6) then + allocate(v_t_initial, mold=v_t_mag) + L_orb_mag = norm2(L_frag_orb(:)) + v_t_initial(1:6) = 0.0_DP + do i = 7, nfrag + v_t_initial(i) = L_orb_mag / (m_frag(i) * rmag(i) * nfrag) + end do + v_t_mag(1:nfrag) = solve_fragment_tangential_velocities(v_t_mag_input=v_t_initial(7:nfrag), lerr=lerr) + else + v_t_mag(1:nfrag) = solve_fragment_tangential_velocities(lerr) + end if + vb_frag(:,:) = vmag_to_vb(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom) + ke_frag = 0.0_DP + do i = 1, nfrag + ke_frag = ke_frag + 0.5_DP * m_frag(i) * dot_product(vb_frag(:, i), vb_frag(:, i)) + end do + ke_offset = ke_frag - ke_target + if ((ke_offset > 0.0_DP) .and. (nfrag >6)) then + objective_function = lambda_obj(tangential_objective_function, lerr) + v_t_mag(7:nfrag) = util_minimize_bfgs(objective_function, nfrag - 6, v_t_initial(7:nfrag), TOL, lerr) + if (lerr) v_t_mag(7:nfrag) = objective_function%lastarg(:) + v_t_mag(:) = solve_fragment_tangential_velocities(v_t_mag_input=v_t_mag(7:nfrag), lerr=lerr) + ke_offset = objective_function%lastval + end if + lerr = (ke_offset > 0.0_DP) + if (lerr) then + lreduce_fragment_number = .false. + ! No solution exists for this case, so we need to indicate that this should be a merge + ! This may happen due to setting the tangential velocities too high when setting the angular momentum constraint + v_frag(:,:) = 0.0_DP + do i = 1, nfrag + vb_frag(:, i) = vcom(:) + end do + else + ! Shift the radial velocity vectors to align with the center of mass of the collisional system (the origin) + vb_frag(:,1:nfrag) = vmag_to_vb(v_r_mag(1:nfrag), v_r_unit(:,1:nfrag), v_t_mag(1:nfrag), v_t_unit(:,1:nfrag), m_frag(1:nfrag), vcom(:)) + do i = 1, nfrag + v_frag(:, i) = vb_frag(:, i) - vcom(:) + end do + end if + + return + + end subroutine set_fragment_tangential_velocities + + function tangential_objective_function(v_t_mag_input, lerr) result(ke_offset) + !! Author: David A. Minton + !! + !! Objective function for evaluating how close our fragment velocities get to minimizing KE error from our required value + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: v_t_mag_input !! Unknown tangential component of velocity vector set previously by angular momentum constraint + logical, intent(out) :: lerr !! Error flag + ! Result + real(DP) :: ke_offset + ! Internals + integer(I4B) :: i + real(DP), dimension(:,:), allocatable :: v_shift + real(DP), dimension(:), allocatable :: v_t_mag_output + + lerr = .false. + allocate(v_t_mag_output, mold=v_t_mag) + v_t_mag_output(:) = solve_fragment_tangential_velocities(v_t_mag_input=v_t_mag_input(:), lerr=lerr) + if (lerr) then + v_t_mag_output(1:6) = 0.0_DP + if (nfrag > 6) v_t_mag_output(7:nfrag) = v_t_mag_input(:) + end if + + allocate(v_shift, mold=vb_frag) + v_shift(:,:) = vmag_to_vb(v_r_mag, v_r_unit, v_t_mag_output, v_t_unit, m_frag, vcom) + + ke_frag = 0.0_DP + do i = 1, nfrag + ke_frag = ke_frag + 0.5_DP * m_frag(i) * dot_product(v_shift(:, i), v_shift(:, i)) + end do + ke_offset = ke_frag - ke_target + lerr = (ke_offset > 0.0_DP) + return + end function tangential_objective_function + + function solve_fragment_tangential_velocities(lerr, v_t_mag_input) result(v_t_mag_output) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Adjusts the positions, velocities, and spins of a collection of fragments such that they conserve angular momentum + implicit none + ! Arguments + logical, intent(out) :: lerr !! Error flag + real(DP), dimension(:), optional, intent(in) :: v_t_mag_input !! Unknown tangential velocities for fragments 7:nfrag + ! Internals + integer(I4B) :: i + ! Result + real(DP), dimension(:), allocatable :: v_t_mag_output + + real(DP), dimension(2 * NDIM, 2 * NDIM) :: A ! LHS of linear equation used to solve for momentum constraint in Gauss elimination code + 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, v + + v_frag(:,:) = 0.0_DP + lerr = .false. - L_orb_mag = norm2(L_frag_orb(:)) ! 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 ! The second 3 are that the sum of the angular momentum of the fragments is conserved from the pre-impact state @@ -618,30 +610,23 @@ subroutine set_fragment_tangential_velocities() A(1:3, i) = m_frag(i) * v_t_unit(:, i) call util_crossproduct(v_r_unit(:, i), v_t_unit(:, i), L(:)) A(4:6, i) = m_frag(i) * rmag(i) * L(:) - else !For the remining bodies, distribute the angular momentum equally amongs them - v_t_mag(i) = L_orb_mag / (m_frag(i) * rmag(i) * nfrag) - v_frag(:, i) = v_t_mag(i) * v_t_unit(:, i) - L_lin_others(:) = L_lin_others(:) + m_frag(i) * v_frag(:, i) - call util_crossproduct(x_frag(:, i), v_frag(:, i), L(:)) + else if (present(v_t_mag_input)) then + v(:) = v_t_mag_input(i - 6) * v_t_unit(:, i) + L_lin_others(:) = L_lin_others(:) + m_frag(i) * v(:) + call util_crossproduct(x_frag(:, i), v(:), L(:)) L_orb_others(:) = L_orb_others(:) + m_frag(i) * L(:) end if end do b(1:3) = -L_lin_others(:) b(4:6) = L_frag_orb(:) - L_orb_others(:) - v_t_mag(1:6) = util_solve_linear_system(A, b, 6, lmerge) - if (lmerge) return - do i = 1, 6 - v_frag(:, i) = v_t_mag(i) * v_t_unit(:, i) - end do - - do i = 1, nfrag - vb_frag(:,i) = v_frag(:,i) + vcom(:) - end do + allocate(v_t_mag_output, mold=v_t_mag) + v_t_mag_output(1:6) = util_solve_linear_system(A, b, 6, lerr) + if (present(v_t_mag_input)) v_t_mag_output(7:nfrag) = v_t_mag_input(:) return - end subroutine set_fragment_tangential_velocities + end function solve_fragment_tangential_velocities - subroutine set_fragment_radial_velocities(lmerge) + subroutine set_fragment_radial_velocities(lerr) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! @@ -651,27 +636,18 @@ subroutine set_fragment_radial_velocities(lmerge) !! It takes in the initial "guess" of velocities and solve for the a scaling factor applied to the radial component wrt the !! center of mass frame needed to correct the kinetic energy of the fragments in the system barycenter frame to match that of !! the target kinetic energy required to satisfy the constraints. - use symba_frag_pos_lambda_implementation implicit none ! Arguments - logical, intent(out) :: lmerge + logical, intent(out) :: lerr ! Internals - real(DP), parameter :: TOL = 1e-9_DP - real(DP), dimension(:), allocatable :: vflat - logical :: lerr + real(DP), parameter :: TOL = 1e-6_DP integer(I4B) :: i real(DP), dimension(:), allocatable :: v_r_initial, v_r_sigma real(DP), dimension(:,:), allocatable :: v_r + type(lambda_obj) :: objective_function ! Set the "target" ke_orbit_after (the value of the orbital kinetic energy that the fragments ought to have) - if ((ke_target < 0.0_DP) .or. (ke_offset > 0.0_DP)) then - !if (ke_target < 0.0_DP) write(*,*) 'Negative ke_target: ',ke_target - !if (ke_offset > 0.0_DP) write(*,*) 'Positive ke_offset: ',ke_offset - lmerge = .true. - return - end if - allocate(v_r_initial, source=v_r_mag) ! Initialize radial velocity magnitudes with a random value that is approximately 10% of that found by distributing the kinetic energy equally allocate(v_r_sigma, source=v_r_mag) @@ -681,18 +657,17 @@ subroutine set_fragment_radial_velocities(lmerge) ! Initialize the lambda function using a structure constructor that calls the init method ! Minimize the ke objective function using the BFGS optimizer - v_r_mag(1:nfrag) = util_minimize_bfgs(ke_constraint(ke_objective_function, v_r_unit(:,1:nfrag), v_t_mag(1:nfrag), v_t_unit(:,1:nfrag), m_frag(1:nfrag), vcom(:), ke_target), & - nfrag, v_r_initial(1:nfrag), TOL, lerr) + lerr = .false. + objective_function = lambda_obj(radial_objective_function) + v_r_mag = util_minimize_bfgs(objective_function, nfrag, v_r_initial, TOL, lerr) if (lerr) then ! No solution exists for this case, so we need to indicate that this should be a merge ! This may happen due to setting the tangential velocities too high when setting the angular momentum constraint - lmerge = .true. v_frag(:,:) = 0.0_DP do i = 1, nfrag vb_frag(:, i) = vcom(:) end do else - lmerge = .false. ! Shift the radial velocity vectors to align with the center of mass of the collisional system (the origin) vb_frag(:,1:nfrag) = vmag_to_vb(v_r_mag(1:nfrag), v_r_unit(:,1:nfrag), v_t_mag(1:nfrag), v_t_unit(:,1:nfrag), m_frag(1:nfrag), vcom(:)) do i = 1, nfrag @@ -703,18 +678,13 @@ subroutine set_fragment_radial_velocities(lmerge) return end subroutine set_fragment_radial_velocities - function ke_objective_function(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom, ke_target) result(fval) + function radial_objective_function(v_r_mag_input) result(fval) !! Author: David A. Minton !! !! Objective function for evaluating how close our fragment velocities get to minimizing KE error from our required value implicit none ! Arguments - real(DP), dimension(:), intent(in) :: v_r_mag !! Unknown radial component of fragment velocity vector - real(DP), dimension(:), intent(in) :: v_t_mag !! Tangential component of velocity vector set previously by angular momentum constraint - real(DP), dimension(:,:), intent(in) :: v_r_unit, v_t_unit !! Radial and tangential unit vectors for each fragment - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:), intent(in) :: vcom !! Barycentric velocity of collisional system center of mass - real(DP), intent(in) :: ke_target !! Target kinetic energy + real(DP), dimension(:), intent(in) :: v_r_mag_input !! Unknown radial component of fragment velocity vector ! Result real(DP) :: fval !! The objective function result, which is the square of the difference between the calculated fragment kinetic energy and our target !! Minimizing this brings us closer to our objective @@ -724,7 +694,7 @@ function ke_objective_function(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_frag, vco real(DP), dimension(:,:), allocatable :: v_shift allocate(v_shift, mold=vb_frag) - v_shift(:,:) = vmag_to_vb(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom) + v_shift(:,:) = vmag_to_vb(v_r_mag_input, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom) fval = -ke_target do i = 1, nfrag fval = fval + 0.5_DP * m_frag(i) * dot_product(v_shift(:, i), v_shift(:, i)) @@ -734,7 +704,7 @@ function ke_objective_function(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_frag, vco return - end function ke_objective_function + end function radial_objective_function function vmag_to_vb(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom) result(vb) !! Author: David A. Minton @@ -775,7 +745,7 @@ subroutine restructure_failed_fragments() implicit none integer(I4B) :: i - if (ke_offset < 0.0_DP) then + if (lreduce_fragment_number) then if (iflip == 1) then iflip = 2 else @@ -787,7 +757,7 @@ subroutine restructure_failed_fragments() rad_frag(nfrag) = 0.0_DP nfrag = nfrag - 1 end if - r_max_start = r_max_start + 1.0_DP ! The larger lever arm can help if the problem is in the angular momentum step + r_max_start = r_max_start + 1.00_DP ! The larger lever arm can help if the problem is in the angular momentum step end subroutine restructure_failed_fragments diff --git a/src/user/user_dump_param.f90 b/src/user/user_dump_param.f90 index e464a149c..a1505b978 100644 --- a/src/user/user_dump_param.f90 +++ b/src/user/user_dump_param.f90 @@ -39,6 +39,7 @@ module subroutine user_dump_param(param,t) !! due to compiler incompatabilities !write(LUN,'(DT)') param_dump call param_dump%udio_writer(LUN, iotype="none",v_list=(/0/),iostat=ierr,iomsg=error_message) + call random_seed(put = param_dump%seed) if (idx == 2) then idx = 1 diff --git a/src/user/user_udio_reader.f90 b/src/user/user_udio_reader.f90 index 9814b98fe..027eadefb 100644 --- a/src/user/user_udio_reader.f90 +++ b/src/user/user_udio_reader.f90 @@ -161,18 +161,22 @@ module subroutine user_udio_reader(param, unit, iotype, v_list, iostat, iomsg) ! Because the number of seeds can vary between compilers/systems, we need to make sure we can handle cases in which the input file has a different ! number of seeds than the current system. If the number of seeds in the file is smaller than required, we will use them as a source to fill in the missing elements. ! If the number of seeds in the file is larger than required, we will truncate the seed array. - if (nseeds_from_file >= nseeds) then + if (nseeds_from_file > nseeds) then + nseeds = nseeds_from_file + deallocate(param%seed) + allocate(param%seed(nseeds)) do i = 1, nseeds ifirst = ilast + 1 param_value = user_get_token(line, ifirst, ilast, iostat) read(param_value, *) param%seed(i) end do - else if (nseeds_from_file < nseeds) then ! Seed array in file is too small + else ! Seed array in file is too small do i = 1, nseeds_from_file ifirst = ilast + 1 param_value = user_get_token(line, ifirst, ilast, iostat) read(param_value, *) param%seed(i) end do + param%seed(nseeds_from_file+1:nseeds) = [(param%seed(1) - param%seed(nseeds_from_file) + i, i=nseeds_from_file+1, nseeds)] end if seed_set = .true. case default diff --git a/src/util/util_minimize_bfgs.f90 b/src/util/util_minimize_bfgs.f90 index 45000b5f3..846cee0a8 100644 --- a/src/util/util_minimize_bfgs.f90 +++ b/src/util/util_minimize_bfgs.f90 @@ -20,7 +20,7 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) implicit none ! Arguments integer(I4B), intent(in) :: N - class(lambda_obj), intent(in) :: f + class(lambda_obj), intent(inout) :: f real(DP), dimension(:), intent(in) :: x0 real(DP), intent(in) :: eps logical, intent(out) :: lerr @@ -50,7 +50,11 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) ! Initialize approximate Hessian with the identity matrix (i.e. begin with method of steepest descent) ! Get initial gradient and initialize arrays for updated values of gradient and x H(:,:) = reshape([((0._DP, i=1, j-1), 1._DP, (0._DP, i=j+1, N), j=1, N)], [N,N]) - grad0 = gradf(f, N, x0(:), gradeps) + grad0 = gradf(f, N, x0(:), gradeps, lerr) + if (lerr) then + call ieee_set_status(original_fpe_status) + return + end if grad1(:) = grad0(:) do i = 1, MAXLOOP !check for convergence @@ -70,7 +74,7 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) x1(:) = x1(:) + P(:) ! Calculate new gradient grad0(:) = grad1(:) - grad1 = gradf(f, N, x1, gradeps) + grad1 = gradf(f, N, x1, gradeps, lerr) y(:) = grad1(:) - grad0(:) Py = sum(P(:) * y(:)) ! set up factors for H matrix update @@ -119,7 +123,7 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) contains - function gradf(f, N, x1, dx) result(grad) + function gradf(f, N, x1, dx, lerr) result(grad) !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !! Purpose: Estimates the gradient of a function using a central difference !! approximation @@ -128,21 +132,25 @@ function gradf(f, N, x1, dx) result(grad) !! N : number of variables N !! x1 : x value array !! dx : step size to use when calculating derivatives + !! Outputs: + !! lerr : .true. if an error occurred. Otherwise returns .false. !! Returns !! grad : N sized array containing estimated gradient of f at x1 !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! Arguments integer(I4B), intent(in) :: N - class(lambda_obj), intent(in) :: f + class(lambda_obj), intent(inout) :: f real(DP), dimension(:), intent(in) :: x1 real(DP), intent(in) :: dx + logical, intent(out) :: lerr ! Result real(DP), dimension(N) :: grad ! Internals integer(I4B) :: i, j real(DP), dimension(N) :: xp, xm real(DP) :: fp, fm + logical :: lerrp, lerrm do i = 1, N do j = 1, N @@ -154,9 +162,20 @@ function gradf(f, N, x1, dx) result(grad) xm(j) = x1(j) end if end do - fp = f%eval(xp) - fm = f%eval(xm) + select type (f) + class is (lambda_obj_err) + fp = f%eval(xp) + lerrp = f%lerr + fm = f%eval(xm) + lerrm = f%lerr + lerr = lerrp .or. lerrm + class is (lambda_obj) + fp = f%eval(xp) + fm = f%eval(xm) + lerr = .false. + end select grad(i) = (fp - fm) / (2 * dx) + if (lerr) return end do return end function gradf @@ -182,7 +201,7 @@ function minimize1D(f, x0, S, N, eps, lerr) result(astar) implicit none ! Arguments integer(I4B), intent(in) :: N - class(lambda_obj), intent(in) :: f + class(lambda_obj), intent(inout) :: f real(DP), dimension(:), intent(in) :: x0, S real(DP), intent(in) :: eps logical, intent(out) :: lerr @@ -234,13 +253,15 @@ function minimize1D(f, x0, S, N, eps, lerr) result(astar) return end function minimize1D - function n2one(f, x0, S, N, a) result(fnew) + function n2one(f, x0, S, N, a, lerr) result(fnew) implicit none ! Arguments integer(I4B), intent(in) :: N - class(lambda_obj), intent(in) :: f + class(lambda_obj), intent(inout) :: f real(DP), dimension(:), intent(in) :: x0, S real(DP), intent(in) :: a + logical, intent(out) :: lerr + ! Return real(DP) :: fnew ! Internals @@ -249,6 +270,12 @@ function n2one(f, x0, S, N, a) result(fnew) xnew(:) = x0(:) + a * S(:) fnew = f%eval(xnew(:)) + select type(f) + class is (lambda_obj_err) + lerr = f%lerr + class is (lambda_obj) + lerr = .false. + end select return end function n2one @@ -272,7 +299,7 @@ subroutine bracket(f, x0, S, N, gam, step, lo, hi, lerr) implicit none ! Arguments integer(I4B), intent(in) :: N - class(lambda_obj), intent(in) :: f + class(lambda_obj), intent(inout) :: f real(DP), dimension(:), intent(in) :: x0, S real(DP), intent(in) :: gam, step real(DP), intent(inout) :: lo @@ -281,20 +308,22 @@ subroutine bracket(f, x0, S, N, gam, step, lo, hi, lerr) ! Internals real(DP) :: a0, a1, a2, a3, da real(DP) :: f0, f1, f2, fcon - integer(I4B) :: i - integer(I4B), parameter :: MAXLOOP = 1000 ! maximum number of loops before method is determined to have failed - real(DP), parameter :: eps = epsilon(lo) ! small number precision to test floating point equality + integer(I4B) :: i, j + integer(I4B), parameter :: MAXLOOP = 10000 ! maximum number of loops before method is determined to have failed + real(DP), parameter :: eps = 2 * epsilon(lo) ! small number precision to test floating point equality real(DP), parameter :: dela = 2.0324_DP ! arbitrary number to test if function is constant ! set up initial bracket points - lerr = .false. a0 = lo da = step a1 = a0 + da a2 = a0 + 2 * da - f0 = n2one(f, x0, S, N, a0) - f1 = n2one(f, x0, S, N, a1) - f2 = n2one(f, x0, S, N, a2) + f0 = n2one(f, x0, S, N, a0, lerr) + if (lerr) return + f1 = n2one(f, x0, S, N, a1, lerr) + if (lerr) return + f2 = n2one(f, x0, S, N, a2, lerr) + if (lerr) return ! loop over bracket method until either min is bracketed method fails do i = 1, MAXLOOP if ((f0 > f1) .and. (f2 > f1)) then ! minimum was found @@ -309,7 +338,7 @@ subroutine bracket(f, x0, S, N, gam, step, lo, hi, lerr) a2 = a3 f0 = f1 f1 = f2 - f2 = n2one(f, x0, S, N, a2) + f2 = n2one(f, x0, S, N, a2, lerr) else if ((f0 < f1) .and. (f1 < f2)) then ! function appears to increase da = da * gam a3 = a0 - da @@ -317,7 +346,7 @@ subroutine bracket(f, x0, S, N, gam, step, lo, hi, lerr) a1 = a0 a0 = a3 f2 = f1 - f0 = n2one(f, x0, S, N, a0) + f0 = n2one(f, x0, S, N, a0, lerr) else if ((f0 < f1) .and. (f1 > f2)) then ! more than one minimum present, so in this case we arbitrarily choose the RHS min da = da * gam a3 = a2 + da @@ -326,25 +355,22 @@ subroutine bracket(f, x0, S, N, gam, step, lo, hi, lerr) a2 = a3 f0 = f1 f1 = f2 - f2 = n2one(f, x0, S, N, a2) + f2 = n2one(f, x0, S, N, a2, lerr) else if ((f0 > f1) .and. (abs(f2 - f1) <= eps)) then ! RHS equal da = da * gam a3 = a2 + da a2 = a3 - f2 = n2one(f, x0, S, N, a2) + f2 = n2one(f, x0, S, N, a2, lerr) else if ((abs(f0 - f1) < eps) .and. (f2 > f1)) then ! LHS equal da = da * gam a3 = a0 - da a0 = a3 - f0 = n2one(f, x0, S, N, a0) + f0 = n2one(f, x0, S, N, a0, lerr) else ! all values equal stops if there is no minimum or takes RHS min if it exists ! test if function itself is constant - fcon = n2one(f, x0, S, N, a2 * dela) !change a2 by an arbitrary number to see if constant - if (abs(f2 - fcon) < eps) then - lerr = .true. - !write(*,*) 'bracket determined function is constant' - return ! function is constant - end if + fcon = n2one(f, x0, S, N, a2 * dela, lerr) !change a2 by an arbitrary number to see if constant + lerr = lerr .or. (abs(f2 - fcon) < eps) + if (lerr) return a3 = a0 + 0.5_DP * (a1 - a0) a0 = a1 a1 = a2 @@ -355,8 +381,9 @@ subroutine bracket(f, x0, S, N, gam, step, lo, hi, lerr) a1 = a2 a2 = a3 f1 = f2 - f2 = n2one(f, x0, S, N, a2) + f2 = n2one(f, x0, S, N, a2, lerr) end if + if (lerr) exit ! An error occurred while evaluating the function end do lerr = .true. return ! no minimum found @@ -384,7 +411,7 @@ subroutine golden(f, x0, S, N, eps, lo, hi, lerr) implicit none ! Arguments integer(I4B), intent(in) :: N - class(lambda_obj), intent(in) :: f + class(lambda_obj), intent(inout) :: f real(DP), dimension(:), intent(in) :: x0, S real(DP), intent(in) :: eps real(DP), intent(inout) :: lo @@ -401,9 +428,10 @@ subroutine golden(f, x0, S, N, eps, lo, hi, lerr) i0 = hi - lo a1 = hi - tau * i0 a2 = lo + tau * i0 - f1 = n2one(f, x0, S, N, a1) - f2 = n2one(f, x0, S, N, a2) - lerr = .false. + f1 = n2one(f, x0, S, N, a1, lerr) + if (lerr) return + f2 = n2one(f, x0, S, N, a2, lerr) + if (lerr) return do i = 1, MAXLOOP if (abs((hi - lo) / i0) <= eps) return ! interval reduced to input amount if (f2 > f1) then @@ -411,14 +439,15 @@ subroutine golden(f, x0, S, N, eps, lo, hi, lerr) a2 = a1 f2 = f1 a1 = hi - tau * (hi - lo) - f1 = n2one(f, x0, S, N, a1) + f1 = n2one(f, x0, S, N, a1, lerr) else lo = a1 a1 = a2 f2 = f1 a2 = hi - (1.0_DP - tau) * (hi - lo) - f2 = n2one(f, x0, S, N, a2) + f2 = n2one(f, x0, S, N, a2, lerr) end if + if (lerr) exit end do lerr = .true. return ! search took too many iterations - no minimum found @@ -443,7 +472,7 @@ subroutine quadfit(f, x0, S, N, eps, lo, hi, lerr) implicit none ! Arguments integer(I4B), intent(in) :: N - class(lambda_obj), intent(in) :: f + class(lambda_obj), intent(inout) :: f real(DP), dimension(:), intent(in) :: x0, S real(DP), intent(in) :: eps real(DP), intent(inout) :: lo @@ -465,9 +494,12 @@ subroutine quadfit(f, x0, S, N, eps, lo, hi, lerr) a3 = hi aold = a1 astar = a2 - f1 = n2one(f, x0, S, N, a1) - f2 = n2one(f, x0, S, N, a2) - f3 = n2one(f, x0, S, N, a3) + f1 = n2one(f, x0, S, N, a1, lerr) + if (lerr) return + f2 = n2one(f, x0, S, N, a2, lerr) + if (lerr) return + f3 = n2one(f, x0, S, N, a3, lerr) + if (lerr) return do i = 1, MAXLOOP ! check to see if convergence is reached and exit errval = abs((astar - aold) / astar) @@ -498,12 +530,12 @@ subroutine quadfit(f, x0, S, N, eps, lo, hi, lerr) call ieee_set_halting_mode(ieee_divide_by_zero, .false.) if (lerr) then !write(*,*) 'quadfit fpe:' - !write(*,*) 'uttil_solve_linear_system failed' + !write(*,*) 'util_solve_linear_system failed' exit end if aold = astar if (soln(2) == soln(3)) then ! Handles the case where they are both 0. 0/0 is an unhandled exception - astar = 0.5_DP + astar = -0.5_DP else astar = -soln(2) / (2 * soln(3)) end if @@ -511,10 +543,15 @@ subroutine quadfit(f, x0, S, N, eps, lo, hi, lerr) if (any(fpe_flag)) then !write(*,*) 'quadfit fpe' !write(*,*) 'soln(2:3): ',soln(2:3) + !write(*,*) 'a1, a2, a3' + !write(*,*) a1, a2, a3 + !write(*,*) 'f1, f2, f3' + !write(*,*) f1, f2, f3 lerr = .true. exit end if - fstar = n2one(f, x0, S, N, astar) + fstar = n2one(f, x0, S, N, astar, lerr) + if (lerr) exit ! keep the three closest a values to astar and discard the fourth d1 = abs(a1 - astar) d2 = abs(a2 - astar) diff --git a/src/util/util_solve_linear_system.f90 b/src/util/util_solve_linear_system.f90 index 59a72b119..73dc49764 100644 --- a/src/util/util_solve_linear_system.f90 +++ b/src/util/util_solve_linear_system.f90 @@ -32,6 +32,7 @@ function util_solve_linear_system_d(A,b,n,lerr) result(x) lerr = any(fpe_flag) if (lerr) then x = 0.0_DP + write(*,*) 'fpe in util_solve_linear_system' else x = real(qx, kind=DP) end if @@ -91,10 +92,11 @@ function solve_wbs(u) result(x) ! solve with backward substitution n = size(u, 1) if (allocated(x)) deallocate(x) if (.not.allocated(x)) allocate(x(n)) - call ieee_set_halting_mode(ieee_divide_by_zero, .false.) - do i = n, 1, -1 - x(i) = (u(i, n + 1) - sum(u(i, i + 1:n) * x(i + 1:n))) / u(i, i) - end do + if (any(abs(u) < tiny(1._DP))) then + x(:) = 0._DP + return + end if + forall (i=n:1:-1) x(i) = (u(i, n + 1) - sum(u(i, i + 1:n) * x(i + 1:n))) / u(i, i) return end function solve_wbs