From 1416a5fedc029647a3ccc28422e9bd9b7189abf8 Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 14 May 2021 23:04:20 -0400 Subject: [PATCH] Scaled the quantities in symba_frag_pos to keep things of O(~1) --- src/modules/module_interfaces.f90 | 18 +-- src/symba/symba_casedisruption.f90 | 4 +- src/symba/symba_casehitandrun.f90 | 4 +- src/symba/symba_casesupercatastrophic.f90 | 4 +- src/symba/symba_frag_pos.f90 | 149 +++++++++++----------- src/util/util_minimize_bfgs.f90 | 17 +-- src/util/util_solve_linear_system.f90 | 6 +- 7 files changed, 98 insertions(+), 104 deletions(-) diff --git a/src/modules/module_interfaces.f90 b/src/modules/module_interfaces.f90 index aa4d9704b..c0905e2d8 100644 --- a/src/modules/module_interfaces.f90 +++ b/src/modules/module_interfaces.f90 @@ -723,8 +723,8 @@ function symba_casedisruption (symba_plA, family, nmergeadd, mergeadd_list, x, v integer(I4B), dimension(:), intent(in) :: family integer(I4B), intent(inout) :: nmergeadd type(symba_merger), intent(inout) :: mergeadd_list - real(DP), dimension(:,:), intent(in) :: x, v, L_spin, Ip - real(DP), dimension(:), intent(in) :: mass, radius, mass_res + real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip + real(DP), dimension(:), intent(inout) :: mass, radius, mass_res type(user_input_parameters),intent(inout) :: param real(DP), intent(in) :: Qloss integer(I4B) :: status @@ -741,8 +741,8 @@ function symba_casehitandrun (symba_plA, family, nmergeadd, mergeadd_list, name, integer(I4B), intent(inout) :: nmergeadd type(symba_merger), intent(inout) :: mergeadd_list integer(I4B), dimension(:), intent(in) :: name - real(DP), dimension(:,:), intent(in) :: x, v, Lspin, Ip - real(DP), dimension(:), intent(in) :: mass, radius, mass_res + real(DP), dimension(:,:), intent(inout) :: x, v, Lspin, Ip + real(DP), dimension(:), intent(inout) :: mass, radius, mass_res type(user_input_parameters),intent(inout) :: param real(DP), intent(in) :: Qloss integer(I4B) :: status @@ -757,8 +757,8 @@ function symba_casemerge (symba_plA, family, nmergeadd, mergeadd_list, x, v, mas integer(I4B), dimension(:), intent(in) :: family integer(I4B), intent(inout) :: nmergeadd type(symba_merger), intent(inout) :: mergeadd_list - real(DP), dimension(:,:), intent(in) :: x, v, lspin, Ip - real(DP), dimension(:), intent(in) :: mass, radius + real(DP), dimension(:,:), intent(inout) :: x, v, lspin, Ip + real(DP), dimension(:), intent(inout) :: mass, radius type(user_input_parameters),intent(inout) :: param integer(I4B) :: status end function symba_casemerge @@ -930,9 +930,9 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi type(symba_pl), intent(inout) :: symba_plA integer(I4B), dimension(:), intent(in) :: family real(DP), intent(in) :: Qloss - real(DP), dimension(:,:), intent(in) :: x, v, L_spin, Ip - real(DP), dimension(:), intent(in) :: mass, radius, m_frag, rad_frag - real(DP), dimension(:,:), intent(in) :: Ip_frag + real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip + real(DP), dimension(:), intent(inout) :: mass, radius, m_frag, rad_frag + real(DP), dimension(:,:), intent(inout) :: Ip_frag real(DP), dimension(:,:), intent(out) :: xb_frag, vb_frag, rot_frag logical, intent(out) :: lmerge end subroutine symba_frag_pos diff --git a/src/symba/symba_casedisruption.f90 b/src/symba/symba_casedisruption.f90 index 0559e1ed2..eb679daa6 100644 --- a/src/symba/symba_casedisruption.f90 +++ b/src/symba/symba_casedisruption.f90 @@ -16,8 +16,8 @@ function symba_casedisruption (symba_plA, family, nmergeadd, mergeadd_list, x, v integer(I4B), dimension(:), intent(in) :: family integer(I4B), intent(inout) :: nmergeadd type(symba_merger), intent(inout) :: mergeadd_list - real(DP), dimension(:), intent(in) :: mass, radius, mass_res - real(DP), dimension(:,:), intent(in) :: x, v, L_spin, Ip + real(DP), dimension(:), intent(inout) :: mass, radius, mass_res + real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip type(user_input_parameters),intent(inout) :: param real(DP), intent(inout) :: Qloss ! Result diff --git a/src/symba/symba_casehitandrun.f90 b/src/symba/symba_casehitandrun.f90 index 08f0682cb..d32c8d730 100644 --- a/src/symba/symba_casehitandrun.f90 +++ b/src/symba/symba_casehitandrun.f90 @@ -16,8 +16,8 @@ function symba_casehitandrun(symba_plA, family, nmergeadd, mergeadd_list, id, x, integer(I4B), intent(inout) :: nmergeadd type(symba_merger), intent(inout) :: mergeadd_list integer(I4B), dimension(:), intent(in) :: id - real(DP), dimension(:,:), intent(in) :: x, v, L_spin, Ip - real(DP), dimension(:), intent(in) :: mass, radius, mass_res + real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip + real(DP), dimension(:), intent(inout) :: mass, radius, mass_res type(user_input_parameters),intent(inout) :: param real(DP), intent(inout) :: Qloss ! Result diff --git a/src/symba/symba_casesupercatastrophic.f90 b/src/symba/symba_casesupercatastrophic.f90 index e14d8564e..6b1b77180 100644 --- a/src/symba/symba_casesupercatastrophic.f90 +++ b/src/symba/symba_casesupercatastrophic.f90 @@ -15,8 +15,8 @@ function symba_casesupercatastrophic (symba_plA, family, nmergeadd, mergeadd_lis integer(I4B), dimension(:), intent(in) :: family integer(I4B), intent(inout) :: nmergeadd type(symba_merger), intent(inout) :: mergeadd_list - real(DP), dimension(:,:), intent(in) :: x, v, L_spin, Ip - real(DP), dimension(:), intent(in) :: mass, radius, mass_res + real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip + real(DP), dimension(:), intent(inout) :: mass, radius, mass_res type(user_input_parameters),intent(inout) :: param real(DP), intent(inout) :: Qloss ! Result diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index 39d9157ec..7c7578675 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -12,24 +12,25 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad use module_interfaces, EXCEPT_THIS_ONE => symba_frag_pos implicit none ! Arguments - type(user_input_parameters), intent(in) :: param - type(symba_pl), intent(inout) :: symba_plA - integer(I4B), dimension(:), intent(in) :: family - real(DP), intent(in) :: Qloss - real(DP), dimension(:,:), intent(in) :: x, v, L_spin, Ip - real(DP), dimension(:), intent(in) :: mass, radius, m_frag, rad_frag - real(DP), dimension(:,:), intent(in) :: Ip_frag - real(DP), dimension(:,:), intent(out) :: xb_frag, vb_frag, rot_frag - logical, intent(out) :: lmerge ! Answers the question: Should this have been a merger instead? + type(user_input_parameters), intent(in) :: param + type(symba_pl), intent(inout) :: symba_plA + integer(I4B), dimension(:), intent(in) :: family + real(DP), intent(in) :: Qloss + real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip + real(DP), dimension(:), intent(inout) :: mass, radius, m_frag, rad_frag + real(DP), dimension(:,:), intent(inout) :: Ip_frag + real(DP), dimension(:,:), intent(inout) :: xb_frag, vb_frag, rot_frag + logical, intent(out) :: lmerge ! Answers the question: Should this have been a merger instead? ! Internals - real(DP), dimension(:,:), allocatable :: x_frag, v_frag ! Fragment positions and velocities in the collision center of mass frame + real(DP) :: mscale, rscale, vscale, Lscale, tscale , Qloss_scaled ! Scale factors that reduce quantities to O(~1) in the collisional system real(DP), dimension(NDIM, 2) :: rot, L_orb integer(I4B) :: i, j, nfrag, fam_size, istart + real(DP), dimension(:,:), allocatable :: x_frag, v_frag real(DP), dimension(NDIM) :: xcom, vcom, Ltot_before, Ltot_after, L_residual, L_spin_frag, L_target, L_frag real(DP) :: mtot, Lmag_before, Lmag_after real(DP) :: Etot_before, Etot_after, ke_before, pe_before real(DP) :: pe_after, ke_spin_before, ke_spin_after, ke_after, ke_family, ke_target, ke_frag - real(DP), dimension(NDIM) :: h, dx + real(DP), dimension(NDIM) :: h, dx, xc, vc, vcom_scaled real(DP) :: rmag logical, dimension(:), allocatable :: lexclude character(len=*), parameter :: fmtlabel = "(A14,10(ES9.2,1X,:))" @@ -40,21 +41,45 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad Mpl => symba_plA%helio%swiftest%mass, Ippl => symba_plA%helio%swiftest%Ip, radpl => symba_plA%helio%swiftest%radius, & rotpl => symba_plA%helio%swiftest%rot, status => symba_plA%helio%swiftest%status, npl => symba_plA%helio%swiftest%nbody, name => symba_plA%helio%swiftest%id) + mtot = 1.0_DP + + ! Set scale factors + mscale = sum(mass(:)) + rscale = norm2(x(:,2) - x(:,1)) + tscale = rscale / norm2(v(:,2) - v(:,1)) + vscale = rscale / tscale + Lscale = mscale * rscale * vscale + + mass = mass / mscale + radius = radius / rscale + x = x / rscale + v = v / vscale + L_spin = L_spin / Lscale + + xb_frag = xb_frag / rscale + vb_frag = vb_frag / rscale + m_frag = m_frag / mscale + rot_frag = rot_frag * tscale + rad_frag = rad_frag / rscale + Qloss_scaled = Qloss / (mscale * vscale**2) + allocate(x_frag, source=xb_frag) allocate(v_frag, source=vb_frag) allocate(v_r, mold=v_frag) allocate(v_t, mold=v_frag) + fam_size = size(family) ! Find the center of mass of the collisional system - mtot = sum(mass(:)) - xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot - vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot + xcom(:) = mass(1) * x(:,1) + mass(2) * x(:,2) + vcom(:) = mass(1) * v(:,1) + mass(2) * v(:,2) L_orb(:, :) = 0.0_DP ! Compute orbital angular momentum of pre-impact system do j = 1, 2 - call utiL_crossproduct(x(:, j) - xcom(:), v(:, j) - vcom(:), x_cross_v(:)) + xc(:) = x(:, j) - xcom(:) + vc(:) = v(:, j) - vcom(:) + call util_crossproduct(xc(:), vc(:), x_cross_v(:)) L_orb(:, j) = mass(j) * x_cross_v(:) end do @@ -72,14 +97,14 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad ! We need the original kinetic energy of just the pre-impact family members in order to balance the energy later ke_family = 0.0_DP do i = 1, fam_size - ke_family = ke_family + Mpl(family(i)) * dot_product(vbpl(:,family(i)), vbpl(:,family(i))) ! + ke_family = ke_family + Mpl(family(i)) / mscale * dot_product(vbpl(:,family(i)), vbpl(:,family(i))) / vscale**2 lexclude(family(i)) = .true. ! For all subsequent energy calculations the pre-impact family members will be replaced by the fragments end do ke_family = 0.5_DP * ke_family nfrag = size(m_frag) ! Initialize positions and velocities of fragments that conserve angular momentum - call symba_frag_pos_initialize_fragments(nfrag, xcom, vcom, x, v, L_orb, L_spin, mass, radius, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag, lmerge) + call symba_frag_pos_initialize_fragments(x, v, L_orb, L_spin, mass, radius, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag, lmerge) if (lmerge) return ! Energy calculation requires the fragments to be in the system barcyentric frame, so we need to temporarily shift them @@ -109,12 +134,12 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad write(*, "(' ---------------------------------------------------------------------------')") write(*, "(' Second pass to get energy ')") write(*, "(' ---------------------------------------------------------------------------')") - write(*,fmtlabel) ' Q_loss |',-Qloss / abs(Etot_before) + write(*,fmtlabel) ' Q_loss |',-Qloss_scaled / abs(Etot_before) write(*, "(' ---------------------------------------------------------------------------')") ! Set the "target" ke_after (the value of the orbital kinetic energy that the fragments ought to have) - ke_target = ke_family + (ke_spin_before - ke_spin_after) + (pe_before - pe_after) - Qloss + ke_target = ke_family + (ke_spin_before - ke_spin_after) + (pe_before - pe_after) - Qloss_scaled L_target(:) = 0.0_DP do i = 1, nfrag call util_crossproduct(x_frag(:, i), v_frag(:, i), L_frag(:)) @@ -155,40 +180,25 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad norm2(Ltot_after - Ltot_before) / Lmag_before lmerge = lmerge .or. ((Etot_after - Etot_before) / abs(Etot_before) > 0._DP) - !if (.not.lmerge) then - ! L_residual(:) = Ltot_before(:) - Ltot_after(:) - ! L_spin_frag(:) = L_residual(:) / nfrag - ! do i = 1, nfrag - ! rot_frag(:,i) = rot_frag(:,i) + L_spin_frag(:) / (Ip_frag(:, i) * m_frag(i) * rad_frag(i)**2) - ! end do - !end if - -! call symba_frag_pos_energy_calc(npl, symba_plA, lexclude, ke_after, ke_spin_after, pe_after, Ltot_after,& -! nfrag=nfrag, Ip_frag=Ip_frag, m_frag=m_frag, rad_frag=rad_frag, xb_frag=xb_frag, vb_frag=vb_frag, rot_frag=rot_frag) -! Etot_after = ke_after + ke_spin_after + pe_after -! L_residual(:) = Ltot_before(:) - Ltot_after(:) -! Lmag_after = norm2(Ltot_after(:)) -! -! write(*, "(' ---------------------------------------------------------------------------')") -! write(*, "(' Third pass for correcting any residual angular momentum ')") -! write(*, "(' ---------------------------------------------------------------------------')") -! !write(*, "(' | T_orb T_spin T pe Etot Ltot')") -! !write(*, "(' ---------------------------------------------------------------------------')") -! write(*,fmtlabel) ' change |',(ke_after - ke_before) / abs(Etot_before), & -! (ke_spin_after - ke_spin_before)/ abs(Etot_before), & -! (ke_after + ke_spin_after - ke_before - ke_spin_before)/ abs(Etot_before), & -! (pe_after - pe_before) / abs(Etot_before), & -! (Etot_after - Etot_before) / abs(Etot_before), & -! norm2(Ltot_after - Ltot_before) / Lmag_before -! write(*, "(' ---------------------------------------------------------------------------')") - !**************************************************************************************************************** + + mass = mass * mscale + radius = radius * rscale + x = x * rscale + v = v * vscale + L_spin = L_spin * Lscale + + xb_frag = xb_frag * rscale + vb_frag = vb_frag * rscale + m_frag = m_frag * mscale + rot_frag = rot_frag / tscale + rad_frag = rad_frag * rscale end associate return contains - subroutine symba_frag_pos_initialize_fragments(nfrag, xcom, vcom, x, v, L_orb, L_spin, mass, radius, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag, lmerge) + subroutine symba_frag_pos_initialize_fragments(x, v, L_orb, L_spin, mass, radius, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag, lmerge) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Initializes the orbits of the fragments around the center of mass. The fragments are initially placed on a plane defined by the @@ -196,8 +206,6 @@ subroutine symba_frag_pos_initialize_fragments(nfrag, xcom, vcom, x, v, L_orb, L !! The initial positions do not conserve energy or momentum, so these need to be adjusted later. implicit none ! Arguments - integer(I4B), intent(in) :: nfrag !! Number of collisional fragments - real(DP), dimension(:), intent(in) :: xcom, vcom !! Center of mass position and velocity in the system barycenter frame real(DP), dimension(:,:), intent(in) :: x, v, L_orb, L_spin !! Pre-impact angular momentum vectors real(DP), dimension(:), intent(in) :: mass, radius !! Pre-impact masses and radii real(DP), dimension(:), intent(in) :: m_frag, rad_frag !! Fragment masses and radii @@ -205,13 +213,12 @@ subroutine symba_frag_pos_initialize_fragments(nfrag, xcom, vcom, x, v, L_orb, L real(DP), dimension(:,:), intent(out) :: x_frag, v_frag, rot_frag !! Fragment position, velocities, and spin states logical, intent(out) :: lmerge ! Internals - real(DP) :: mtot, theta, v_frag_norm, r_frag_norm, v_col_norm, r_col_norm + real(DP) :: theta, v_frag_norm, r_frag_norm, v_col_norm, r_col_norm real(DP), dimension(NDIM) :: Ltot, delta_r, delta_v real(DP), dimension(NDIM) :: x_col_unit, y_col_unit, z_col_unit integer(I4B) :: i ! Now create the fragment distribution - mtot = sum(mass(:)) ! Compute orbital angular momentum of pre-impact system. This will be the normal vector to the collision fragment plane Ltot = L_spin(:,1) + L_spin(:,2) + L_orb(:, 1) + L_orb(:, 2) @@ -243,19 +250,17 @@ subroutine symba_frag_pos_initialize_fragments(nfrag, xcom, vcom, x, v, L_orb, L call symba_frag_pos_com_adjust(xcom, m_frag, x_frag) v_frag(:,:) = 0._DP - call symba_frag_pos_ang_mtm(nfrag, xcom, vcom, L_orb, L_spin, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag, lmerge) + call symba_frag_pos_ang_mtm(L_orb, L_spin, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag, lmerge) return end subroutine symba_frag_pos_initialize_fragments - subroutine symba_frag_pos_ang_mtm(nfrag, xcom, vcom, L_orb, L_spin, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag, lmerge) + subroutine symba_frag_pos_ang_mtm(L_orb, L_spin, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag, lmerge) !! 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 - integer(I4B), intent(in) :: nfrag !! Number of collisional fragments - real(DP), dimension(:), intent(in) :: xcom, vcom !! Center of mass position and velocity in the system barycenter frame real(DP), dimension(:,:), intent(in) :: L_orb, L_spin !! Pre-impact position, velocity, and spin states, Ip_frag real(DP), dimension(:), intent(in) :: m_frag, rad_frag !! Fragment masses and radii real(DP), dimension(:,:), intent(in) :: Ip_frag, x_frag !! Fragment prinicpal moments of inertia and position vectors @@ -345,16 +350,11 @@ subroutine symba_frag_pos_kinetic_energy(m_frag, ke_target, L_target, x_frag, v_ logical, intent(out) :: lmerge ! Internals - real(DP) :: mtot !! Total mass of fragments - integer(I4B) :: i, nfrag, neval - real(DP), parameter :: TOL = 1e-9_DP + integer(I4B) :: i, neval + real(DP), parameter :: TOL = 1e-12_DP real(DP), dimension(:), allocatable :: vflat logical :: lerr - - nfrag = size(m_frag) - mtot = sum(m_frag) - ! Initialize the lambda function using a structure constructor that calls the init method ! Minimize error using the BFGS optimizer vflat = util_minimize_bfgs(symba_frag_lambda(lambda=symba_frag_pos_ke_objective_function, v_frag=v_frag, x_frag=x_frag, m_frag=m_frag, L_target=L_target, ke_target=ke_target), & @@ -430,11 +430,8 @@ subroutine symba_frag_pos_com_adjust(vec_com, m_frag, vec_frag) ! Internals real(DP), dimension(NDIM) :: mvec_frag, COM_offset - real(DP) :: mtot - integer(I4B) :: i, nfrag + integer(I4B) :: i - mtot = sum(m_frag(:)) - nfrag = size(m_frag(:)) mvec_frag(:) = 0.0_DP do i = 1, nfrag @@ -492,14 +489,14 @@ subroutine symba_frag_pos_energy_calc(npl, symba_plA, lexclude, ke_orbit, ke_spi ! Copy over old data symba_plwksp%helio%swiftest%id(1:npl) = symba_plA%helio%swiftest%id(1:npl) symba_plwksp%helio%swiftest%status(1:npl) = symba_plA%helio%swiftest%status(1:npl) - symba_plwksp%helio%swiftest%mass(1:npl) = symba_plA%helio%swiftest%mass(1:npl) - symba_plwksp%helio%swiftest%radius(1:npl) = symba_plA%helio%swiftest%radius(1:npl) - symba_plwksp%helio%swiftest%xh(:,1:npl) = symba_plA%helio%swiftest%xh(:,1:npl) - symba_plwksp%helio%swiftest%vh(:,1:npl) = symba_plA%helio%swiftest%vh(:,1:npl) - symba_plwksp%helio%swiftest%rhill(1:npl) = symba_plA%helio%swiftest%rhill(1:npl) - symba_plwksp%helio%swiftest%xb(:,1:npl) = symba_plA%helio%swiftest%xb(:,1:npl) - symba_plwksp%helio%swiftest%vb(:,1:npl) = symba_plA%helio%swiftest%vb(:,1:npl) - symba_plwksp%helio%swiftest%rot(:,1:npl) = symba_plA%helio%swiftest%rot(:,1:npl) + symba_plwksp%helio%swiftest%mass(1:npl) = symba_plA%helio%swiftest%mass(1:npl) / mscale + symba_plwksp%helio%swiftest%radius(1:npl) = symba_plA%helio%swiftest%radius(1:npl) / rscale + symba_plwksp%helio%swiftest%xh(:,1:npl) = symba_plA%helio%swiftest%xh(:,1:npl) / rscale + symba_plwksp%helio%swiftest%vh(:,1:npl) = symba_plA%helio%swiftest%vh(:,1:npl) / vscale + symba_plwksp%helio%swiftest%rhill(1:npl) = symba_plA%helio%swiftest%rhill(1:npl) / rscale + symba_plwksp%helio%swiftest%xb(:,1:npl) = symba_plA%helio%swiftest%xb(:,1:npl) / rscale + symba_plwksp%helio%swiftest%vb(:,1:npl) = symba_plA%helio%swiftest%vb(:,1:npl) / vscale + symba_plwksp%helio%swiftest%rot(:,1:npl) = symba_plA%helio%swiftest%rot(:,1:npl) * tscale symba_plwksp%helio%swiftest%Ip(:,1:npl) = symba_plA%helio%swiftest%Ip(:,1:npl) if (present(nfrag)) then ! Append the fragments if they are included @@ -532,9 +529,5 @@ subroutine symba_frag_pos_energy_calc(npl, symba_plA, lexclude, ke_orbit, ke_spi return end subroutine symba_frag_pos_energy_calc - - - - end subroutine symba_frag_pos diff --git a/src/util/util_minimize_bfgs.f90 b/src/util/util_minimize_bfgs.f90 index 8971cf6e3..dc11dc5f3 100644 --- a/src/util/util_minimize_bfgs.f90 +++ b/src/util/util_minimize_bfgs.f90 @@ -27,7 +27,8 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) real(DP), dimension(:), allocatable :: x1 ! Internals integer(I4B) :: i, j, k, l, conv, num, fnum - integer(I4B), parameter :: MAXLOOP = 2000 !! Maximum number of loops before method is determined to have failed + integer(I4B), parameter :: MAXLOOP = 100 !! Maximum number of loops before method is determined to have failed + real(DP), parameter :: gradeps = 1e-5_DP !! Tolerance for gradient calculations real(DP), dimension(N) :: S !! Direction vectors real(DP), dimension(N) :: Snorm !! normalized direction real(DP), dimension(N,N) :: H !! Approximated inverse Hessian matrix @@ -43,7 +44,7 @@ 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) H(:,:) = reshape([((0._DP, i=1, j-1), 1._DP, (0._DP, i=j+1, N), j=1, N)], [N,N]) ! Get initial gradient and initialize arrays for updated values of gradient and x - fnum = fnum + gradf(f, N, x0(:), grad0, eps) + fnum = fnum + gradf(f, N, x0(:), grad0, gradeps) allocate(x1, source=x0) grad1(:) = grad0(:) do i = 1, MAXLOOP @@ -69,7 +70,7 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) x1(:) = x1(:) + P(:) ! Calculate new gradient grad0(:) = grad1(:) - fnum = fnum + gradf(f, N, x1, grad1, eps) + fnum = fnum + gradf(f, N, x1, grad1, gradeps) y(:) = grad1(:) - grad0(:) Py = sum(P(:) * y(:)) ! set up factors for H matrix update @@ -495,11 +496,11 @@ function quadfit(f, x0, S, N, lo, hi, eps) result(fnum) ! Solve system of equations soln(:) = util_solve_linear_system(lhs, rhs, 3, lerr) if (lerr) then - write(*,*) "Could not solve polynomial on loop ", i - write(*,'("a1 = ",f9.6," f1 = ",f9.6)') a1, f1 - write(*,'("a2 = ",f9.6," f2 = ",f9.6)') a2, f2 - write(*,'("a3 = ",f9.6," f3 = ",f9.6)') a3, f3 - write(*,'("aold = ",f7.4)') aold + !write(*,*) "Could not solve polynomial on loop ", i + !write(*,'("a1 = ",f9.6," f1 = ",f9.6)') a1, f1 + !write(*,'("a2 = ",f9.6," f2 = ",f9.6)') a2, f2 + !write(*,'("a3 = ",f9.6," f3 = ",f9.6)') a3, f3 + !write(*,'("aold = ",f7.4)') aold fnum = 0 return end if diff --git a/src/util/util_solve_linear_system.f90 b/src/util/util_solve_linear_system.f90 index b52f08796..3933d1f68 100644 --- a/src/util/util_solve_linear_system.f90 +++ b/src/util/util_solve_linear_system.f90 @@ -38,9 +38,9 @@ function solve_wbs(u, lerr) result(x) ! solve with backward substitution integer(I4B) :: i,n real(DP), parameter :: epsilon = 10 * tiny(1._DP) - if (allocated(x)) deallocate(x) - allocate(x(n)) - + n = size(u, 1) + if (allocated(x) .and. (size(x) /= n)) deallocate(x) + if (.not.allocated(x)) allocate(x(n)) lerr = any(abs(u(:,:)) < epsilon) if (lerr) then x(:) = 0._DP