app/spectre_poincare.f90 - GCC Code Coverage Report

Directory:	./
File:	app/spectre_poincare.f90
Date:	2026-10-22 05:02:02

	Coverage	Exec	Total
Lines:	96.2%	126	131
Functions:	100.0%	8	8
Branches:	63.5%	99	156

  
      Line
      Branch
      Exec
      Source
    
      5
      program spectre_poincare
    
          !> Symplectic field-line Poincare sections and rotational-transform profiles
    
          !> for a SPECTRE/SPEC equilibrium, one file per volume.
    
          !>
    
          !> Reads a namelist (default spectre_poincare.in), loads the SPECTRE HDF5
    
          !> equilibrium, seeds field lines per volume at theta = 0.2 with s spread
    
          !> across the volume, and integrates each line with the implicit midpoint
    
          !> map in field_line_spectre. Field periodicity puts every
    
          !> section exactly on zeta0 + k*2*pi/Nfp, so no interpolation is needed.
    
          !>
    
          !> Output per volume: poincare_vol<N>.dat (seed, section, s, theta, R, Z)
    
          !> and iota_vol<N>.dat (seed, fitted rotational transform).
    
          use, intrinsic :: iso_fortran_env, only: dp => real64
    
      5
          use spectre_reader, only: spectre_data_t, load_spectre, free_spectre
    
          use field_line_spectre, only: field_line_psi, field_line_step, spectre_rz, &
    
                                        prepare_field_line_spline, cleanup_field_line_spline
    
          implicit none
    
          real(dp), parameter :: pi = 3.141592653589793238462643383279_dp
    
          real(dp), parameter :: seed_theta = 0.2_dp
    
          real(dp), parameter :: seed_s_min = -0.85_dp
    
          real(dp), parameter :: seed_s_max = 0.85_dp
    
          real(dp), parameter :: s_terminate = 1.0_dp - 1.0e-10_dp
    
          integer, parameter :: max_volumes = 128
    
          character(len=1024) :: spectre_file, seed_file
    
          real(dp) :: zeta0
    
          integer :: nturns, nsteps_per_period, nseeds_per_volume
    
          integer :: spline_n_s, spline_n_theta, spline_n_zeta, spline_order
    
          integer :: volumes(max_volumes)
    
          namelist /poincare/ spectre_file, zeta0, nturns, &
    
              nsteps_per_period, nseeds_per_volume, volumes, seed_file, spline_n_s, &
    
              spline_n_theta, spline_n_zeta, spline_order
    
      5
          type(spectre_data_t) :: data
    
          character(len=1024) :: config_file
    
          integer :: ierr, iv, lvol, n_terminated
    
      5
          integer, allocatable :: vol_list(:)
    
          real(dp) :: h
    
      5
          spectre_file = 'spectre.h5'
    
      5
          seed_file = ''
    
      5
          zeta0 = 0.0_dp
    
      5
          nturns = 200
    
      5
          nsteps_per_period = 256
    
      5
          nseeds_per_volume = 4
    
      5
          spline_n_s = 32
    
      5
          spline_n_theta = 32
    
      5
          spline_n_zeta = 32
    
      5
          spline_order = 3
    
        2/2✓ Branch 0 taken 640 times.
✓ Branch 1 taken 5 times.

      645
          volumes = -1
    
      5
          config_file = 'spectre_poincare.in'
    
        1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.

      5
          if (command_argument_count() >= 1) call get_command_argument(1, config_file)
    
      5
          call read_input(config_file)
    
        1/2✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.

      5
          call load_spectre(trim(spectre_file), data, ierr)
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.

      5
          if (ierr /= 0) error stop 'spectre_poincare: load_spectre failed'
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.

      5
          call select_volumes(data%Mvol, volumes, vol_list)
    
      5
          h = (2.0_dp*pi/real(data%Nfp, dp))/real(nsteps_per_period, dp)
    
      5
          write (*, '(A,I0,A,I0,A,I0)') 'spectre_poincare: Nfp = ', data%Nfp, &
    
      10
              ', Mvol = ', data%Mvol, ', tracing volumes: ', size(vol_list)
    
      5
          n_terminated = 0
    
        2/2✓ Branch 0 taken 7 times.
