From 6706e16e7662cf8802d28af1a589f70247106809 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 6 Aug 2021 17:05:01 -0400 Subject: [PATCH 1/8] Added missing G term to potential energy calculations --- src/fragmentation/fragmentation.f90 | 40 +++++++++++++-------------- src/util/util_get_energy_momentum.f90 | 4 +-- 2 files changed, 22 insertions(+), 22 deletions(-) diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index ad28edf6c..e1b69dd44 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -525,21 +525,21 @@ subroutine set_fragment_tan_vel(lerr) call calculate_system_energy(linclude_fragments=.true.) ke_frag_budget = -dEtot - Qloss - !write(*,*) '***************************************************' - !write(*,*) 'Original dis : ',norm2(x(:,2) - x(:,1)) - !write(*,*) 'r_max : ',r_max - !write(*,*) 'f_spin : ',f_spin - !write(*,*) '***************************************************' - !write(*,*) 'Energy balance so far: ' - !write(*,*) 'ke_frag_budget : ',ke_frag_budget - !write(*,*) 'ke_orbit_before: ',ke_orbit_before - !write(*,*) 'ke_orbit_after : ',ke_orbit_after - !write(*,*) 'ke_spin_before : ',ke_spin_before - !write(*,*) 'ke_spin_after : ',ke_spin_after - !write(*,*) 'pe_before : ',pe_before - !write(*,*) 'pe_after : ',pe_after - !write(*,*) 'Qloss : ',Qloss - !write(*,*) '***************************************************' + write(*,*) '***************************************************' + write(*,*) 'Original dis : ',norm2(x(:,2) - x(:,1)) + write(*,*) 'r_max : ',r_max + write(*,*) 'f_spin : ',f_spin + write(*,*) '***************************************************' + write(*,*) 'Energy balance so far: ' + write(*,*) 'ke_frag_budget : ',ke_frag_budget + write(*,*) 'ke_orbit_before: ',ke_orbit_before + write(*,*) 'ke_orbit_after : ',ke_orbit_after + write(*,*) 'ke_spin_before : ',ke_spin_before + write(*,*) 'ke_spin_after : ',ke_spin_after + write(*,*) 'pe_before : ',pe_before + write(*,*) 'pe_after : ',pe_after + write(*,*) 'Qloss : ',Qloss + write(*,*) '***************************************************' if (ke_frag_budget < 0.0_DP) then write(*,*) 'Negative ke_frag_budget: ',ke_frag_budget r_max_start = r_max_start / 2 @@ -604,11 +604,11 @@ subroutine set_fragment_tan_vel(lerr) ! If we are over the energy budget, flag this as a failure so we can try again lerr = (ke_radial < 0.0_DP) - !write(*,*) 'Tangential' - !write(*,*) 'ke_frag_budget: ',ke_frag_budget - !write(*,*) 'ke_frag_orbit : ',ke_frag_orbit - !write(*,*) 'ke_frag_spin : ',ke_frag_spin - !write(*,*) 'ke_radial : ',ke_radial + write(*,*) 'Tangential' + write(*,*) 'ke_frag_budget: ',ke_frag_budget + write(*,*) 'ke_frag_orbit : ',ke_frag_orbit + write(*,*) 'ke_frag_spin : ',ke_frag_spin + write(*,*) 'ke_radial : ',ke_radial return diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index fd331bcc3..4ec9c6c36 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -76,7 +76,7 @@ module subroutine util_get_energy_momentum_system(self, param) !$omp simd do i = 1, npl associate(px => pl%xh(1,i), py => pl%xh(2,i), pz => pl%xh(3,i)) - pecb(i) = -cb%mass * pl%mass(i) / sqrt(px**2 + py**2 + pz**2) + pecb(i) = -param%GU * cb%mass * pl%mass(i) / sqrt(px**2 + py**2 + pz**2) end associate end do @@ -100,7 +100,7 @@ module subroutine util_get_energy_momentum_system(self, param) irh(i) = 1.0_DP / norm2(pl%xh(:,i)) end do call obl_pot(npl, cb%mass, pl%mass, cb%j2rp2, cb%j4rp4, pl%xh, irh, oblpot) - system%pe = system%pe + oblpot + system%pe = system%pe + param%GU * oblpot end if system%Lorbit(1) = sum(Lplorbitx(1:npl), lstatus(1:npl)) From a72bdc548c22ca4f3cb287f5090e1751226e6fdc Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 6 Aug 2021 18:30:00 -0400 Subject: [PATCH 2/8] Moved G out from sqrt --- src/fragmentation/fragmentation.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index e1b69dd44..270f48a7a 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -167,7 +167,7 @@ subroutine set_scale_factors() ! Set scale factors Escale = 0.5_DP * (mass(1) * dot_product(v(:,1), v(:,1)) + mass(2) * dot_product(v(:,2), v(:,2))) rscale = sum(radius(:)) - mscale = sqrt(Escale * rscale / param%GU) + mscale = sqrt(Escale * rscale) / param%GU vscale = sqrt(Escale / mscale) tscale = rscale / vscale Lscale = mscale * rscale * vscale From e880213553aa62d9aeeeeba9e0fba77247beaca7 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 7 Aug 2021 08:53:08 -0400 Subject: [PATCH 3/8] Improved scaling and unit consistency --- src/fragmentation/fragmentation.f90 | 18 ++++++++++++++++-- src/util/util_get_energy_momentum.f90 | 8 ++++---- 2 files changed, 20 insertions(+), 6 deletions(-) diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index 270f48a7a..688ac1a6e 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -167,7 +167,7 @@ subroutine set_scale_factors() ! Set scale factors Escale = 0.5_DP * (mass(1) * dot_product(v(:,1), v(:,1)) + mass(2) * dot_product(v(:,2), v(:,2))) rscale = sum(radius(:)) - mscale = sqrt(Escale * rscale) / param%GU + mscale = sqrt(Escale * rscale / param%GU) vscale = sqrt(Escale / mscale) tscale = rscale / vscale Lscale = mscale * rscale * vscale @@ -341,11 +341,25 @@ subroutine calculate_system_energy(linclude_fragments) ltmp(1:npl) = .true. call tmpsys%pl%setup(npl_new, param) call tmpsys%pl%fill(pl, ltmp) + tmpsys%pl%mass(1:npl) = tmpsys%pl%mass(1:npl) / mscale + tmpsys%pl%Gmass(1:npl) = tmpsys%pl%Gmass(1:npl) * tscale**2 / rscale**3 + tmpsys%cb%mass = tmpsys%cb%mass / mscale + tmpsys%cb%Gmass = tmpsys%cb%Gmass * tscale**2 / rscale**3 + tmpsys%cb%radius = tmpsys%cb%radius / rscale + tmpsys%pl%radius(1:npl) = tmpsys%pl%radius(1:npl) / rscale + tmpsys%pl%xh(:,1:npl) = tmpsys%pl%xh(:,1:npl) / rscale + tmpsys%pl%vh(:,1:npl) = tmpsys%pl%vh(:,1:npl) / vscale + tmpsys%pl%xb(:,1:npl) = tmpsys%pl%xb(:,1:npl) / rscale + tmpsys%pl%vb(:,1:npl) = tmpsys%pl%vb(:,1:npl) / vscale + tmpsys%pl%rot(:,1:npl) = tmpsys%pl%rot(:,1:npl) * tscale + tmpsys%cb%xb(:) = tmpsys%cb%xb(:) / rscale + tmpsys%cb%vb(:) = tmpsys%cb%vb(:) / vscale + if (linclude_fragments) then ! Append the fragments if they are included ! Energy calculation requires the fragments to be in the system barcyentric frame, s tmpsys%pl%mass(npl+1:npl_new) = m_frag(1:nfrag) - tmpsys%pl%Gmass(npl+1:npl_new) = param%GU * m_frag(1:nfrag) + tmpsys%pl%Gmass(npl+1:npl_new) = m_frag(1:nfrag) * param%GU * rscale**3 * mscale / tscale**2 tmpsys%pl%radius(npl+1:npl_new) = rad_frag(1:nfrag) tmpsys%pl%xb(:,npl+1:npl_new) = xb_frag(:,1:nfrag) tmpsys%pl%vb(:,npl+1:npl_new) = vb_frag(:,1:nfrag) diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index 4ec9c6c36..05c3e2436 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -76,14 +76,14 @@ module subroutine util_get_energy_momentum_system(self, param) !$omp simd do i = 1, npl associate(px => pl%xh(1,i), py => pl%xh(2,i), pz => pl%xh(3,i)) - pecb(i) = -param%GU * cb%mass * pl%mass(i) / sqrt(px**2 + py**2 + pz**2) + pecb(i) = -cb%Gmass * pl%mass(i) / sqrt(px**2 + py**2 + pz**2) end associate end do ! Do the potential energy between pairs of massive bodies do k = 1, pl%nplpl associate(ik => pl%k_plpl(1, k), jk => pl%k_plpl(2, k)) - pepl(k) = -param%GU * pl%mass(ik) * pl%mass(jk) / norm2(pl%xh(:, jk) - pl%xh(:, ik)) + pepl(k) = -pl%Gmass(ik) * pl%mass(jk) / norm2(pl%xh(:, jk) - pl%xh(:, ik)) lstatpl(k) = (lstatus(ik) .and. lstatus(jk)) end associate end do @@ -99,8 +99,8 @@ module subroutine util_get_energy_momentum_system(self, param) do i = 1, npl irh(i) = 1.0_DP / norm2(pl%xh(:,i)) end do - call obl_pot(npl, cb%mass, pl%mass, cb%j2rp2, cb%j4rp4, pl%xh, irh, oblpot) - system%pe = system%pe + param%GU * oblpot + call obl_pot(npl, cb%Gmass, pl%mass, cb%j2rp2, cb%j4rp4, pl%xh, irh, oblpot) + system%pe = system%pe + oblpot end if system%Lorbit(1) = sum(Lplorbitx(1:npl), lstatus(1:npl)) From dcf602a736a72793cb195606609117577f907e2c Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 7 Aug 2021 09:13:33 -0400 Subject: [PATCH 4/8] Deallocate temporary variables during rearray --- src/symba/symba_util.f90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 2ab088d02..028b0678c 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -388,6 +388,10 @@ module subroutine symba_util_rearray_pl(self, system, param) call tmp%setup(0,param) deallocate(tmp) + ! Deallocate any temporary variables + if (allocated(pl%xbeg)) deallocate(pl%xbeg) + if (allocated(pl%xend)) deallocate(pl%xend) + ! Add in any new bodies call pl%append(pl_adds) From 02b692d1e3279a27f849ea7643bb91be6e1c0049 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 7 Aug 2021 12:53:38 -0400 Subject: [PATCH 5/8] Added new rescale method for systems --- Makefile.Defines | 2 +- src/fragmentation/fragmentation.f90 | 18 ++----- src/modules/swiftest_classes.f90 | 73 ++++++++++++++++------------- src/util/util_rescale.f90 | 51 ++++++++++++++++++++ src/util/util_resize.f90 | 1 + 5 files changed, 98 insertions(+), 47 deletions(-) create mode 100644 src/util/util_rescale.f90 diff --git a/Makefile.Defines b/Makefile.Defines index 291f2c604..07126f842 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -71,7 +71,7 @@ FORTRAN = ifort #AR = xiar #FORTRAN = gfortran -#FFLAGS = -ffree-line-length-none $(GDEBUG) $(GMEM) +#FFLAGS = -ffree-line-length-none $(GDEBUG) #$(GMEM) AR = ar # DO NOT include in CFLAGS the "-c" option to compile object only diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index 688ac1a6e..d879eec53 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -314,6 +314,7 @@ subroutine calculate_system_energy(linclude_fragments) logical, dimension(:), allocatable :: ltmp logical :: lk_plpl class(swiftest_nbody_system), allocatable :: tmpsys + class(swiftest_parameters), allocatable :: tmpparam ! Because we're making a copy of symba_pl with the excludes/fragments appended, we need to deallocate the ! big k_plpl array and recreate it when we're done, otherwise we run the risk of blowing up the memory by @@ -334,27 +335,16 @@ subroutine calculate_system_energy(linclude_fragments) npl_new = npl end if call setup_construct_system(tmpsys, param) + call tmpsys%tp%setup(0, param) deallocate(tmpsys%cb) allocate(tmpsys%cb, source=cb) + allocate(tmpparam, source=param) allocate(ltmp(npl_new)) ltmp(:) = .false. ltmp(1:npl) = .true. call tmpsys%pl%setup(npl_new, param) call tmpsys%pl%fill(pl, ltmp) - tmpsys%pl%mass(1:npl) = tmpsys%pl%mass(1:npl) / mscale - tmpsys%pl%Gmass(1:npl) = tmpsys%pl%Gmass(1:npl) * tscale**2 / rscale**3 - tmpsys%cb%mass = tmpsys%cb%mass / mscale - tmpsys%cb%Gmass = tmpsys%cb%Gmass * tscale**2 / rscale**3 - tmpsys%cb%radius = tmpsys%cb%radius / rscale - tmpsys%pl%radius(1:npl) = tmpsys%pl%radius(1:npl) / rscale - tmpsys%pl%xh(:,1:npl) = tmpsys%pl%xh(:,1:npl) / rscale - tmpsys%pl%vh(:,1:npl) = tmpsys%pl%vh(:,1:npl) / vscale - tmpsys%pl%xb(:,1:npl) = tmpsys%pl%xb(:,1:npl) / rscale - tmpsys%pl%vb(:,1:npl) = tmpsys%pl%vb(:,1:npl) / vscale - tmpsys%pl%rot(:,1:npl) = tmpsys%pl%rot(:,1:npl) * tscale - tmpsys%cb%xb(:) = tmpsys%cb%xb(:) / rscale - tmpsys%cb%vb(:) = tmpsys%cb%vb(:) / vscale - + call tmpsys%rescale(tmpparam, mscale, rscale, tscale) if (linclude_fragments) then ! Append the fragments if they are included ! Energy calculation requires the fragments to be in the system barcyentric frame, s diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 7f4fd140e..4d0e98704 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -306,6 +306,7 @@ module swiftest_classes procedure :: step_spin => tides_step_spin_system !! Steps the spins of the massive & central bodies due to tides. procedure :: set_msys => util_set_msys !! Sets the value of msys from the masses of system bodies. procedure :: get_energy_and_momentum => util_get_energy_momentum_system !! Calculates the total system energy and momentum + procedure :: rescale => util_rescale_system !! Rescales the system into a new set of units procedure :: validate_ids => util_valid_id_system !! Validate the numerical ids passed to the system and save the maximum value end type swiftest_nbody_system @@ -984,6 +985,13 @@ end subroutine util_fill_arr_logical end interface interface + module subroutine util_rescale_system(self, param, mscale, dscale, tscale) + implicit none + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters. Returns with new values of the scale vactors and GU + real(DP), intent(in) :: mscale, dscale, tscale !! Scale factors for mass, distance, and time units, respectively. + end subroutine util_rescale_system + module function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) use lambda_function implicit none @@ -1003,38 +1011,6 @@ module subroutine util_peri_tp(self, system, param) end subroutine util_peri_tp end interface - interface util_solve_linear_system - module function util_solve_linear_system_d(A,b,n,lerr) result(x) - implicit none - integer(I4B), intent(in) :: n - real(DP), dimension(:,:), intent(in) :: A - real(DP), dimension(:), intent(in) :: b - logical, intent(out) :: lerr - real(DP), dimension(n) :: x - end function util_solve_linear_system_d - - module function util_solve_linear_system_q(A,b,n,lerr) result(x) - implicit none - integer(I4B), intent(in) :: n - real(QP), dimension(:,:), intent(in) :: A - real(QP), dimension(:), intent(in) :: b - logical, intent(out) :: lerr - real(QP), dimension(n) :: x - end function util_solve_linear_system_q - end interface - - interface - module function util_solve_rkf45(f, y0in, t1, dt0, tol) result(y1) - use lambda_function - implicit none - class(lambda_obj), intent(inout) :: f !! lambda function object that has been initialized to be a function of derivatives. The object will return with components lastarg and lasteval set - real(DP), dimension(:), intent(in) :: y0in !! Initial value at t=0 - real(DP), intent(in) :: t1 !! Final time - real(DP), intent(in) :: dt0 !! Initial step size guess - real(DP), intent(in) :: tol !! Tolerance on solution - real(DP), dimension(:), allocatable :: y1 !! Final result - end function util_solve_rkf45 - end interface interface util_resize module subroutine util_resize_arr_char_string(arr, nnew) @@ -1142,6 +1118,39 @@ module subroutine util_set_rhill_approximate(self,cb) end subroutine util_set_rhill_approximate end interface + interface util_solve_linear_system + module function util_solve_linear_system_d(A,b,n,lerr) result(x) + implicit none + integer(I4B), intent(in) :: n + real(DP), dimension(:,:), intent(in) :: A + real(DP), dimension(:), intent(in) :: b + logical, intent(out) :: lerr + real(DP), dimension(n) :: x + end function util_solve_linear_system_d + + module function util_solve_linear_system_q(A,b,n,lerr) result(x) + implicit none + integer(I4B), intent(in) :: n + real(QP), dimension(:,:), intent(in) :: A + real(QP), dimension(:), intent(in) :: b + logical, intent(out) :: lerr + real(QP), dimension(n) :: x + end function util_solve_linear_system_q + end interface + + interface + module function util_solve_rkf45(f, y0in, t1, dt0, tol) result(y1) + use lambda_function + implicit none + class(lambda_obj), intent(inout) :: f !! lambda function object that has been initialized to be a function of derivatives. The object will return with components lastarg and lasteval set + real(DP), dimension(:), intent(in) :: y0in !! Initial value at t=0 + real(DP), intent(in) :: t1 !! Final time + real(DP), intent(in) :: dt0 !! Initial step size guess + real(DP), intent(in) :: tol !! Tolerance on solution + real(DP), dimension(:), allocatable :: y1 !! Final result + end function util_solve_rkf45 + end interface + interface util_sort module subroutine util_sort_i4b(arr) implicit none diff --git a/src/util/util_rescale.f90 b/src/util/util_rescale.f90 new file mode 100644 index 000000000..061ecf9a5 --- /dev/null +++ b/src/util/util_rescale.f90 @@ -0,0 +1,51 @@ +submodule (swiftest_classes) s_util_rescale + use swiftest +contains + module subroutine util_rescale_system(self, param, mscale, dscale, tscale) + !! author: David A. Minton + !! + !! Rescales an nbody system to a new set of units. Inputs are the multipliers on the mass (mscale), distance (dscale), and time units (tscale). + !! Rescales all united quantities in the system, as well as the mass conversion factors, gravitational constant, and Einstein's constant in the parameter object. + implicit none + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters. Returns with new values of the scale vactors and GU + real(DP), intent(in) :: mscale, dscale, tscale !! Scale factors for mass, distance, and time units, respectively. + ! Internals + real(DP) :: vscale + + param%MU2KG = param%MU2KG * mscale + param%DU2M = param%DU2M * dscale + param%TU2S = param%TU2S * tscale + + ! Calculate the G for the system units + param%GU = GC / (param%DU2M**3 / (param%MU2KG * param%TU2S**2)) + + ! Calculate the inverse speed of light in the system units + param%inv_c2 = einsteinC * param%TU2S / param%DU2M + param%inv_c2 = (param%inv_c2)**(-2) + + vscale = dscale / tscale + + associate(cb => self%cb, pl => self%pl, npl => self%pl%nbody, tp => self%tp, ntp => self%tp%nbody) + + cb%mass = cb%mass / mscale + cb%Gmass = param%GU * cb%mass + cb%radius = cb%radius / dscale + cb%xb(:) = cb%xb(:) / dscale + cb%vb(:) = cb%vb(:) / vscale + pl%mass(1:npl) = pl%mass(1:npl) / mscale + pl%Gmass(1:npl) = param%GU * pl%mass(1:npl) + pl%radius(1:npl) = pl%radius(1:npl) / dscale + pl%xh(:,1:npl) = pl%xh(:,1:npl) / dscale + pl%vh(:,1:npl) = pl%vh(:,1:npl) / vscale + pl%xb(:,1:npl) = pl%xb(:,1:npl) / dscale + pl%vb(:,1:npl) = pl%vb(:,1:npl) / vscale + pl%rot(:,1:npl) = pl%rot(:,1:npl) * tscale + + end associate + + + return + end subroutine util_rescale_system + +end submodule s_util_rescale \ No newline at end of file diff --git a/src/util/util_resize.f90 b/src/util/util_resize.f90 index c6d5aa34f..80d87209c 100644 --- a/src/util/util_resize.f90 +++ b/src/util/util_resize.f90 @@ -1,6 +1,7 @@ submodule (swiftest_classes) s_util_resize use swiftest contains + module subroutine util_resize_arr_char_string(arr, nnew) !! author: David A. Minton !! From fb38873a6bb188c876f9519c314082a8bacd2489 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 7 Aug 2021 13:50:07 -0400 Subject: [PATCH 6/8] Refactored rscale to dscale for consistency --- src/fragmentation/fragmentation.f90 | 38 ++++++++++++++--------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index d879eec53..5d92a1a3e 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -22,7 +22,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, real(DP), intent(inout) :: Qloss !! Energy lost during the collision logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? ! Internals - real(DP) :: mscale, rscale, vscale, tscale, Lscale, Escale ! Scale factors that reduce quantities to O(~1) in the collisional system + real(DP) :: mscale, dscale, 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 @@ -46,6 +46,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, integer(I4B), parameter :: MAXTRY = 3000 integer(I4B), parameter :: TANTRY = 3 logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes + class(swiftest_parameters), allocatable :: tmpparam if (nfrag < NFRAG_MIN) then write(*,*) "symba_frag_pos needs at least ",NFRAG_MIN," fragments, but only ",nfrag," were given." @@ -59,7 +60,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, f_spin = 0.05_DP mscale = 1.0_DP - rscale = 1.0_DP + dscale = 1.0_DP vscale = 1.0_DP tscale = 1.0_DP Lscale = 1.0_DP @@ -166,19 +167,19 @@ subroutine set_scale_factors() ! Set scale factors Escale = 0.5_DP * (mass(1) * dot_product(v(:,1), v(:,1)) + mass(2) * dot_product(v(:,2), v(:,2))) - rscale = sum(radius(:)) - mscale = sqrt(Escale * rscale / param%GU) + dscale = sum(radius(:)) + mscale = sqrt(Escale * dscale / param%GU) vscale = sqrt(Escale / mscale) - tscale = rscale / vscale - Lscale = mscale * rscale * vscale + tscale = dscale / vscale + Lscale = mscale * dscale * vscale - xcom(:) = xcom(:) / rscale + xcom(:) = xcom(:) / dscale vcom(:) = vcom(:) / vscale mtot = mtot / mscale mass = mass / mscale - radius = radius / rscale - x = x / rscale + radius = radius / dscale + x = x / dscale v = v / vscale L_spin = L_spin / Lscale do i = 1, 2 @@ -186,7 +187,7 @@ subroutine set_scale_factors() end do m_frag = m_frag / mscale - rad_frag = rad_frag / rscale + rad_frag = rad_frag / dscale Qloss = Qloss / Escale return @@ -201,13 +202,13 @@ subroutine restore_scale_factors() call ieee_set_halting_mode(IEEE_ALL,.false.) ! Restore scale factors - xcom(:) = xcom(:) * rscale + xcom(:) = xcom(:) * dscale vcom(:) = vcom(:) * vscale mtot = mtot * mscale mass = mass * mscale - radius = radius * rscale - x = x * rscale + radius = radius * dscale + x = x * dscale v = v * vscale L_spin = L_spin * Lscale do i = 1, 2 @@ -215,9 +216,9 @@ subroutine restore_scale_factors() end do m_frag = m_frag * mscale - rad_frag = rad_frag * rscale + rad_frag = rad_frag * dscale rot_frag = rot_frag / tscale - x_frag = x_frag * rscale + x_frag = x_frag * dscale v_frag = v_frag * vscale Qloss = Qloss * Escale @@ -240,7 +241,7 @@ subroutine restore_scale_factors() Lmag_after = Lmag_after * Lscale mscale = 1.0_DP - rscale = 1.0_DP + dscale = 1.0_DP vscale = 1.0_DP tscale = 1.0_DP Lscale = 1.0_DP @@ -314,7 +315,6 @@ subroutine calculate_system_energy(linclude_fragments) logical, dimension(:), allocatable :: ltmp logical :: lk_plpl class(swiftest_nbody_system), allocatable :: tmpsys - class(swiftest_parameters), allocatable :: tmpparam ! Because we're making a copy of symba_pl with the excludes/fragments appended, we need to deallocate the ! big k_plpl array and recreate it when we're done, otherwise we run the risk of blowing up the memory by @@ -344,12 +344,12 @@ subroutine calculate_system_energy(linclude_fragments) ltmp(1:npl) = .true. call tmpsys%pl%setup(npl_new, param) call tmpsys%pl%fill(pl, ltmp) - call tmpsys%rescale(tmpparam, mscale, rscale, tscale) + call tmpsys%rescale(tmpparam, mscale, dscale, tscale) if (linclude_fragments) then ! Append the fragments if they are included ! Energy calculation requires the fragments to be in the system barcyentric frame, s tmpsys%pl%mass(npl+1:npl_new) = m_frag(1:nfrag) - tmpsys%pl%Gmass(npl+1:npl_new) = m_frag(1:nfrag) * param%GU * rscale**3 * mscale / tscale**2 + tmpsys%pl%Gmass(npl+1:npl_new) = m_frag(1:nfrag) * param%GU * dscale**3 * mscale / tscale**2 tmpsys%pl%radius(npl+1:npl_new) = rad_frag(1:nfrag) tmpsys%pl%xb(:,npl+1:npl_new) = xb_frag(:,1:nfrag) tmpsys%pl%vb(:,npl+1:npl_new) = vb_frag(:,1:nfrag) From 6b04d411504aa5423ce99eaad4a767ffdafb56b5 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 8 Aug 2021 16:53:22 -0400 Subject: [PATCH 7/8] Fixed bad indexing on potential energy calculation. Set mass scaling to be total mass in fragmentation --- src/fragmentation/fragmentation.f90 | 50 ++++++++++++++------------- src/util/util_get_energy_momentum.f90 | 6 ++-- 2 files changed, 29 insertions(+), 27 deletions(-) diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index 5d92a1a3e..0b7f2721a 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -168,7 +168,7 @@ subroutine set_scale_factors() ! Set scale factors Escale = 0.5_DP * (mass(1) * dot_product(v(:,1), v(:,1)) + mass(2) * dot_product(v(:,2), v(:,2))) dscale = sum(radius(:)) - mscale = sqrt(Escale * dscale / param%GU) + mscale = mtot vscale = sqrt(Escale / mscale) tscale = dscale / vscale Lscale = mscale * dscale * vscale @@ -338,6 +338,7 @@ subroutine calculate_system_energy(linclude_fragments) call tmpsys%tp%setup(0, param) deallocate(tmpsys%cb) allocate(tmpsys%cb, source=cb) + if (allocated(tmpparam)) deallocate(tmpparam) allocate(tmpparam, source=param) allocate(ltmp(npl_new)) ltmp(:) = .false. @@ -347,13 +348,14 @@ subroutine calculate_system_energy(linclude_fragments) call tmpsys%rescale(tmpparam, mscale, dscale, tscale) if (linclude_fragments) then ! Append the fragments if they are included - ! Energy calculation requires the fragments to be in the system barcyentric frame, s + ! Energy calculation requires the fragments to be in the system barcyentric frame + tmpsys%pl%mass(npl+1:npl_new) = m_frag(1:nfrag) - tmpsys%pl%Gmass(npl+1:npl_new) = m_frag(1:nfrag) * param%GU * dscale**3 * mscale / tscale**2 + tmpsys%pl%Gmass(npl+1:npl_new) = m_frag(1:nfrag) * tmpparam%GU tmpsys%pl%radius(npl+1:npl_new) = rad_frag(1:nfrag) tmpsys%pl%xb(:,npl+1:npl_new) = xb_frag(:,1:nfrag) tmpsys%pl%vb(:,npl+1:npl_new) = vb_frag(:,1:nfrag) - tmpsys%pl%status(npl+1:npl_new) = COLLISION + tmpsys%pl%status(npl+1:npl_new) = ACTIVE if (param%lrotation) then tmpsys%pl%Ip(:,npl+1:npl_new) = Ip_frag(:,1:nfrag) tmpsys%pl%rot(:,npl+1:npl_new) = rot_frag(:,1:nfrag) @@ -529,21 +531,21 @@ subroutine set_fragment_tan_vel(lerr) call calculate_system_energy(linclude_fragments=.true.) ke_frag_budget = -dEtot - Qloss - write(*,*) '***************************************************' - write(*,*) 'Original dis : ',norm2(x(:,2) - x(:,1)) - write(*,*) 'r_max : ',r_max - write(*,*) 'f_spin : ',f_spin - write(*,*) '***************************************************' - write(*,*) 'Energy balance so far: ' - write(*,*) 'ke_frag_budget : ',ke_frag_budget - write(*,*) 'ke_orbit_before: ',ke_orbit_before - write(*,*) 'ke_orbit_after : ',ke_orbit_after - write(*,*) 'ke_spin_before : ',ke_spin_before - write(*,*) 'ke_spin_after : ',ke_spin_after - write(*,*) 'pe_before : ',pe_before - write(*,*) 'pe_after : ',pe_after - write(*,*) 'Qloss : ',Qloss - write(*,*) '***************************************************' + !write(*,*) '***************************************************' + !write(*,*) 'Original dis : ',norm2(x(:,2) - x(:,1)) + !write(*,*) 'r_max : ',r_max + !write(*,*) 'f_spin : ',f_spin + !write(*,*) '***************************************************' + !write(*,*) 'Energy balance so far: ' + !write(*,*) 'ke_frag_budget : ',ke_frag_budget + !write(*,*) 'ke_orbit_before: ',ke_orbit_before + !write(*,*) 'ke_orbit_after : ',ke_orbit_after + !write(*,*) 'ke_spin_before : ',ke_spin_before + !write(*,*) 'ke_spin_after : ',ke_spin_after + !write(*,*) 'pe_before : ',pe_before + !write(*,*) 'pe_after : ',pe_after + !write(*,*) 'Qloss : ',Qloss + !write(*,*) '***************************************************' if (ke_frag_budget < 0.0_DP) then write(*,*) 'Negative ke_frag_budget: ',ke_frag_budget r_max_start = r_max_start / 2 @@ -608,11 +610,11 @@ subroutine set_fragment_tan_vel(lerr) ! If we are over the energy budget, flag this as a failure so we can try again lerr = (ke_radial < 0.0_DP) - write(*,*) 'Tangential' - write(*,*) 'ke_frag_budget: ',ke_frag_budget - write(*,*) 'ke_frag_orbit : ',ke_frag_orbit - write(*,*) 'ke_frag_spin : ',ke_frag_spin - write(*,*) 'ke_radial : ',ke_radial + !write(*,*) 'Tangential' + !write(*,*) 'ke_frag_budget: ',ke_frag_budget + !write(*,*) 'ke_frag_orbit : ',ke_frag_orbit + !write(*,*) 'ke_frag_spin : ',ke_frag_spin + !write(*,*) 'ke_radial : ',ke_radial return diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index 05c3e2436..5cf8d20c7 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -75,7 +75,7 @@ module subroutine util_get_energy_momentum_system(self, param) ! Do the central body potential energy component first !$omp simd do i = 1, npl - associate(px => pl%xh(1,i), py => pl%xh(2,i), pz => pl%xh(3,i)) + associate(px => pl%xb(1,i), py => pl%xb(2,i), pz => pl%xb(3,i)) pecb(i) = -cb%Gmass * pl%mass(i) / sqrt(px**2 + py**2 + pz**2) end associate end do @@ -83,7 +83,7 @@ module subroutine util_get_energy_momentum_system(self, param) ! Do the potential energy between pairs of massive bodies do k = 1, pl%nplpl associate(ik => pl%k_plpl(1, k), jk => pl%k_plpl(2, k)) - pepl(k) = -pl%Gmass(ik) * pl%mass(jk) / norm2(pl%xh(:, jk) - pl%xh(:, ik)) + pepl(k) = -pl%Gmass(ik) * pl%mass(jk) / norm2(pl%xb(:, jk) - pl%xb(:, ik)) lstatpl(k) = (lstatus(ik) .and. lstatus(jk)) end associate end do @@ -91,7 +91,7 @@ module subroutine util_get_energy_momentum_system(self, param) system%ke_orbit = 0.5_DP * sum(kepl(1:npl), lstatus(:)) if (param%lrotation) system%ke_spin = 0.5_DP * sum(kespinpl(1:npl), lstatus(:)) - system%pe = sum(pepl(:), lstatpl(:)) + sum(pecb(2:npl), lstatus(2:npl)) + system%pe = sum(pepl(:), lstatpl(:)) + sum(pecb(1:npl), lstatus(1:npl)) ! Potential energy from the oblateness term if (param%loblatecb) then From 14d37ee8fca33b563ceae5559e104f98b171be7a Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 8 Aug 2021 17:15:50 -0400 Subject: [PATCH 8/8] Added missing central body terms to energy and momentum calculations --- src/util/util_get_energy_momentum.f90 | 40 +++++++++++++++++++++------ 1 file changed, 31 insertions(+), 9 deletions(-) diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index 5cf8d20c7..90f0d2242 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -16,6 +16,7 @@ module subroutine util_get_energy_momentum_system(self, param) integer(I4B) :: i, j integer(I8B) :: k real(DP) :: rmag, v2, rot2, oblpot, hx, hy, hz, hsx, hsy, hsz + real(DP) :: kecb, kespincb, Lcborbitx, Lcborbity, Lcborbitz, Lcbspinx, Lcbspiny, Lcbspinz 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 @@ -38,6 +39,13 @@ module subroutine util_get_energy_momentum_system(self, param) Lplspinz(:) = 0.0_DP lstatus(1:npl) = pl%status(1:npl) /= INACTIVE call pl%h2b(cb) + kecb = cb%mass * dot_product(cb%vb(:), cb%vb(:)) + hx = cb%xb(2) * cb%vb(3) - cb%xb(3) * cb%vb(2) + hy = cb%xb(3) * cb%vb(1) - cb%xb(1) * cb%vb(3) + hz = cb%xb(1) * cb%vb(2) - cb%xb(2) * cb%vb(1) + Lcborbitx = cb%mass * hx + Lcborbity = cb%mass * hy + Lcborbitz = cb%mass * hz !!$omp simd private(v2, rot2, hx, hy, hz) do i = 1, npl v2 = dot_product(pl%vb(:,i), pl%vb(:,i)) @@ -55,6 +63,17 @@ module subroutine util_get_energy_momentum_system(self, param) end do if (param%lrotation) then + kespincb = cb%mass * cb%Ip(3) * cb%radius**2 * dot_product(cb%rot(:), cb%rot(:)) + ! For simplicity, we always assume that the rotation pole is the 3rd principal axis + hsx = cb%Ip(3) * cb%radius**2 * cb%rot(1) + hsy = cb%Ip(3) * cb%radius**2 * cb%rot(2) + hsz = cb%Ip(3) * cb%radius**2 * cb%rot(3) + + ! Angular momentum from spin + Lcbspinx = cb%mass * hsx + Lcbspiny = cb%mass * hsy + Lcbspinz = cb%mass * hsz + do i = 1, npl rot2 = dot_product(pl%rot(:,i), pl%rot(:,i)) ! For simplicity, we always assume that the rotation pole is the 3rd principal axis @@ -69,6 +88,7 @@ module subroutine util_get_energy_momentum_system(self, param) kespinpl(i) = pl%mass(i) * pl%Ip(3, i) * pl%radius(i)**2 * rot2 end do else + kespincb = 0.0_DP kespinpl(:) = 0.0_DP end if @@ -88,8 +108,8 @@ module subroutine util_get_energy_momentum_system(self, param) end associate end do - system%ke_orbit = 0.5_DP * sum(kepl(1:npl), lstatus(:)) - if (param%lrotation) system%ke_spin = 0.5_DP * sum(kespinpl(1:npl), lstatus(:)) + 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%pe = sum(pepl(:), lstatpl(:)) + sum(pecb(1:npl), lstatus(1:npl)) @@ -103,13 +123,15 @@ module subroutine util_get_energy_momentum_system(self, param) system%pe = system%pe + oblpot end if - system%Lorbit(1) = sum(Lplorbitx(1:npl), lstatus(1:npl)) - system%Lorbit(2) = sum(Lplorbity(1:npl), lstatus(1:npl)) - system%Lorbit(3) = sum(Lplorbitz(1:npl), lstatus(1:npl)) - - system%Lspin(1) = sum(Lplspinx(1:npl), lstatus(1:npl)) - system%Lspin(2) = sum(Lplspiny(1:npl), lstatus(1:npl)) - system%Lspin(3) = sum(Lplspinz(1:npl), lstatus(1:npl)) + system%Lorbit(1) = Lcborbitx + sum(Lplorbitx(1:npl), lstatus(1:npl)) + system%Lorbit(2) = Lcborbity + sum(Lplorbity(1:npl), lstatus(1:npl)) + system%Lorbit(3) = Lcborbitz + sum(Lplorbitz(1:npl), lstatus(1:npl)) + + if (param%lrotation) then + system%Lspin(1) = Lcbspinx + sum(Lplspinx(1:npl), lstatus(1:npl)) + system%Lspin(2) = Lcbspiny + sum(Lplspiny(1:npl), lstatus(1:npl)) + system%Lspin(3) = Lcbspinz + sum(Lplspinz(1:npl), lstatus(1:npl)) + end if end associate return