Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Merge branch 'debug'
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Aug 9, 2021
2 parents d2384a7 + 14d37ee commit 945b95d
Show file tree
Hide file tree
Showing 7 changed files with 159 additions and 66 deletions.
2 changes: 1 addition & 1 deletion Makefile.Defines
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
44 changes: 25 additions & 19 deletions src/fragmentation/fragmentation.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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."
Expand All @@ -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
Expand Down Expand Up @@ -166,27 +167,27 @@ 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 = mtot
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
rot(:,i) = L_spin(:,i) / (mass(i) * radius(i)**2 * Ip(3, i))
end do

m_frag = m_frag / mscale
rad_frag = rad_frag / rscale
rad_frag = rad_frag / dscale
Qloss = Qloss / Escale

return
Expand All @@ -201,23 +202,23 @@ 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
rot(:,i) = L_spin(:,i) * (mass(i) * radius(i)**2 * Ip(3, i))
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

Expand All @@ -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
Expand Down Expand Up @@ -334,22 +335,27 @@ 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)
if (allocated(tmpparam)) deallocate(tmpparam)
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)
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) = param%GU * m_frag(1:nfrag)
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)
Expand Down
73 changes: 41 additions & 32 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down
4 changes: 4 additions & 0 deletions src/symba/symba_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
50 changes: 36 additions & 14 deletions src/util/util_get_energy_momentum.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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))
Expand All @@ -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
Expand All @@ -69,47 +88,50 @@ 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

! 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))
pecb(i) = -cb%mass * pl%mass(i) / sqrt(px**2 + py**2 + pz**2)
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

! 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%xb(:, jk) - pl%xb(:, ik))
lstatpl(k) = (lstatus(ik) .and. lstatus(jk))
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(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
!$omp simd
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)
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))
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
Expand Down
Loading

0 comments on commit 945b95d

Please sign in to comment.