✓ Branch 1 taken 5 times.

      12
          do iv = 1, size(vol_list)
    
      7
              lvol = vol_list(iv)
    
      12
              call trace_volume(data, lvol, h, n_terminated)
    
          end do
    
      5
          write (*, '(A,I0,A)') 'spectre_poincare: ', n_terminated, &
    
      10
              ' seed(s) terminated at |s| > 1 - 1e-10'
    
      5
          call free_spectre(data)
    
      contains
    
      5
          subroutine read_input(fname)
    
              character(len=*), intent(in) :: fname
    
              integer :: unit, ios
    
              logical :: exists
    
      5
              inquire (file=fname, exist=exists)
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.

      5
              if (.not. exists) then
    
      ✗
                  write (*, '(A)') 'spectre_poincare: input file not found: '//trim(fname)
    
      ✗
                  error stop 'spectre_poincare: missing input file'
    
              end if
    
      5
              open (newunit=unit, file=fname, status='old', action='read')
    
      5
              read (unit, nml=poincare, iostat=ios)
    
      5
              close (unit)
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.

      5
              if (ios /= 0) error stop 'spectre_poincare: failed to read namelist'
    
      5
          end subroutine read_input
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.

      5
          subroutine select_volumes(mvol, requested, list)
    
              integer, intent(in) :: mvol
    
              integer, intent(in) :: requested(:)
    
              integer, allocatable, intent(out) :: list(:)
    
              integer :: i, k, count
    
      5
              count = 0
    
        2/2✓ Branch 0 taken 640 times.
✓ Branch 1 taken 5 times.

      645
              do i = 1, size(requested)
    
        2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 636 times.

      645
                  if (requested(i) >= 1) count = count + 1
    
              end do
    
        2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 4 times.

      5
              if (count == 0) then
    
        7/14✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 time.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 time.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 time.
✗ Branch 14 not taken.
✓ Branch 15 taken 1 time.

      3
                  allocate (list(mvol))
    
        2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3 times.

      4
                  do i = 1, mvol
    
      4
                      list(i) = i
    
                  end do
    
                  return
    
              end if
    
        7/14✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 4 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 4 times.

      12
              allocate (list(count))
    
      4
              k = 0
    
        2/2✓ Branch 0 taken 512 times.
✓ Branch 1 taken 4 times.

      516
              do i = 1, size(requested)
    
        2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 508 times.

      516
                  if (requested(i) >= 1) then
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.

      4
                      if (requested(i) > mvol) then
    
      ✗
                          error stop 'spectre_poincare: requested volume exceeds Mvol'
    
                      end if
    
      4
                      k = k + 1
    
      4
                      list(k) = requested(i)
    
                  end if
    
              end do
    
          end subroutine select_volumes
    
      24
          pure function seed_s(iseed, nseeds) result(s)
    
              integer, intent(in) :: iseed, nseeds
    
              real(dp) :: s
    
      24
              if (nseeds <= 1) then
    
                  s = 0.5_dp*(seed_s_min + seed_s_max)
    
              else
    
                  s = seed_s_min + (seed_s_max - seed_s_min) &
    
      24
                      *real(iseed - 1, dp)/real(nseeds - 1, dp)
    
              end if
    
          end function seed_s
    
      7
          subroutine trace_volume(data, lvol, h, n_terminated)
    
              type(spectre_data_t), intent(in) :: data
    
              integer, intent(in) :: lvol
    
              real(dp), intent(in) :: h
    
              integer, intent(inout) :: n_terminated
    
              integer :: nsec, iseed, ksec, istep, uposec, uiota
    
              real(dp) :: s, theta, zeta, psi, R, Z, iota
    
      7
              real(dp), allocatable :: zeta_sec(:), theta_sec(:), seeds(:, :)
    
              character(len=64) :: poi_name, iota_name
    
              logical :: terminated
    
      7
              nsec = nturns*data%Nfp
    
              call volume_seeds(lvol, seeds)
    
              call prepare_field_line_spline(data, lvol, spline_n_s, spline_n_theta, &
    
      7
                                             spline_n_zeta, spline_order)
    
        10/20✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 7 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 7 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 7 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 7 times.
