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

Commit

Permalink
Made some changes that allow the code to compile in gfortran
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Jan 13, 2023
1 parent 7e44513 commit e760cb4
Show file tree
Hide file tree
Showing 7 changed files with 130 additions and 126 deletions.
2 changes: 2 additions & 0 deletions CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,8 @@ INCLUDE(${CMAKE_MODULE_PATH}/SetMKL.cmake)
# RELEASE, PROFILING, and TESTING.
INCLUDE(${CMAKE_MODULE_PATH}/SetFortranFlags.cmake)

INCLUDE_DIRECTORIES($ENV{NETCDF_FORTRAN_HOME}/include;$ENV{NETCDF_HOME}/include)


# There is an error in CMAKE with this flag for pgf90. Unset it
GET_FILENAME_COMPONENT(FCNAME ${CMAKE_Fortran_COMPILER} NAME)
Expand Down
2 changes: 1 addition & 1 deletion src/collision/collision_resolve.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -331,7 +331,7 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status)
class is (swiftest_nbody_system)
select type(param)
class is (swiftest_parameters)
associate(pl => nbody_system%pl, pl_discards => nbody_system%pl_discards, info => nbody_system%pl%info, pl_adds => nbody_system%pl_adds, cb => nbody_system%cb, npl => pl%nbody, &
associate(pl => nbody_system%pl, pl_discards => nbody_system%pl_discards, info => nbody_system%pl%info, pl_adds => nbody_system%pl_adds, cb => nbody_system%cb, npl => nbody_system%pl%nbody, &
collider => nbody_system%collider, impactors => nbody_system%collider%impactors,fragments => nbody_system%collider%fragments)

! Add the impactors%id bodies to the subtraction list
Expand Down
50 changes: 26 additions & 24 deletions src/collision/collision_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -818,37 +818,39 @@ module subroutine collision_util_snapshot(self, param, nbody_system, t, arg)
select type(before_snap => snapshot%collider%before )
class is (swiftest_nbody_system)
select type(before_orig => nbody_system%collider%before)
class is (swiftest_nbody_system)
select type(plsub => before_orig%pl)
class is (swiftest_pl)
! Log the properties of the old and new bodies
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Removing bodies:")
do i = 1, plsub%nbody
write(message,*) trim(adjustl(plsub%info(i)%name)), " (", trim(adjustl(plsub%info(i)%particle_type)),")"
call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
end do

call move_alloc(before_orig%pl, before_snap%pl)
end select
class is (swiftest_nbody_system)
select type(plsub => before_orig%pl)
class is (swiftest_pl)
! Log the properties of the old and new bodies
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Removing bodies:")
do i = 1, plsub%nbody
write(message,*) trim(adjustl(plsub%info(i)%name)), " (", trim(adjustl(plsub%info(i)%particle_type)),")"
call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
end do

allocate(before_snap%pl, source=plsub)
end select
deallocate(before_orig%pl)
end select
end select


select type(after_snap => snapshot%collider%after )
class is (swiftest_nbody_system)
select type(after_orig => nbody_system%collider%after)
class is (swiftest_nbody_system)
select type(plnew => after_orig%pl)
class is (swiftest_pl)
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Adding bodies:")
do i = 1, plnew%nbody
write(message,*) trim(adjustl(plnew%info(i)%name)), " (", trim(adjustl(plnew%info(i)%particle_type)),")"
call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
end do
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "***********************************************************" // &
"***********************************************************")

call move_alloc(after_orig%pl, after_snap%pl)
end select
select type(plnew => after_orig%pl)
class is (swiftest_pl)
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Adding bodies:")
do i = 1, plnew%nbody
write(message,*) trim(adjustl(plnew%info(i)%name)), " (", trim(adjustl(plnew%info(i)%particle_type)),")"
call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
end do
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "***********************************************************" // &
"***********************************************************")
allocate(after_snap%pl, source=plnew)
end select
deallocate(after_orig%pl)
end select
end select

