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