✓ Branch 16 taken 7 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 7 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 7 times.

      28
              allocate (zeta_sec(0:nsec), theta_sec(0:nsec))
    
      7
              write (poi_name, '(A,I0,A)') 'poincare_vol', lvol, '.dat'
    
      7
              write (iota_name, '(A,I0,A)') 'iota_vol', lvol, '.dat'
    
      7
              open (newunit=uposec, file=trim(poi_name), status='replace', action='write')
    
      7
              open (newunit=uiota, file=trim(iota_name), status='replace', action='write')
    
      7
              write (uposec, '(A)') '# seed  section  s  theta  R  Z'
    
      7
              write (uiota, '(A)') '# seed  iota'
    
        2/2✓ Branch 0 taken 26 times.
✓ Branch 1 taken 7 times.

      33
              do iseed = 1, size(seeds, 1)
    
      26
                  s = seeds(iseed, 1)
    
      26
                  theta = seeds(iseed, 2)
    
      26
                  zeta = zeta0
    
      26
                  psi = field_line_psi(data, lvol, s, theta, zeta)
    
      26
                  terminated = .false.
    
      26
                  call spectre_rz(data, lvol, s, theta, zeta0, R, Z)
    
      26
                  call write_section(uposec, iseed, 0, s, theta, R, Z)
    
      26
                  zeta_sec(0) = zeta
    
      26
                  theta_sec(0) = theta
    
        2/2✓ Branch 0 taken 3790 times.
✓ Branch 1 taken 26 times.

      3816
                  do ksec = 1, nsec
    
        2/2✓ Branch 0 taken 1029760 times.
✓ Branch 1 taken 3790 times.

      1033550
                      do istep = 1, nsteps_per_period
    
      1029760
                          call field_line_step(data, lvol, h, s, theta, zeta, psi)
    
        1/2✓ Branch 0 taken 1029760 times.
✗ Branch 1 not taken.

      1033550
                          if (abs(s) > s_terminate) then
    
                              terminated = .true.
    
                              exit
    
                          end if
    
                      end do
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3790 times.

      3790
                      if (terminated) exit
    
      3790
                      call spectre_rz(data, lvol, s, theta, zeta0, R, Z)
    
      3790
                      call write_section(uposec, iseed, ksec, s, theta, R, Z)
    
      3790
                      zeta_sec(ksec) = zeta
    
      7606
                      theta_sec(ksec) = theta
    
                  end do
    
        1/2✗ Branch 1 not taken.
✓ Branch 2 taken 26 times.

      59
                  if (terminated) then
    
      ✗
                      n_terminated = n_terminated + 1
    
      ✗
                      write (uiota, '(I6,1x,A)') iseed, 'terminated'
    
                  else
    
      26
                      iota = fit_iota(zeta_sec, theta_sec)
    
      26
                      write (uiota, '(I6,1x,ES23.15)') iseed, iota
    
                  end if
    
              end do
    
      7
              close (uposec)
    
      7
              close (uiota)
    
        2/4✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 7 times.

      7
              deallocate (zeta_sec, theta_sec)
    
      7
              call cleanup_field_line_spline
    
        1/2✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.

      14
          end subroutine trace_volume
    
      7
          subroutine volume_seeds(lvol, seeds)
    
              integer, intent(in) :: lvol
    
              real(dp), allocatable, intent(out) :: seeds(:, :)
    
              integer :: unit, ios, ivol, count, i
    
              real(dp) :: sval, theta
    
        2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 1 time.

      7
              if (len_trim(seed_file) == 0) then
    
        7/14✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 6 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 6 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 6 times.

      18
                  allocate (seeds(nseeds_per_volume, 2))
    
        2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 6 times.

      30
                  do i = 1, nseeds_per_volume
    
        3/4✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 48 times.
✓ Branch 3 taken 24 times.

      102
                      seeds(i, :) = [seed_s(i, nseeds_per_volume), seed_theta]
    
                  end do
    
      6
                  return
    
              end if
    
      1
              open (newunit=unit, file=trim(seed_file), status='old', action='read')
    
      1
              count = 0
    
              do
    
      3
                  read (unit, *, iostat=ios) ivol, sval, theta
    
        2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.

      3
                  if (ios /= 0) exit
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.

      2
                  if (ivol == lvol - 1) count = count + 1
    
              end do
    
      1
              rewind (unit)
    
        7/14✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 time.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 time.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 time.