Expand Down
8 changes: 4 additions & 4 deletions src/encounter/encounter_check.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -490,7 +490,7 @@ pure subroutine encounter_check_all_sweep_one(i, n, xi, yi, zi, vxi, vyi, vzi, x
end subroutine encounter_check_all_sweep_one


pure subroutine encounter_check_all_triangular_one(i, n, xi, yi, zi, vxi, vyi, vzi, x, y, z, vx, vy, vz, renci, renc, &
subroutine encounter_check_all_triangular_one(i, n, xi, yi, zi, vxi, vyi, vzi, x, y, z, vx, vy, vz, renci, renc, &
dt, ind_arr, lenci)
!! author: David A. Minton
!!
Expand Down Expand Up @@ -1023,8 +1023,8 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt
integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for the other body in each encounter candidate pair
logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical array indicating which pairs are approaching
! Internals
integer(I4B) :: i, nbox, dim
integer(I8B) :: k, itmp
integer(I4B) :: i, dim, itmp, nbox
integer(I8B) :: k
logical, dimension(n) :: loverlap
logical, dimension(2*n) :: lencounteri
real(DP), dimension(2*n) :: xind, yind, zind, vxind, vyind, vzind, rencind
Expand Down Expand Up @@ -1062,7 +1062,7 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt
if (loverlap(i)) then
ibeg = self%aabb(dim)%ibeg(i) + 1_I8B
iend = self%aabb(dim)%iend(i) - 1_I8B
nbox = iend - ibeg + 1
nbox = int(iend - ibeg + 1, kind=I4B)
lencounteri(ibeg:iend) = .true.
call encounter_check_all_sweep_one(i, nbox, x(1,i), x(2,i), x(3,i), v(1,i), v(2,i), v(3,i), &
xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),&
Expand Down
4 changes: 2 additions & 2 deletions src/fraggle/fraggle_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -312,9 +312,9 @@ module subroutine fraggle_generate_pos_vec(collider)
istart = 2
end if

call random_number(fragments%rc(:,istart:nfrag))
do concurrent(i = istart:nfrag, loverlap(i))
! Make a random cloud
call random_number(fragments%rc(:,i))

! Make the fragment cloud symmertic about 0
fragments%rc(:,i) = 2 *(fragments%rc(:,i) - 0.5_DP)
Expand Down Expand Up @@ -407,7 +407,7 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param)
fragments%rotmag(1) = .mag.fragments%rot(:,1)
! Add in some random spin noise. The magnitude will be scaled by the before-after amount and the direction will be random
call random_number(fragments%rotmag(2:nfrag))
do concurrent(i = 2:nfrag)
do i = 2,nfrag
call random_number(rotdir)
rotdir = rotdir - 0.5_DP
rotdir = .unit. rotdir
Expand Down
186 changes: 93 additions & 93 deletions src/misc/solver_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,17 +23,17 @@ module solver
module procedure solve_linear_system_qp
end interface

interface
module function solve_rkf45(f, y0in, t1, dt0, tol) result(y1)
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 solve_rkf45
end interface
! interface
! module function solve_rkf45(f, y0in, t1, dt0, tol) result(y1)
! 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 solve_rkf45
! end interface


contains
Expand Down Expand Up @@ -171,93 +171,93 @@ function ge_wpp(A, b) result(u) ! gaussian eliminate with partial pivoting
end function ge_wpp


module function solve_rkf45(f, y0in, t1, dt0, tol) result(y1)
!! author: David A. Minton
!!
!! Implements the 4th order Runge-Kutta-Fehlberg ODE solver for initial value problems of the form f=dy/dt, y0 = y(t=0), solving for y1 = y(t=t1). Uses a 5th order adaptive step size control.
!! Uses a lambda function object as defined in the lambda_function module
implicit none
! Arguments
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
! Result
real(DP), dimension(:), allocatable :: y1 !! Final result
! Internals
integer(I4B), parameter :: MAXREDUX = 1000 !! Maximum number of times step size can be reduced
real(DP), parameter :: DTFAC = 0.95_DP !! Step size reduction safety factor (Value just under 1.0 to prevent adaptive step size control from discarding steps too aggressively)
integer(I4B), parameter :: RKS = 6 !! Number of RK stages
real(DP), dimension(RKS, RKS - 1), parameter :: rkf45_btab = reshape( & !! Butcher tableau for Runge-Kutta-Fehlberg method
(/ 1./4., 1./4., 0., 0., 0., 0.,&
3./8., 3./32., 9./32., 0., 0., 0.,&
12./13., 1932./2197., -7200./2197., 7296./2197., 0., 0.,&
1., 439./216., -8., 3680./513., -845./4104., 0.,&
1./2., -8./27., 2., -3544./2565., 1859./4104., -11./40./), shape(rkf45_btab))
real(DP), dimension(RKS), parameter :: rkf4_coeff = (/ 25./216., 0., 1408./2565. , 2197./4104. , -1./5., 0. /)
real(DP), dimension(RKS), parameter :: rkf5_coeff = (/ 16./135., 0., 6656./12825., 28561./56430., -9./50., 2./55. /)
real(DP), dimension(:, :), allocatable :: k !! Runge-Kutta coefficient vector
real(DP), dimension(:), allocatable :: ynorm !! Normalized y value used for adaptive step size control
real(DP), dimension(:), allocatable :: y0 !! Value of y at the beginning of each substep
integer(I4B) :: Nvar !! Number of variables in problem
integer(I4B) :: rkn !! Runge-Kutta loop index
real(DP) :: t, x1, dt, trem !! Current time, step size and total time remaining
real(DP) :: s, yerr, yscale !! Step size reduction factor, error in dependent variable, and error scale factor
integer(I4B) :: i
! module function solve_rkf45(f, y0in, t1, dt0, tol) result(y1)
! !! author: David A. Minton
! !!
! !! Implements the 4th order Runge-Kutta-Fehlberg ODE solver for initial value problems of the form f=dy/dt, y0 = y(t=0), solving for y1 = y(t=t1). Uses a 5th order adaptive step size control.
! !! Uses a lambda function object as defined in the lambda_function module
! implicit none
! ! Arguments
! 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
! ! Result
! real(DP), dimension(:), allocatable :: y1 !! Final result
! ! Internals
! integer(I4B), parameter :: MAXREDUX = 1000 !! Maximum number of times step size can be reduced
! real(DP), parameter :: DTFAC = 0.95_DP !! Step size reduction safety factor (Value just under 1.0 to prevent adaptive step size control from discarding steps too aggressively)
! integer(I4B), parameter :: RKS = 6 !! Number of RK stages
! real(DP), dimension(RKS, RKS - 1), parameter :: rkf45_btab = reshape( & !! Butcher tableau for Runge-Kutta-Fehlberg method
! (/ 1./4., 1./4., 0., 0., 0., 0.,&
! 3./8., 3./32., 9./32., 0., 0., 0.,&
! 12./13., 1932./2197., -7200./2197., 7296./2197., 0., 0.,&
! 1., 439./216., -8., 3680./513., -845./4104., 0.,&
! 1./2., -8./27., 2., -3544./2565., 1859./4104., -11./40./), shape(rkf45_btab))
! real(DP), dimension(RKS), parameter :: rkf4_coeff = (/ 25./216., 0., 1408./2565. , 2197./4104. , -1./5., 0. /)
! real(DP), dimension(RKS), parameter :: rkf5_coeff = (/ 16./135., 0., 6656./12825., 28561./56430., -9./50., 2./55. /)
! real(DP), dimension(:, :), allocatable :: k !! Runge-Kutta coefficient vector
! real(DP), dimension(:), allocatable :: ynorm !! Normalized y value used for adaptive step size control
! real(DP), dimension(:), allocatable :: y0 !! Value of y at the beginning of each substep
! integer(I4B) :: Nvar !! Number of variables in problem
! integer(I4B) :: rkn !! Runge-Kutta loop index
! real(DP) :: t, x1, dt, trem !! Current time, step size and total time remaining
! real(DP) :: s, yerr, yscale !! Step size reduction factor, error in dependent variable, and error scale factor
! integer(I4B) :: i

allocate(y0, source=y0in)
allocate(y1, mold=y0)
allocate(ynorm, mold=y0)
Nvar = size(y0)
allocate(k(Nvar, RKS))
! allocate(y0, source=y0in)
! allocate(y1, mold=y0)
! allocate(ynorm, mold=y0)
! Nvar = size(y0)
! allocate(k(Nvar, RKS))

dt = dt0
! dt = dt0

trem = t1
t = 0._DP
do
yscale = norm2(y0(:))
do i = 1, MAXREDUX
select type(f)
class is (lambda_obj_tvar)
do rkn = 1, RKS
y1(:) = y0(:) + matmul(k(:, 1:rkn - 1), rkf45_btab(2:rkn, rkn - 1))
if (rkn == 1) then
x1 = t
else
x1 = t + rkf45_btab(1,rkn-1)
end if
k(:, rkn) = dt * f%evalt(y1(:), t)
end do
class is (lambda_obj)
do rkn = 1, RKS
y1(:) = y0(:) + matmul(k(:, 1:rkn - 1), rkf45_btab(2:rkn, rkn - 1))
k(:, rkn) = dt * f%eval(y1(:))
end do
end select
! Now determine if the step size needs adjusting
ynorm(:) = matmul(k(:,:), (rkf5_coeff(:) - rkf4_coeff(:))) / yscale
yerr = norm2(ynorm(:))
s = (tol / (2 * yerr))**(0.25_DP)
dt = min(s * DTFAC * dt, trem) ! Alter step size either up or down, but never bigger than the remaining time
if (s >= 1.0_DP) exit ! Good step!
if (i == MAXREDUX) then
write(*,*) "Something has gone wrong in util_solve_rkf45!! Step size reduction has gone too far this time!"
call base_util_exit(FAILURE)
end if
end do
! trem = t1
! t = 0._DP
! do
! yscale = norm2(y0(:))
! do i = 1, MAXREDUX
! select type(f)
! class is (lambda_obj_tvar)
! do rkn = 1, RKS
! y1(:) = y0(:) + matmul(k(:, 1:rkn - 1), rkf45_btab(2:rkn, rkn - 1))
! if (rkn == 1) then
! x1 = t
! else
! x1 = t + rkf45_btab(1,rkn-1)
! end if
! k(:, rkn) = dt * f%evalt(y1(:), t)
! end do
! class is (lambda_obj)
! do rkn = 1, RKS
! y1(:) = y0(:) + matmul(k(:, 1:rkn - 1), rkf45_btab(2:rkn, rkn - 1))
! k(:, rkn) = dt * f%eval(y1(:))
! end do
! end select
! ! Now determine if the step size needs adjusting
! ynorm(:) = matmul(k(:,:), (rkf5_coeff(:) - rkf4_coeff(:))) / yscale
! yerr = norm2(ynorm(:))
! s = (tol / (2 * yerr))**(0.25_DP)
! dt = min(s * DTFAC * dt, trem) ! Alter step size either up or down, but never bigger than the remaining time
! if (s >= 1.0_DP) exit ! Good step!
! if (i == MAXREDUX) then
! write(*,*) "Something has gone wrong in util_solve_rkf45!! Step size reduction has gone too far this time!"
! call base_util_exit(FAILURE)
! end if
! end do

! Compute new value then step ahead in time
y1(:) = y0(:) + matmul(k(:, :), rkf4_coeff(:))
trem = trem - dt
t = t + dt
if (trem <= 0._DP) exit
y0(:) = y1(:)
end do
! ! Compute new value then step ahead in time
! y1(:) = y0(:) + matmul(k(:, :), rkf4_coeff(:))
! trem = trem - dt
! t = t + dt
! if (trem <= 0._DP) exit
! y0(:) = y1(:)
! end do

return
end function solve_rkf45
! return
! end function solve_rkf45


end module solver
Loading

0 comments on commit e760cb4

Please sign in to comment.