From e760cb489613ccf1ec34478470c651aa6ba07204 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 12 Jan 2023 19:44:33 -0500 Subject: [PATCH] Made some changes that allow the code to compile in gfortran --- CMakeLists.txt | 2 + src/collision/collision_resolve.f90 | 2 +- src/collision/collision_util.f90 | 50 ++++---- src/encounter/encounter_check.f90 | 8 +- src/fraggle/fraggle_generate.f90 | 4 +- src/misc/solver_module.f90 | 186 ++++++++++++++-------------- src/swiftest/swiftest_io.f90 | 4 +- 7 files changed, 130 insertions(+), 126 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 5d4c8637f..9ce5af9cf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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) diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index 8c4c6505a..ae928b39b 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -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 diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index cb2b922a9..fe1544b61 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -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 diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index d609fa1ca..52b6dfea0 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -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 !! @@ -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 @@ -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),& diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index c53d155ae..71190d64d 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -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) @@ -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 diff --git a/src/misc/solver_module.f90 b/src/misc/solver_module.f90 index d8fc85c89..d0aa38e62 100644 --- a/src/misc/solver_module.f90 +++ b/src/misc/solver_module.f90 @@ -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 @@ -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 \ No newline at end of file diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index 83b9ca4ce..b58e3c20a 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -248,7 +248,7 @@ module subroutine swiftest_io_display_run_information(self, param, integration_t if (param%display_style == "PROGRESS") then write(pbarmessage,fmt=pbarfmt) self%t, param%tstop - call pbar%update(1,message=pbarmessage) + call pbar%update(1_I8B,message=pbarmessage) else if (param%display_style == "COMPACT") then call self%compact_output(param,integration_timer) end if @@ -2894,7 +2894,7 @@ module subroutine swiftest_io_set_display_param(self, display_style) self%display_unit = OUTPUT_UNIT !! stdout from iso_fortran_env self%log_output = .false. case ('COMPACT', 'PROGRESS') - inquire(SWIFTEST_LOG_FILE, exist=fileExists) + inquire(file=SWIFTEST_LOG_FILE, exist=fileExists) if (self%lrestart.and.fileExists) then open(unit=SWIFTEST_LOG_OUT, file=SWIFTEST_LOG_FILE, status="OLD", position="APPEND", err = 667, iomsg = errmsg) else