✗ Branch 14 not taken.
✓ Branch 15 taken 1 time.

      3
              allocate (seeds(count, 2))
    
      1
              i = 0
    
              do
    
      3
                  read (unit, *, iostat=ios) ivol, sval, theta
    
        2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.

      3
                  if (ios /= 0) exit
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.

      3
                  if (ivol == lvol - 1) then
    
      2
                      i = i + 1
    
        2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.

      6
                      seeds(i, :) = [sval, theta]
    
                  end if
    
              end do
    
      1
              close (unit)
    
          end subroutine volume_seeds
    
      3816
          subroutine write_section(unit, iseed, ksec, s, theta, R, Z)
    
              integer, intent(in) :: unit, iseed, ksec
    
              real(dp), intent(in) :: s, theta, R, Z
    
      3816
              write (unit, '(I6,1x,I8,4(1x,ES23.15))') iseed, ksec, s, theta, R, Z
    
      3816
          end subroutine write_section
    
        2/4✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 26 times.

      26
          pure function fit_iota(zeta_sec, theta_sec) result(iota)
    
              real(dp), intent(in) :: zeta_sec(:), theta_sec(:)
    
              real(dp) :: iota
    
              integer :: n
    
              real(dp) :: zbar, tbar, szz, szt
    
      26
              n = size(zeta_sec)
    
        2/2✓ Branch 0 taken 3816 times.
✓ Branch 1 taken 26 times.

      3842
              zbar = sum(zeta_sec)/real(n, dp)
    
        2/2✓ Branch 0 taken 3816 times.
✓ Branch 1 taken 26 times.

      3842
              tbar = sum(theta_sec)/real(n, dp)
    
        2/2✓ Branch 0 taken 26 times.
✓ Branch 1 taken 3816 times.

      3842
              szz = sum((zeta_sec - zbar)**2)
    
        2/2✓ Branch 0 taken 3816 times.
✓ Branch 1 taken 26 times.

      3842
              szt = sum((zeta_sec - zbar)*(theta_sec - tbar))
    
      26
              iota = szt/szz
    
      26
          end function fit_iota
    
      end program spectre_poincare

