diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index 3c7a7aa4a..0dd7ddd64 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -34,7 +34,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, real(DP), dimension(NDIM) :: Ltot_after real(DP) :: Etot_before, ke_orbit_before, ke_spin_before, pe_before, Lmag_before real(DP) :: Etot_after, ke_orbit_after, ke_spin_after, pe_after, Lmag_after, dEtot, dLmag - real(DP), dimension(NDIM) :: L_frag_tot, L_frag_orb, L_frag_spin, L_frag_budget + real(DP), dimension(NDIM) :: L_frag_tot, L_frag_orb, L_frag_spin, L_frag_budget, Lorbit_before, Lorbit_after, Lspin_before, Lspin_after real(DP) :: ke_frag_budget, ke_frag_orbit, ke_radial, ke_frag_spin, ke_avg_deficit, ke_avg_deficit_old real(DP), dimension(NDIM) :: x_col_unit, y_col_unit, z_col_unit character(len=*), parameter :: fmtlabel = "(A14,10(ES11.4,1X,:))" @@ -125,12 +125,12 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, v_r_mag(:) = 0.0_DP call set_fragment_position_vectors() + do concurrent (ii = 1:nfrag) + vb_frag(:, ii) = vcom(:) + end do + call calculate_system_energy(linclude_fragments=.true.) ke_frag_budget = -dEtot - Qloss - L_frag_budget(:) = Ltot_after(:) - Ltot_before(:) - do ii = 1, nfrag - L_frag_budget(:) = L_frag_budget(:) + m_frag(ii) * (xb_frag(:, ii) .cross. vcom(:)) - end do call define_coordinate_system() call set_fragment_tan_vel(lfailure) @@ -168,11 +168,11 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, try = try + 1 end do call restore_scale_factors() - call calculate_system_energy(linclude_fragments=.true.) write(*, "(' -------------------------------------------------------------------------------------')") write(*, "(' Final diagnostic')") write(*, "(' -------------------------------------------------------------------------------------')") + call calculate_system_energy(linclude_fragments=.true.) if (lfailure) then write(*,*) "symba_frag_pos failed after: ",try," tries" do ii = 1, nfrag @@ -210,12 +210,12 @@ subroutine set_scale_factors() vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot ! 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 = mtot - vscale = sqrt(0.5_DP) * (norm2(v(:,1)) + norm2(v(:,2))) + mscale = mtot + vscale = sqrt(Escale / mscale) tscale = dscale / vscale Lscale = mscale * dscale * vscale - Escale = mscale * vscale**2 xcom(:) = xcom(:) / dscale vcom(:) = vcom(:) / vscale @@ -471,6 +471,8 @@ subroutine calculate_system_energy(linclude_fragments) ! Calculate the current fragment energy and momentum balances if (linclude_fragments) then + Lorbit_after(:) = tmpsys%Lorbit + Lspin_after(:) = tmpsys%Lspin Ltot_after(:) = tmpsys%Lorbit(:) + tmpsys%Lspin(:) Lmag_after = norm2(Ltot_after(:)) ke_orbit_after = tmpsys%ke_orbit @@ -480,6 +482,8 @@ subroutine calculate_system_energy(linclude_fragments) dEtot = Etot_after - Etot_before dLmag = norm2(Ltot_after(:) - Ltot_before(:)) else + Lorbit_before(:) = tmpsys%Lorbit + Lspin_before(:) = tmpsys%Lspin Ltot_before(:) = tmpsys%Lorbit(:) + tmpsys%Lspin(:) Lmag_before = norm2(Ltot_before(:)) ke_orbit_before = tmpsys%ke_orbit @@ -504,7 +508,7 @@ subroutine calculate_fragment_ang_mtm() L_frag_spin(:) = 0.0_DP do i = 1, nfrag - L_frag_orb(:) = L_frag_orb(:) + m_frag(i) * (x_frag(:, i) .cross. v_frag(:, i)) + L_frag_orb(:) = L_frag_orb(:) + m_frag(i) * (x_frag(:, i) .cross. v_frag(:, i)) L_frag_spin(:) = L_frag_spin(:) + m_frag(i) * rad_frag(i)**2 * Ip_frag(:, i) * rot_frag(:, i) end do @@ -631,6 +635,22 @@ subroutine set_fragment_tan_vel(lerr) lerr = .false. + ! 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 @@ -715,7 +735,7 @@ subroutine set_fragment_tan_vel(lerr) lerr = (ke_radial < 0.0_DP) write(*,*) 'Tangential' write(*,*) 'Failure? ',lerr - write(*,*) '|L_remainder| : ',.mag.L_remainder(:) + write(*,*) '|L_remainder| : ',.mag.L_remainder(:) / Lmag_before write(*,*) 'ke_frag_budget: ',ke_frag_budget write(*,*) 'ke_frag_spin : ',ke_frag_spin write(*,*) 'ke_tangential : ',ke_frag_orbit @@ -784,17 +804,16 @@ function solve_fragment_tan_vel(lerr, v_t_mag_input) result(v_t_mag_output) do i = 1, nfrag if (i <= 2 * NDIM) then ! The tangential velocities of the first set of bodies will be the unknowns we will solve for to satisfy the constraints A(1:3, i) = m_frag(i) * v_t_unit(:, i) - L(:) = v_r_unit(:, i) .cross. v_t_unit(:, i) - A(4:6, i) = m_frag(i) * rmag(i) * L(:) + A(4:6, i) = m_frag(i) * rmag(i) * (v_r_unit(:, i) .cross. v_t_unit(:, i)) else if (present(v_t_mag_input)) then vtmp(:) = v_t_mag_input(i - 6) * v_t_unit(:, i) L_lin_others(:) = L_lin_others(:) + m_frag(i) * vtmp(:) - L(:) = m_frag(i) * (x_frag(:, i) .cross. vtmp(:)) + L(:) = m_frag(i) * (x_frag(:, i) .cross. vtmp(:)) L_orb_others(:) = L_orb_others(:) + L(:) end if end do b(1:3) = -L_lin_others(:) - b(4:6) = L_frag_budget(:) - L_frag_spin(:) - L_orb_others(:) + b(4:6) = L_frag_orb(:) - L_orb_others(:) allocate(v_t_mag_output(nfrag)) 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(:) diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 0e6c69440..80a9550bb 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -454,11 +454,11 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec do k = 1, nenc if (lcollision(k)) self%status(k) = COLLISION self%t(k) = t - self%x1(:,k) = pl%xh(:,ind1(k)) - self%v1(:,k) = pl%vb(:,ind1(k)) - system%cb%vb(:) + self%x1(:,k) = pl%xh(:,ind1(k)) + system%cb%xb(:) + self%v1(:,k) = pl%vb(:,ind1(k)) if (isplpl) then - self%x2(:,k) = pl%xh(:,ind2(k)) - self%v2(:,k) = pl%vb(:,ind2(k)) - system%cb%vb(:) + self%x2(:,k) = pl%xh(:,ind2(k)) + system%cb%xb(:) + self%v2(:,k) = pl%vb(:,ind2(k)) if (lcollision(k)) then ! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional family if (pl%lcollision(ind1(k)) .or. pl%lcollision(ind2(k))) call pl%make_family([ind1(k),ind2(k)]) @@ -469,8 +469,8 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec pl%status([ind1(k), ind2(k)]) = COLLISION end if else - self%x2(:,k) = tp%xh(:,ind2(k)) - self%v2(:,k) = tp%vb(:,ind2(k)) - system%cb%vb(:) + self%x2(:,k) = tp%xh(:,ind2(k)) + system%cb%xb(:) + self%v2(:,k) = tp%vb(:,ind2(k)) if (lcollision(k)) then tp%status(ind2(k)) = DISCARDED_PLR tp%ldiscard(ind2(k)) = .true. @@ -539,7 +539,7 @@ pure elemental function symba_collision_check_one(xr, yr, zr, vxr, vyr, vzr, Gmt end function symba_collision_check_one - function symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v, mass, radius, L_spin, Ip) result(lflag) + function symba_collision_consolidate_familes(pl, cb, param, idx_parent, family, x, v, mass, radius, L_spin, Ip) result(lflag) !! author: David A. Minton !! !! Loops through the pl-pl collision list and groups families together by index. Outputs the indices of all family members, @@ -547,6 +547,7 @@ function symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v implicit none ! Arguments class(symba_pl), intent(inout) :: pl !! SyMBA massive body object + class(symba_cb), intent(inout) :: cb !! SyMBA central body object class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(2), intent(inout) :: idx_parent !! Index of the two bodies considered the "parents" of the collision integer(I4B), dimension(:), allocatable, intent(out) :: family !! List of indices of all bodies inovlved in the collision @@ -611,7 +612,7 @@ function symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v ! Find the barycenter of each body along with its children, if it has any do j = 1, 2 - x(:, j) = pl%xh(:, idx_parent(j)) + x(:, j) = pl%xh(:, idx_parent(j)) + cb%xb(:) v(:, j) = pl%vb(:, idx_parent(j)) ! Assume principal axis rotation about axis corresponding to highest moment of inertia (3rd Ip) if (param%lrotation) then @@ -624,7 +625,7 @@ function symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v idx_child = parent_child_index_array(j)%idx(i + 1) if (.not. pl%lcollision(idx_child)) cycle mchild = pl%mass(idx_child) - xchild(:) = pl%xh(:, idx_child) + xchild(:) = pl%xh(:, idx_child) + cb%xb(:) vchild(:) = pl%vb(:, idx_child) volchild = (4.0_DP / 3.0_DP) * PI * pl%radius(idx_child)**3 volume(j) = volume(j) + volchild @@ -944,68 +945,71 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) integer(I4B), parameter :: NRES = 3 !! Number of collisional product results real(DP), dimension(NRES) :: mass_res - associate(plpl_collisions => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, cb => system%cb) + associate(plpl_collisions => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2) select type(pl => system%pl) class is (symba_pl) - do i = 1, ncollisions - idx_parent(1) = pl%kin(idx1(i))%parent - idx_parent(2) = pl%kin(idx2(i))%parent - lgoodcollision = symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v, mass, radius, L_spin, Ip) - if (.not. lgoodcollision) cycle - if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved - - ! Convert from DH to barycentric - x(:,1) = x(:,1) + cb%xb(:) - x(:,2) = x(:,2) + cb%xb(:) - - ! Convert all quantities to SI units and determine which of the pair is the projectile vs. target before sending them - ! to symba_regime - if (mass(1) > mass(2)) then - jtarg = 1 - jproj = 2 - else - jtarg = 2 - jproj = 1 - end if - mass_si(:) = (mass(:)) * param%MU2KG !! The collective mass of the parent and its children - radius_si(:) = radius(:) * param%DU2M !! The collective radius of the parent and its children - x1_si(:) = plpl_collisions%x1(:,i) * param%DU2M !! The position of the parent from inside the step (at collision) - v1_si(:) = plpl_collisions%v1(:,i) * param%DU2M / param%TU2S !! The velocity of the parent from inside the step (at collision) - x2_si(:) = plpl_collisions%x2(:,i) * param%DU2M !! The position of the parent from inside the step (at collision) - v2_si(:) = plpl_collisions%v2(:,i) * param%DU2M / param%TU2S !! The velocity of the parent from inside the step (at collision) - density_si(:) = mass_si(:) / (4.0_DP / 3._DP * PI * radius_si(:)**3) !! The collective density of the parent and its children - Mcb_si = cb%mass * param%MU2KG - mtiny_si = (param%GMTINY / param%GU) * param%MU2KG + select type (cb => system%cb) + class is (symba_cb) + do i = 1, ncollisions + idx_parent(1) = pl%kin(idx1(i))%parent + idx_parent(2) = pl%kin(idx2(i))%parent + lgoodcollision = symba_collision_consolidate_familes(pl, cb, param, idx_parent, family, x, v, mass, radius, L_spin, Ip) + if (.not. lgoodcollision) cycle + if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved + + ! Convert from DH to barycentric + x(:,1) = x(:,1) + cb%xb(:) + x(:,2) = x(:,2) + cb%xb(:) + + ! Convert all quantities to SI units and determine which of the pair is the projectile vs. target before sending them + ! to symba_regime + if (mass(1) > mass(2)) then + jtarg = 1 + jproj = 2 + else + jtarg = 2 + jproj = 1 + end if + mass_si(:) = (mass(:)) * param%MU2KG !! The collective mass of the parent and its children + radius_si(:) = radius(:) * param%DU2M !! The collective radius of the parent and its children + x1_si(:) = plpl_collisions%x1(:,i) * param%DU2M !! The position of the parent from inside the step (at collision) + v1_si(:) = plpl_collisions%v1(:,i) * param%DU2M / param%TU2S !! The velocity of the parent from inside the step (at collision) + x2_si(:) = plpl_collisions%x2(:,i) * param%DU2M !! The position of the parent from inside the step (at collision) + v2_si(:) = plpl_collisions%v2(:,i) * param%DU2M / param%TU2S !! The velocity of the parent from inside the step (at collision) + density_si(:) = mass_si(:) / (4.0_DP / 3._DP * PI * radius_si(:)**3) !! The collective density of the parent and its children + Mcb_si = cb%mass * param%MU2KG + mtiny_si = (param%GMTINY / param%GU) * param%MU2KG + + mass_res(:) = 0.0_DP - mass_res(:) = 0.0_DP - - mtot = sum(mass_si(:)) - dentot = sum(mass_si(:) * density_si(:)) / mtot - - !! Use the positions and velocities of the parents from indside the step (at collision) to calculate the collisional regime - call fragmentation_regime(Mcb_si, mass_si(jtarg), mass_si(jproj), radius_si(jtarg), radius_si(jproj), x1_si(:), x2_si(:),& - v1_si(:), v2_si(:), density_si(jtarg), density_si(jproj), regime, mlr, mslr, mtiny_si, Qloss) - - mass_res(1) = min(max(mlr, 0.0_DP), mtot) - mass_res(2) = min(max(mslr, 0.0_DP), mtot) - mass_res(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot) - mass_res(:) = (mass_res(:) / param%MU2KG) - Qloss = Qloss * (param%TU2S / param%DU2M)**2 / param%MU2KG - - select case (regime) - case (COLLRESOLVE_REGIME_DISRUPTION) - status = symba_collision_casedisruption(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) - case (COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - status = symba_collision_casesupercatastrophic(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) - case (COLLRESOLVE_REGIME_HIT_AND_RUN) - status = symba_collision_casehitandrun(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) - case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - status = symba_collision_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) - case default - write(*,*) "Error in symba_collision, unrecognized collision regime" - call util_exit(FAILURE) - end select - end do + mtot = sum(mass_si(:)) + dentot = sum(mass_si(:) * density_si(:)) / mtot + + !! Use the positions and velocities of the parents from indside the step (at collision) to calculate the collisional regime + call fragmentation_regime(Mcb_si, mass_si(jtarg), mass_si(jproj), radius_si(jtarg), radius_si(jproj), x1_si(:), x2_si(:),& + v1_si(:), v2_si(:), density_si(jtarg), density_si(jproj), regime, mlr, mslr, mtiny_si, Qloss) + + mass_res(1) = min(max(mlr, 0.0_DP), mtot) + mass_res(2) = min(max(mslr, 0.0_DP), mtot) + mass_res(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot) + mass_res(:) = (mass_res(:) / param%MU2KG) + Qloss = Qloss * (param%TU2S / param%DU2M)**2 / param%MU2KG + + select case (regime) + case (COLLRESOLVE_REGIME_DISRUPTION) + status = symba_collision_casedisruption(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) + case (COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + status = symba_collision_casesupercatastrophic(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) + case (COLLRESOLVE_REGIME_HIT_AND_RUN) + status = symba_collision_casehitandrun(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) + case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) + status = symba_collision_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) + case default + write(*,*) "Error in symba_collision, unrecognized collision regime" + call util_exit(FAILURE) + end select + end do + end select end select end associate @@ -1031,22 +1035,25 @@ module subroutine symba_collision_resolve_mergers(self, system, param) logical :: lgoodcollision integer(I4B) :: i, status - associate(plpl_collisions => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, cb => system%cb) + associate(plpl_collisions => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2) select type(pl => system%pl) class is (symba_pl) - do i = 1, ncollisions - idx_parent(1) = pl%kin(idx1(i))%parent - idx_parent(2) = pl%kin(idx2(i))%parent - lgoodcollision = symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v, mass, radius, L_spin, Ip) - if (.not. lgoodcollision) cycle - if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved - - ! Convert from DH to barycentric - x(:,1) = x(:,1) + cb%xb(:) - x(:,2) = x(:,2) + cb%xb(:) - - status = symba_collision_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) - end do + select type(cb => system%cb) + class is (symba_cb) + do i = 1, ncollisions + idx_parent(1) = pl%kin(idx1(i))%parent + idx_parent(2) = pl%kin(idx2(i))%parent + lgoodcollision = symba_collision_consolidate_familes(pl, cb, param, idx_parent, family, x, v, mass, radius, L_spin, Ip) + if (.not. lgoodcollision) cycle + if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved + + ! Convert from DH to barycentric + x(:,1) = x(:,1) + cb%xb(:) + x(:,2) = x(:,2) + cb%xb(:) + + status = symba_collision_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) + end do + end select end select end associate diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 70c0898a5..f6f7c0e6b 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -97,10 +97,10 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body) Ltot(:) = 0.0_DP do i = 1, pl%nbody - Lpl(:) = pL%mass(i) * pl%xb(:,i) .cross. pl%vb(:, i) + Lpl(:) = pL%mass(i) * (pl%xb(:,i) .cross. pl%vb(:, i)) Ltot(:) = Ltot(:) + Lpl(:) end do - Ltot(:) = Ltot(:) + cb%mass * cb%xb(:) .cross. cb%vb(:) + Ltot(:) = Ltot(:) + cb%mass * (cb%xb(:) .cross. cb%vb(:)) call pl%b2h(cb) oldstat = pl%status(ipl) pl%status(ipl) = INACTIVE @@ -108,20 +108,20 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body) pl%status(ipl) = oldstat do i = 1, pl%nbody if (i == ipl) cycle - Lpl(:) = pl%mass(i) * pl%xb(:,i) .cross. pl%vb(:, i) + Lpl(:) = pl%mass(i) * (pl%xb(:,i) .cross. pl%vb(:, i)) Ltot(:) = Ltot(:) - Lpl(:) end do - Ltot(:) = Ltot(:) - cb%mass * cb%xb(:) .cross. cb%vb(:) + Ltot(:) = Ltot(:) - cb%mass * (cb%xb(:) .cross. cb%vb(:)) system%Lescape(:) = system%Lescape(:) + Ltot(:) if (param%lrotation) system%Lescape(:) = system%Lescape + pl%mass(ipl) * pl%radius(ipl)**2 * pl%Ip(3, ipl) * pl%rot(:, ipl) else xcom(:) = (pl%mass(ipl) * pl%xb(:, ipl) + cb%mass * cb%xb(:)) / (cb%mass + pl%mass(ipl)) vcom(:) = (pl%mass(ipl) * pl%vb(:, ipl) + cb%mass * cb%vb(:)) / (cb%mass + pl%mass(ipl)) - Lpl(:) = (pl%xb(:,ipl) - xcom(:)) .cross. pL%vb(:,ipl) - vcom(:) + Lpl(:) = (pl%xb(:,ipl) - xcom(:)) .cross. (pL%vb(:,ipl) - vcom(:)) if (param%lrotation) Lpl(:) = pl%mass(ipl) * (Lpl(:) + pl%radius(ipl)**2 * pl%Ip(3,ipl) * pl%rot(:, ipl)) - Lcb(:) = cb%mass * (cb%xb(:) - xcom(:)) .cross. (cb%vb(:) - vcom(:)) + Lcb(:) = cb%mass * ((cb%xb(:) - xcom(:)) .cross. (cb%vb(:) - vcom(:))) ke_orbit = ke_orbit + 0.5_DP * cb%mass * dot_product(cb%vb(:), cb%vb(:)) if (param%lrotation) ke_spin = ke_spin + 0.5_DP * cb%mass * cb%radius**2 * cb%Ip(3) * dot_product(cb%rot(:), cb%rot(:)) diff --git a/src/tides/tides_getacch_pl.f90 b/src/tides/tides_getacch_pl.f90 index f0bf64cc7..4feb76221 100644 --- a/src/tides/tides_getacch_pl.f90 +++ b/src/tides/tides_getacch_pl.f90 @@ -48,8 +48,8 @@ module subroutine tides_kick_getacch_pl(self, system) Ftr = -3 / rmag**7 * (r5cbterm + r5plterm) - 3 * vmag / rmag * (Ptocb + Ptopl) F_T(:) = (Ftr + (Ptocb + Ptopl) * dot_product(v_unit, r_unit) / rmag) * r_unit(:) & - + Ptopl * (pl%rot(:,i) - theta_dot(:)) .cross. r_unit(:) & - + Ptocb * (cb%rot(:) - theta_dot(:)) .cross. r_unit(:) + + Ptopl * ((pl%rot(:,i) - theta_dot(:)) .cross. r_unit(:)) & + + Ptocb * ((cb%rot(:) - theta_dot(:)) .cross. r_unit(:)) cb%atide(:) = cb%atide(:) + F_T(:) / cb%Gmass pl%atide(:,i) = F_T(:) / pl%Gmass(i) end do diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index fa7cda43d..7535ea67b 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -40,20 +40,20 @@ module subroutine util_get_energy_momentum_system(self, param) lstatus(1:npl) = pl%status(1:npl) /= INACTIVE kecb = cb%mass * dot_product(cb%vb(:), cb%vb(:)) - Lcborbit(:) = cb%mass * cb%xb(:) .cross. cb%vb(:) + Lcborbit(:) = cb%mass * (cb%xb(:) .cross. cb%vb(:)) do concurrent (i = 1:npl, lstatus(i)) block ! We use a block construct to prevent generating temporary arrays for local variables - real(DP) :: v2, hx, hy, hz + real(DP) :: v2 + real(DP), dimension(NDIM) :: L + v2 = dot_product(pl%vb(:,i), pl%vb(:,i)) - hx = pl%xb(2,i) * pl%vb(3,i) - pl%xb(3,i) * pl%vb(2,i) - hy = pl%xb(3,i) * pl%vb(1,i) - pl%xb(1,i) * pl%vb(3,i) - hz = pl%xb(1,i) * pl%vb(2,i) - pl%xb(2,i) * pl%vb(1,i) + L(:) = pl%mass(i) * (pl%xb(:,i) .cross. pl%vb(:,i)) ! Angular momentum from orbit - Lplorbitx(i) = pl%mass(i) * hx - Lplorbity(i) = pl%mass(i) * hy - Lplorbitz(i) = pl%mass(i) * hz + Lplorbitx(i) = L(1) + Lplorbity(i) = L(2) + Lplorbitz(i) = L(3) ! Kinetic energy from orbit and spin kepl(i) = pl%mass(i) * v2 @@ -67,21 +67,12 @@ module subroutine util_get_energy_momentum_system(self, param) Lcbspin(:) = cb%Ip(3) * cb%mass * cb%radius**2 * cb%rot(:) do concurrent (i = 1:npl, lstatus(i)) - block - real(DP) :: rot2, hsx, hsy, hsz - - rot2 = dot_product(pl%rot(:,i), pl%rot(:,i)) - ! For simplicity, we always assume that the rotation pole is the 3rd principal axis - hsx = pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(1,i) - hsy = pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(2,i) - hsz = pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(3,i) - - ! Angular momentum from spin - Lplspinx(i) = pl%mass(i) * hsx - Lplspiny(i) = pl%mass(i) * hsy - Lplspinz(i) = pl%mass(i) * hsz - kespinpl(i) = pl%mass(i) * pl%Ip(3, i) * pl%radius(i)**2 * rot2 - end block + ! Currently we assume that the rotation pole is the 3rd principal axis + ! Angular momentum from spin + Lplspinx(i) = pl%mass(i) * pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(1,i) + Lplspiny(i) = pl%mass(i) * pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(2,i) + Lplspinz(i) = pl%mass(i) * pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(3,i) + kespinpl(i) = pl%mass(i) * pl%Ip(3, i) * pl%radius(i)**2 * dot_product(pl%rot(:,i), pl%rot(:,i)) end do else kespincb = 0.0_DP