diff --git a/app/simple_diag_orbit.f90 b/app/simple_diag_orbit.f90 index de31dfff..4e138248 100644 --- a/app/simple_diag_orbit.f90 +++ b/app/simple_diag_orbit.f90 @@ -19,15 +19,11 @@ program diag_orbit_main implicit none -! Define dp kind parameter -integer, parameter :: dp = kind(1.0d0) - character(256) :: config_file, particle_arg type(tracer_t) :: norb type(symplectic_integrator_t) :: si type(field_can_t) :: field_can integer :: particle_number -real(dp), dimension(5) :: z ! Initialize timing call init_timer() @@ -143,4 +139,4 @@ program diag_orbit_main print *, "Generated detailed trajectory plots showing real particle" print *, "motion in the magnetic field using symplectic integration." -end program diag_orbit_main \ No newline at end of file +end program diag_orbit_main diff --git a/cmake/Util.cmake b/cmake/Util.cmake index 83d5ccae..7dbf7f17 100644 --- a/cmake/Util.cmake +++ b/cmake/Util.cmake @@ -61,14 +61,6 @@ function(find_or_fetch DEPENDENCY) endif() endif() - if(NOT _use_local AND DEFINED ENV{CODE}) - get_filename_component(_candidate "$ENV{CODE}/${_dep_lower}" ABSOLUTE) - if(IS_DIRECTORY "${_candidate}") - set(_source_dir "${_candidate}") - set(_use_local TRUE) - endif() - endif() - if(_use_local) if("${_dep_lower}" STREQUAL "libneo") set(LIBNEO_ENABLE_TESTS OFF CACHE BOOL "" FORCE) diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index b8ac740f..2110387f 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -86,11 +86,7 @@ set(FILES_TO_WRAP ) # Add libneo files to wrap -if(DEFINED ENV{CODE}) - set(LIBNEO_SOURCE_DIR $ENV{CODE}/libneo) -else() - set(LIBNEO_SOURCE_DIR ${libneo_SOURCE_DIR}) -endif() +set(LIBNEO_SOURCE_DIR ${libneo_SOURCE_DIR}) list(APPEND FILES_TO_WRAP ${LIBNEO_SOURCE_DIR}/src/canonical_coordinates_mod.f90) diff --git a/src/boozer_converter.F90 b/src/boozer_converter.F90 index 209e01f8..5d61a4f5 100644 --- a/src/boozer_converter.F90 +++ b/src/boozer_converter.F90 @@ -82,7 +82,7 @@ subroutine get_boozer_coordinates_impl use vector_potentail_mod, only: ns, hs use new_vmec_stuff_mod, only: n_theta, n_phi, h_theta, h_phi, ns_s, ns_tp use boozer_coordinates_mod, only: ns_s_B, ns_tp_B, ns_B, n_theta_B, n_phi_B, & - hs_B, h_theta_B, h_phi_B, use_B_r + hs_B, h_theta_B, h_phi_B implicit none @@ -784,7 +784,7 @@ end subroutine compute_boozer_data !> Compute radial covariant magnetic field B_rho from symmetry flux coordinates subroutine compute_br_from_symflux(rho_tor, aiota_arr, Gfunc, Bcovar_symfl) - use boozer_coordinates_mod, only: ns_B, n_theta_B, n_phi_B + use boozer_coordinates_mod, only: ns_B, n_phi_B use plag_coeff_sub, only: plag_coeff real(dp), intent(in) :: rho_tor(:) @@ -798,7 +798,7 @@ subroutine compute_br_from_symflux(rho_tor, aiota_arr, Gfunc, Bcovar_symfl) integer :: i_rho, i_phi, ibeg, iend, nshift real(dp) :: coef(0:NDER, NPOILAG) - nshift = NPOILAG/2 + nshift = (NPOILAG - 1)/2 do i_rho = 1, ns_B ibeg = i_rho - nshift diff --git a/src/classification.f90 b/src/classification.f90 index 2836bb14..3e29da14 100644 --- a/src/classification.f90 +++ b/src/classification.f90 @@ -1,4 +1,5 @@ module classification +use, intrinsic :: iso_fortran_env, only: dp => real64 use omp_lib use params, only: zstart, zend, times_lost, trap_par, perp_inv, iclass, & ntimstep, confpart_trap, confpart_pass, notrace_passing, contr_pp, & @@ -18,9 +19,6 @@ module classification implicit none - ! Define real(dp) kind parameter - integer, parameter :: dp = kind(1.0d0) - ! Classification result type - separates data from I/O ! Note: 0=unclassified means the classification was not computed ! This depends on orbit type (trapped/passing) and class_plot flag @@ -65,7 +63,7 @@ subroutine trace_orbit_with_classifiers(anorb, ipart, class_result) type(tracer_t), intent(inout) :: anorb integer, intent(in) :: ipart type(classification_result_t), intent(out) :: class_result - integer :: ierr, ierr_coll + integer :: ierr real(dp), dimension(5) :: z real(dp) :: bmod,sqrtg real(dp), dimension(3) :: bder, hcovar, hctrvr, hcurl diff --git a/src/cut_detector.f90 b/src/cut_detector.f90 index 7223a592..b727bd27 100644 --- a/src/cut_detector.f90 +++ b/src/cut_detector.f90 @@ -1,4 +1,5 @@ module cut_detector + use, intrinsic :: iso_fortran_env, only: dp => real64 use util, only: twopi use simple, only: tstep use orbit_symplectic, only: symplectic_integrator_t @@ -9,8 +10,6 @@ module cut_detector implicit none - ! Define real(dp) kind parameter - integer, parameter :: dp = kind(1.0d0) save integer, parameter :: n_tip_vars = 6 @@ -86,33 +85,29 @@ subroutine trace_to_cut(self, si, f, z, var_cut, cut_type, ierr) integer :: i real(dp) :: phiper = 0.0d0 - do i=1, nstep_max + do i = 1, nplagr call tstep(si, f, z, ierr) - if(ierr.ne.0) exit - - self%par_inv = self%par_inv+z(5)**2 ! parallel adiabatic invariant - - if(i.le.nplagr) then !<=first nplagr points to initialize stencil - ! Note: i is guaranteed to be <= nplagr here, so array access is safe - self%orb_sten(1:5,i)=z - self%orb_sten(6,i)=self%par_inv - else !<=normal case, shift stencil - self%orb_sten(1:5,self%ipoi(1))=z - self%orb_sten(6,self%ipoi(1))=self%par_inv - self%ipoi=cshift(self%ipoi,1) - endif + if (ierr.ne.0) exit + + self%par_inv = self%par_inv + z(5)**2 ! parallel adiabatic invariant + + self%orb_sten(1:5, i) = z + self%orb_sten(6, i) = self%par_inv !------------------------------------------------------------------------- ! Tip detection and interpolation - if(self%alam_prev.lt.0.d0.and.z(5).gt.0.d0) self%itip=0 !<=tip has been passed + if(self%alam_prev.lt.0.d0.and.z(5).gt.0.d0) then + self%itip=0 !<=tip has been passed + end if self%itip=self%itip+1 self%alam_prev=z(5) !<=use only initialized stencil if(i.gt.nplagr .and. self%itip.eq.nplagr/2) then cut_type = 0 - call plag_coeff(nplagr, nder, 0d0, self%orb_sten(5, self%ipoi), self%coef) + call plag_coeff(nplagr, nder, 0d0, self%orb_sten(5, self%ipoi), & + self%coef) var_cut = matmul(self%orb_sten(:, self%ipoi), self%coef(0,:)) var_cut(2) = modulo(var_cut(2), twopi) var_cut(3) = modulo(var_cut(3), twopi) @@ -140,7 +135,70 @@ subroutine trace_to_cut(self, si, f, z, var_cut, cut_type, ierr) if(i.gt.nplagr .and. self%iper.eq.nplagr/2) then cut_type = 1 !<=stencil around periodic boundary is complete, interpolate - call plag_coeff(nplagr, nder, 0d0, self%orb_sten(3, self%ipoi) - phiper, self%coef) + call plag_coeff(nplagr, nder, 0d0, self%orb_sten(3, self%ipoi) - & + phiper, self%coef) + var_cut = matmul(self%orb_sten(:, self%ipoi), self%coef(0,:)) + var_cut(2) = modulo(var_cut(2), twopi) + var_cut(3) = modulo(var_cut(3), twopi) + return + end if + + ! End periodic boundary footprint detection and interpolation + !------------------------------------------------------------------------- + + end do + + do i = nplagr + 1, nstep_max + call tstep(si, f, z, ierr) + if (ierr.ne.0) exit + + self%par_inv = self%par_inv + z(5)**2 ! parallel adiabatic invariant + + self%orb_sten(1:5, self%ipoi(1)) = z + self%orb_sten(6, self%ipoi(1)) = self%par_inv + self%ipoi = cshift(self%ipoi, 1) + + !------------------------------------------------------------------------- + ! Tip detection and interpolation + + if(self%alam_prev.lt.0.d0.and.z(5).gt.0.d0) then + self%itip=0 !<=tip has been passed + end if + self%itip=self%itip+1 + self%alam_prev=z(5) + + if(self%itip.eq.nplagr/2) then + cut_type = 0 + call plag_coeff(nplagr, nder, 0d0, self%orb_sten(5, self%ipoi), & + self%coef) + var_cut = matmul(self%orb_sten(:, self%ipoi), self%coef(0,:)) + var_cut(2) = modulo(var_cut(2), twopi) + var_cut(3) = modulo(var_cut(3), twopi) + self%par_inv = self%par_inv - var_cut(6) + return + end if + + ! End tip detection and interpolation + + !------------------------------------------------------------------------- + ! Periodic boundary footprint detection and interpolation + + if(z(3).gt.dble(self%kper+1)*self%fper) then + self%iper=0 !<=periodic boundary has been passed + phiper=dble(self%kper+1)*self%fper + self%kper=self%kper+1 + elseif(z(3).lt.dble(self%kper)*self%fper) then + self%iper=0 !<=periodic boundary has been passed + phiper=dble(self%kper)*self%fper + self%kper=self%kper-1 + endif + self%iper=self%iper+1 + + if(self%iper.eq.nplagr/2) then + cut_type = 1 + !<=stencil around periodic boundary is complete, interpolate + call plag_coeff(nplagr, nder, 0d0, self%orb_sten(3, self%ipoi) - & + phiper, self%coef) var_cut = matmul(self%orb_sten(:, self%ipoi), self%coef(0,:)) var_cut(2) = modulo(var_cut(2), twopi) var_cut(3) = modulo(var_cut(3), twopi) diff --git a/src/diag/diag_orbit.f90 b/src/diag/diag_orbit.f90 index 1bd2512d..fcac6163 100644 --- a/src/diag/diag_orbit.f90 +++ b/src/diag/diag_orbit.f90 @@ -133,7 +133,7 @@ subroutine integrate_orbit_with_trajectory_debug(si, f, particle_number) real(dp), allocatable :: pphi_traj(:) integer, allocatable :: newton_iter_traj(:) real(dp), dimension(5) :: z, xlast - integer :: it, ktau, point_idx, newton_iters, ierr_orbit + integer :: it, ktau, point_idx, newton_iters integer, parameter :: maxit = 32 real(dp) :: current_time integer :: total_points, field_eval_count @@ -269,4 +269,4 @@ subroutine write_trajectory_data(time_traj, s_traj, theta_traj, phi_traj, pphi_t end subroutine write_trajectory_data -end module diag_orbit \ No newline at end of file +end module diag_orbit diff --git a/src/netcdf_results_output.f90 b/src/netcdf_results_output.f90 index 8b68094f..0e372ed2 100644 --- a/src/netcdf_results_output.f90 +++ b/src/netcdf_results_output.f90 @@ -318,7 +318,7 @@ subroutine compute_cartesian_positions(xstart_cart, xend_cart) ! Check if particle was traced: zend is set to exactly 0 for untraced ! particles (see simple_main.f90:493 and classification.f90:94) - zend_is_zero = all(zend(1:3, i) == 0.0_dp) + zend_is_zero = all(abs(zend(1:3, i)) <= 0.0_dp) if (zend_is_zero) then ! Particle not traced - use start position xend_cart(:, i) = xstart_cart(:, i) diff --git a/src/orbit_symplectic.f90 b/src/orbit_symplectic.f90 index bc9d7f5a..b37d37a2 100644 --- a/src/orbit_symplectic.f90 +++ b/src/orbit_symplectic.f90 @@ -1,5 +1,6 @@ module orbit_symplectic +use, intrinsic :: iso_fortran_env, only: dp => real64 use util, only: pi, twopi use field_can_mod, only: field_can_t, get_val, get_derivatives, get_derivatives2, & eval_field => evaluate @@ -15,9 +16,6 @@ module orbit_symplectic implicit none -! Define real(dp) kind parameter -integer, parameter :: dp = kind(1.0d0) - procedure(orbit_timestep_sympl_i), pointer :: orbit_timestep_sympl => null() contains @@ -26,8 +24,6 @@ module orbit_symplectic ! subroutine orbit_sympl_init(si, f, z, dt, ntau, rtol_init, mode_init) ! - use plag_coeff_sub, only : plag_coeff - type(symplectic_integrator_t), intent(inout) :: si type(field_can_t), intent(inout) :: f real(dp), intent(in) :: z(:) @@ -36,8 +32,6 @@ subroutine orbit_sympl_init(si, f, z, dt, ntau, rtol_init, mode_init) real(dp), intent(in) :: rtol_init integer, intent(in) :: mode_init - integer :: k - si%atol = 1d-15 si%rtol = rtol_init @@ -429,7 +423,6 @@ subroutine newton2(si, f, x, atol, rtol, maxit, xlast) real(dp), intent(out) :: xlast(n) real(dp) :: fvec(n), fjac(n,n), jinv(n,n) - integer :: pivot(n), info real(dp) :: xabs(n), tolref(n), fabs(n) real(dp) :: det @@ -792,8 +785,13 @@ subroutine jac_rk_lobatto(si, fs, s, jac) call coeff_rk_lobatto(s, a, ahat, b, c) jac = 0d0 + dHprime = 0.0d0 - Hprime = fs%dH(1)/fs%dpth(1) + Hprime = 0.0d0 + Hprime(1) = fs(1)%dH(1)/fs(1)%dpth(1) + do k = 2, s + Hprime(k) = fs(k)%dH(1)/fs(k)%dpth(1) + end do dHprime(1) = (fs(1)%d2H(1)-Hprime(1)*fs(1)%d2pth(1))/fs(1)%dpth(1) ! d/dr dHprime(2) = (fs(1)%d2H(7)-Hprime(1)*fs(1)%d2pth(7))/fs(1)%dpth(1) ! d/dpph do k = 2, s @@ -1225,7 +1223,7 @@ subroutine orbit_timestep_sympl_expl_impl_euler(si, f, ierr) integer, parameter :: maxit = 32 real(dp), dimension(n) :: x, xlast - integer :: k, ktau + integer :: ktau ierr = 0 ktau = 0 @@ -1284,7 +1282,7 @@ subroutine orbit_timestep_sympl_impl_expl_euler(si, f, ierr) integer, parameter :: maxit = 32 real(dp), dimension(n) :: x, xlast, dz - integer :: k, ktau + integer :: ktau ierr = 0 ktau = 0 @@ -1346,7 +1344,7 @@ subroutine orbit_timestep_sympl_midpoint(si, f, ierr) integer, parameter :: maxit = 8 real(dp), dimension(n) :: x, xlast - integer :: k, ktau + integer :: ktau ierr = 0 ktau = 0 diff --git a/src/samplers.f90 b/src/samplers.f90 index 27fb5040..dd3ec43f 100644 --- a/src/samplers.f90 +++ b/src/samplers.f90 @@ -1,11 +1,9 @@ module samplers + use, intrinsic :: iso_fortran_env, only: dp => real64 use util implicit none - ! Define real(dp) kind parameter - integer, parameter :: dp = kind(1.0d0) - character(len=*), parameter :: START_FILE = 'start.dat' character(len=*), parameter :: START_FILE_ANTS = 'start_ants.dat' character(len=*), parameter :: START_FILE_BATCH = 'batch.dat' @@ -158,7 +156,7 @@ end subroutine sample_surface_fieldline subroutine sample_grid(zstart, grid_density) use params, only: ntestpart, zstart_dim1, zend, times_lost, & - trap_par, perp_inv, iclass, xstart, sbeg + trap_par, perp_inv, iclass, sbeg use util, only: pi real(dp), dimension(:,:), allocatable, intent(inout) :: zstart diff --git a/src/sub_alpha_lifetime_can.f90 b/src/sub_alpha_lifetime_can.f90 index dab75ae6..fcc82164 100644 --- a/src/sub_alpha_lifetime_can.f90 +++ b/src/sub_alpha_lifetime_can.f90 @@ -1,9 +1,8 @@ module alpha_lifetime_sub -implicit none +use, intrinsic :: iso_fortran_env, only: dp => real64 -! Define real(dp) kind parameter -integer, parameter :: dp = kind(1.0d0) +implicit none contains @@ -22,9 +21,8 @@ subroutine elefie_can(x,derphi) ! Output parameters: ! formal: derphi - array of derivatives ! - integer :: ierr - real(dp), dimension(3) :: x,derphi - real(dp) :: r,phi,z,psi,phi_el,phi_el_pr,phi_el_prpr + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: derphi(3) ! derphi=0.d0 ! @@ -76,7 +74,6 @@ subroutine velo_can(tau,z,vz) real(dp) rmumag,rovsqg,rosqgb,rovbm real(dp) a_phi,a_b,a_c,hstar real(dp) s_hc,hpstar,phidot,blodot,bra - real(dp) pardeb ! dimension x(3),bder(3),hcovar(3),hctrvr(3),hcurl(3) dimension derphi(3) @@ -172,7 +169,7 @@ subroutine orbit_timestep_can(z,dtau,dtaumin,relerr,ierr) ! integer, parameter :: ndim=5, nstepmax=1000000 ! - integer :: ierr,j + integer :: ierr real(dp) :: dtau,dtaumin,phi,tau1,tau2 ! real(dp), dimension(2) :: y @@ -385,7 +382,7 @@ subroutine orbit_timestep_axis(z,dtau,dtaumin,relerr,ierr) real(dp), parameter :: snear_axis=0.01d0 ! logical :: near_axis - integer :: ierr,j + integer :: ierr real(dp) :: dtau,dtaumin,phi,tau1,tau2,z1,z2 real(dp) :: relerr ! diff --git a/test/tests/CMakeLists.txt b/test/tests/CMakeLists.txt index fade65d1..66b54402 100644 --- a/test/tests/CMakeLists.txt +++ b/test/tests/CMakeLists.txt @@ -43,8 +43,6 @@ if(SIMPLE_ENABLE_PYTHON_TOOLS) set(LIBNEO_PYTHON_DIR "${libneo_SOURCE_DIR}/python") elseif(DEFINED LIBNEO_SOURCE_DIR AND EXISTS "${LIBNEO_SOURCE_DIR}/python/libneo/chartmap.py") set(LIBNEO_PYTHON_DIR "${LIBNEO_SOURCE_DIR}/python") - elseif(DEFINED ENV{CODE} AND EXISTS "$ENV{CODE}/libneo/python/libneo/chartmap.py") - set(LIBNEO_PYTHON_DIR "$ENV{CODE}/libneo/python") endif() if(NOT LIBNEO_PYTHON_DIR STREQUAL "")