Line	Branch	Exec	Source
1		5	program spectre_poincare
2			!> Symplectic field-line Poincare sections and rotational-transform profiles
3			!> for a SPECTRE/SPEC equilibrium, one file per volume.
4			!>
5			!> Reads a namelist (default spectre_poincare.in), loads the SPECTRE HDF5
6			!> equilibrium, seeds field lines per volume at theta = 0.2 with s spread
7			!> across the volume, and integrates each line with the implicit midpoint
8			!> map in field_line_spectre. Field periodicity puts every
9			!> section exactly on zeta0 + k2pi/Nfp, so no interpolation is needed.
10			!>
11			!> Output per volume: poincare_vol<N>.dat (seed, section, s, theta, R, Z)
12			!> and iota_vol<N>.dat (seed, fitted rotational transform).
13			use, intrinsic :: iso_fortran_env, only: dp => real64
14		5	use spectre_reader, only: spectre_data_t, load_spectre, free_spectre
15			use field_line_spectre, only: field_line_psi, field_line_step, spectre_rz, &
16			prepare_field_line_spline, cleanup_field_line_spline
17
18			implicit none
19
20			real(dp), parameter :: pi = 3.141592653589793238462643383279_dp
21			real(dp), parameter :: seed_theta = 0.2_dp
22			real(dp), parameter :: seed_s_min = -0.85_dp
23			real(dp), parameter :: seed_s_max = 0.85_dp
24			real(dp), parameter :: s_terminate = 1.0_dp - 1.0e-10_dp
25			integer, parameter :: max_volumes = 128
26
27			character(len=1024) :: spectre_file, seed_file
28			real(dp) :: zeta0
29			integer :: nturns, nsteps_per_period, nseeds_per_volume
30			integer :: spline_n_s, spline_n_theta, spline_n_zeta, spline_order
31			integer :: volumes(max_volumes)
32			namelist /poincare/ spectre_file, zeta0, nturns, &
33			nsteps_per_period, nseeds_per_volume, volumes, seed_file, spline_n_s, &
34			spline_n_theta, spline_n_zeta, spline_order
35
36		5	type(spectre_data_t) :: data
37			character(len=1024) :: config_file
38			integer :: ierr, iv, lvol, n_terminated
39		5	integer, allocatable :: vol_list(:)
40			real(dp) :: h
41
42		5	spectre_file = 'spectre.h5'
43		5	seed_file = ''
44		5	zeta0 = 0.0_dp
45		5	nturns = 200
46		5	nsteps_per_period = 256
47		5	nseeds_per_volume = 4
48		5	spline_n_s = 32
49		5	spline_n_theta = 32
50		5	spline_n_zeta = 32
51		5	spline_order = 3
52	2/2 ✓ Branch 0 taken 640 times. ✓ Branch 1 taken 5 times.	645	volumes = -1
53
54		5	config_file = 'spectre_poincare.in'
55	1/2 ✓ Branch 1 taken 5 times. ✗ Branch 2 not taken.	5	if (command_argument_count() >= 1) call get_command_argument(1, config_file)
56
57		5	call read_input(config_file)
58	1/2 ✓ Branch 2 taken 5 times. ✗ Branch 3 not taken.	5	call load_spectre(trim(spectre_file), data, ierr)
59	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 5 times.	5	if (ierr /= 0) error stop 'spectre_poincare: load_spectre failed'
60
61	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 5 times.	5	call select_volumes(data%Mvol, volumes, vol_list)
62		5	h = (2.0_dp*pi/real(data%Nfp, dp))/real(nsteps_per_period, dp)
63
64		5	write (*, '(A,I0,A,I0,A,I0)') 'spectre_poincare: Nfp = ', data%Nfp, &
65		10	', Mvol = ', data%Mvol, ', tracing volumes: ', size(vol_list)
66
67		5	n_terminated = 0
68	2/2 ✓ Branch 0 taken 7 times. ✓ Branch 1 taken 5 times.	12	do iv = 1, size(vol_list)
69		7	lvol = vol_list(iv)
70		12	call trace_volume(data, lvol, h, n_terminated)
71			end do
72
73		5	write (*, '(A,I0,A)') 'spectre_poincare: ', n_terminated, &
74		10	' seed(s) terminated at \|s\| > 1 - 1e-10'
75
76		5	call free_spectre(data)
77
78			contains
79
80		5	subroutine read_input(fname)
81			character(len=*), intent(in) :: fname
82
83			integer :: unit, ios
84			logical :: exists
85
86		5	inquire (file=fname, exist=exists)
87	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 5 times.	5	if (.not. exists) then
88		✗	write (*, '(A)') 'spectre_poincare: input file not found: '//trim(fname)
89		✗	error stop 'spectre_poincare: missing input file'
90			end if
91
92		5	open (newunit=unit, file=fname, status='old', action='read')
93		5	read (unit, nml=poincare, iostat=ios)
94		5	close (unit)
95	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 5 times.	5	if (ios /= 0) error stop 'spectre_poincare: failed to read namelist'
96		5	end subroutine read_input
97
98	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 5 times.	5	subroutine select_volumes(mvol, requested, list)
99			integer, intent(in) :: mvol
100			integer, intent(in) :: requested(:)
101			integer, allocatable, intent(out) :: list(:)
102
103			integer :: i, k, count
104
105		5	count = 0
106	2/2 ✓ Branch 0 taken 640 times. ✓ Branch 1 taken 5 times.	645	do i = 1, size(requested)
107	2/2 ✓ Branch 0 taken 4 times. ✓ Branch 1 taken 636 times.	645	if (requested(i) >= 1) count = count + 1
108			end do
109
110	2/2 ✓ Branch 0 taken 1 time. ✓ Branch 1 taken 4 times.	5	if (count == 0) then
111	7/14 ✗ Branch 0 not taken. ✓ Branch 1 taken 1 time. ✓ Branch 2 taken 1 time. ✗ Branch 3 not taken. ✓ Branch 4 taken 1 time. ✗ Branch 5 not taken. ✓ Branch 6 taken 1 time. ✗ Branch 7 not taken. ✗ Branch 8 not taken. ✓ Branch 9 taken 1 time. ✗ Branch 11 not taken. ✓ Branch 12 taken 1 time. ✗ Branch 14 not taken. ✓ Branch 15 taken 1 time.	3	allocate (list(mvol))
112	2/2 ✓ Branch 0 taken 1 time. ✓ Branch 1 taken 3 times.	4	do i = 1, mvol
113		4	list(i) = i
114			end do
115			return
116			end if
117
118	7/14 ✗ Branch 0 not taken. ✓ Branch 1 taken 4 times. ✓ Branch 2 taken 4 times. ✗ Branch 3 not taken. ✓ Branch 4 taken 4 times. ✗ Branch 5 not taken. ✓ Branch 6 taken 4 times. ✗ Branch 7 not taken. ✗ Branch 8 not taken. ✓ Branch 9 taken 4 times. ✗ Branch 11 not taken. ✓ Branch 12 taken 4 times. ✗ Branch 14 not taken. ✓ Branch 15 taken 4 times.	12	allocate (list(count))
119		4	k = 0
120	2/2 ✓ Branch 0 taken 512 times. ✓ Branch 1 taken 4 times.	516	do i = 1, size(requested)
121	2/2 ✓ Branch 0 taken 4 times. ✓ Branch 1 taken 508 times.	516	if (requested(i) >= 1) then
122	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 4 times.	4	if (requested(i) > mvol) then
123		✗	error stop 'spectre_poincare: requested volume exceeds Mvol'
124			end if
125		4	k = k + 1
126		4	list(k) = requested(i)
127			end if
128			end do
129			end subroutine select_volumes
130
131		24	pure function seed_s(iseed, nseeds) result(s)
132			integer, intent(in) :: iseed, nseeds
133			real(dp) :: s
134
135		24	if (nseeds <= 1) then
136			s = 0.5_dp*(seed_s_min + seed_s_max)
137			else
138			s = seed_s_min + (seed_s_max - seed_s_min) &
139		24	*real(iseed - 1, dp)/real(nseeds - 1, dp)
140			end if
141			end function seed_s
142
143		7	subroutine trace_volume(data, lvol, h, n_terminated)
144			type(spectre_data_t), intent(in) :: data
145			integer, intent(in) :: lvol
146			real(dp), intent(in) :: h
147			integer, intent(inout) :: n_terminated
148
149			integer :: nsec, iseed, ksec, istep, uposec, uiota
150			real(dp) :: s, theta, zeta, psi, R, Z, iota
151		7	real(dp), allocatable :: zeta_sec(:), theta_sec(:), seeds(:, :)
152			character(len=64) :: poi_name, iota_name
153			logical :: terminated
154
155		7	nsec = nturns*data%Nfp
156			call volume_seeds(lvol, seeds)
157			call prepare_field_line_spline(data, lvol, spline_n_s, spline_n_theta, &
158		7	spline_n_zeta, spline_order)
159	10/20 ✗ Branch 0 not taken. ✓ Branch 1 taken 7 times. ✓ Branch 2 taken 7 times. ✗ Branch 3 not taken. ✓ Branch 4 taken 7 times. ✗ Branch 5 not taken. ✓ Branch 6 taken 7 times. ✗ Branch 7 not taken. ✗ Branch 8 not taken. ✓ Branch 9 taken 7 times. ✗ Branch 11 not taken. ✓ Branch 12 taken 7 times. ✗ Branch 14 not taken. ✓ Branch 15 taken 7 times. ✓ Branch 16 taken 7 times. ✗ Branch 17 not taken. ✗ Branch 18 not taken. ✓ Branch 19 taken 7 times. ✗ Branch 21 not taken. ✓ Branch 22 taken 7 times.	28	allocate (zeta_sec(0:nsec), theta_sec(0:nsec))
160
161		7	write (poi_name, '(A,I0,A)') 'poincare_vol', lvol, '.dat'
162		7	write (iota_name, '(A,I0,A)') 'iota_vol', lvol, '.dat'
163		7	open (newunit=uposec, file=trim(poi_name), status='replace', action='write')
164		7	open (newunit=uiota, file=trim(iota_name), status='replace', action='write')
165		7	write (uposec, '(A)') '# seed section s theta R Z'
166		7	write (uiota, '(A)') '# seed iota'
167
168	2/2 ✓ Branch 0 taken 26 times. ✓ Branch 1 taken 7 times.	33	do iseed = 1, size(seeds, 1)
169		26	s = seeds(iseed, 1)
170		26	theta = seeds(iseed, 2)
171		26	zeta = zeta0
172		26	psi = field_line_psi(data, lvol, s, theta, zeta)
173		26	terminated = .false.
174
175		26	call spectre_rz(data, lvol, s, theta, zeta0, R, Z)
176		26	call write_section(uposec, iseed, 0, s, theta, R, Z)
177		26	zeta_sec(0) = zeta
178		26	theta_sec(0) = theta
179
180	2/2 ✓ Branch 0 taken 3790 times. ✓ Branch 1 taken 26 times.	3816	do ksec = 1, nsec
181	2/2 ✓ Branch 0 taken 1029760 times. ✓ Branch 1 taken 3790 times.	1033550	do istep = 1, nsteps_per_period
182		1029760	call field_line_step(data, lvol, h, s, theta, zeta, psi)
183	1/2 ✓ Branch 0 taken 1029760 times. ✗ Branch 1 not taken.	1033550	if (abs(s) > s_terminate) then
184			terminated = .true.
185			exit
186			end if
187			end do
188	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 3790 times.	3790	if (terminated) exit
189		3790	call spectre_rz(data, lvol, s, theta, zeta0, R, Z)
190		3790	call write_section(uposec, iseed, ksec, s, theta, R, Z)
191		3790	zeta_sec(ksec) = zeta
192		7606	theta_sec(ksec) = theta
193			end do
194
195	1/2 ✗ Branch 1 not taken. ✓ Branch 2 taken 26 times.	59	if (terminated) then
196		✗	n_terminated = n_terminated + 1
197		✗	write (uiota, '(I6,1x,A)') iseed, 'terminated'
198			else
199		26	iota = fit_iota(zeta_sec, theta_sec)
200		26	write (uiota, '(I6,1x,ES23.15)') iseed, iota
201			end if
202			end do
203
204		7	close (uposec)
205		7	close (uiota)
206	2/4 ✗ Branch 0 not taken. ✓ Branch 1 taken 7 times. ✗ Branch 3 not taken. ✓ Branch 4 taken 7 times.	7	deallocate (zeta_sec, theta_sec)
207		7	call cleanup_field_line_spline
208	1/2 ✓ Branch 1 taken 7 times. ✗ Branch 2 not taken.	14	end subroutine trace_volume
209
210		7	subroutine volume_seeds(lvol, seeds)
211			integer, intent(in) :: lvol
212			real(dp), allocatable, intent(out) :: seeds(:, :)
213
214			integer :: unit, ios, ivol, count, i
215			real(dp) :: sval, theta
216
217	2/2 ✓ Branch 0 taken 6 times. ✓ Branch 1 taken 1 time.	7	if (len_trim(seed_file) == 0) then
218	7/14 ✗ Branch 0 not taken. ✓ Branch 1 taken 6 times. ✓ Branch 2 taken 6 times. ✗ Branch 3 not taken. ✓ Branch 4 taken 6 times. ✗ Branch 5 not taken. ✓ Branch 6 taken 6 times. ✗ Branch 7 not taken. ✗ Branch 8 not taken. ✓ Branch 9 taken 6 times. ✗ Branch 11 not taken. ✓ Branch 12 taken 6 times. ✗ Branch 14 not taken. ✓ Branch 15 taken 6 times.	18	allocate (seeds(nseeds_per_volume, 2))
219	2/2 ✓ Branch 0 taken 24 times. ✓ Branch 1 taken 6 times.	30	do i = 1, nseeds_per_volume
220	3/4 ✓ Branch 0 taken 24 times. ✗ Branch 1 not taken. ✓ Branch 2 taken 48 times. ✓ Branch 3 taken 24 times.	102	seeds(i, :) = [seed_s(i, nseeds_per_volume), seed_theta]
221			end do
222		6	return
223			end if
224
225		1	open (newunit=unit, file=trim(seed_file), status='old', action='read')
226		1	count = 0
227			do
228		3	read (unit, *, iostat=ios) ivol, sval, theta
229	2/2 ✓ Branch 0 taken 2 times. ✓ Branch 1 taken 1 time.	3	if (ios /= 0) exit
230	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 2 times.	2	if (ivol == lvol - 1) count = count + 1
231			end do
232		1	rewind (unit)
233	7/14 ✗ Branch 0 not taken. ✓ Branch 1 taken 1 time. ✓ Branch 2 taken 1 time. ✗ Branch 3 not taken. ✓ Branch 4 taken 1 time. ✗ Branch 5 not taken. ✓ Branch 6 taken 1 time. ✗ Branch 7 not taken. ✗ Branch 8 not taken. ✓ Branch 9 taken 1 time. ✗ Branch 11 not taken. ✓ Branch 12 taken 1 time. ✗ Branch 14 not taken. ✓ Branch 15 taken 1 time.	3	allocate (seeds(count, 2))
234		1	i = 0
235			do
236		3	read (unit, *, iostat=ios) ivol, sval, theta
237	2/2 ✓ Branch 0 taken 2 times. ✓ Branch 1 taken 1 time.	3	if (ios /= 0) exit
238	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 2 times.	3	if (ivol == lvol - 1) then
239		2	i = i + 1
240	2/2 ✓ Branch 0 taken 4 times. ✓ Branch 1 taken 2 times.	6	seeds(i, :) = [sval, theta]
241			end if
242			end do
243		1	close (unit)
244			end subroutine volume_seeds
245
246		3816	subroutine write_section(unit, iseed, ksec, s, theta, R, Z)
247			integer, intent(in) :: unit, iseed, ksec
248			real(dp), intent(in) :: s, theta, R, Z
249
250		3816	write (unit, '(I6,1x,I8,4(1x,ES23.15))') iseed, ksec, s, theta, R, Z
251		3816	end subroutine write_section
252
253	2/4 ✗ Branch 0 not taken. ✓ Branch 1 taken 26 times. ✗ Branch 2 not taken. ✓ Branch 3 taken 26 times.	26	pure function fit_iota(zeta_sec, theta_sec) result(iota)
254			real(dp), intent(in) :: zeta_sec(:), theta_sec(:)
255			real(dp) :: iota
256
257			integer :: n
258			real(dp) :: zbar, tbar, szz, szt
259
260		26	n = size(zeta_sec)
261	2/2 ✓ Branch 0 taken 3816 times. ✓ Branch 1 taken 26 times.	3842	zbar = sum(zeta_sec)/real(n, dp)
262	2/2 ✓ Branch 0 taken 3816 times. ✓ Branch 1 taken 26 times.	3842	tbar = sum(theta_sec)/real(n, dp)
263	2/2 ✓ Branch 0 taken 26 times. ✓ Branch 1 taken 3816 times.	3842	szz = sum((zeta_sec - zbar)**2)
264	2/2 ✓ Branch 0 taken 3816 times. ✓ Branch 1 taken 26 times.	3842	szt = sum((zeta_sec - zbar)*(theta_sec - tbar))
265		26	iota = szt/szz
266		26	end function fit_iota
267
268			end program spectre_poincare
269

Function (Line)	Call count	Line coverage	Branch coverage	Block coverage
MAIN__ (line 1)	called 5 times	100.0%	-%	95.0%
fit_iota.0 (line 253)	called 26 times	100.0%	83.3%	88.0%
main (line 14)	called 5 times	100.0%	66.7%	100.0%
read_input.6 (line 80)	called 5 times	80.0%	50.0%	77.0%
select_volumes.5 (line 98)	called 5 times	94.7%	63.6%	84.0%
trace_volume.3 (line 143)	called 7 times	95.0%	57.9%	85.0%
volume_seeds.2 (line 210)	called 7 times	100.0%	63.0%	91.0%
write_section.1 (line 246)	called 3816 times	100.0%	-%	100.0%