src/gray_core: improve error reporting

- Avoid logging the same error over and over

- Make all the gray_errors actually warnings

- Replace `large_npl` error with `unstable_beam`, which is actually
  the root cause of the former

- Use the gray_main error as exit code

Makefile

@@ -132,7 +132,7 @@ check: $(shell python -Bm tests --list-cases)
 
 # Run a test case
 tests.%: $(GRAY)
-	python -Bm tests '$@' --binary '$(GRAY)'
+	python -Bm tests '$@' --binary '$(GRAY)' --buffer
 
 # Install libraries, binaries and documentation
 install-bin: $(BINARIES) $(LIBRARIES)

src/gray_cli.f90

@@ -64,7 +64,7 @@ contains
     print '(a)', '*Exit status*'
     print '(a)', ' 0  if OK,'
     print '(a)', ' 1  if technical problems   (e.g., invalid arguments, missing input files),'
-    print '(a)', ' 2  if serious trouble      (e.g., computation failed).'
+    print '(a)', ' >1 if simulation trouble   (e.g., unstable gradient, negative absorption).'
     print '(a)', ''
     print '(a)', 'Full documentation is available at'
     print '(a)', '<file://'//PREFIX//'/share/doc/manual.html>'

src/gray_core.f90

@@ -11,13 +11,14 @@ contains
     use const_and_precisions, only : zero, one, comp_tiny
     use polarization, only : ellipse_to_field
     use types,        only : contour, table, wrap
-    use gray_params,  only : gray_parameters, gray_results, EQ_VACUUM
+    use gray_params,  only : gray_parameters, gray_results, EQ_VACUUM, ABSORP_OFF
     use gray_equil,   only : abstract_equil
     use gray_plasma,  only : abstract_plasma
     use gray_tables,  only : init_tables, store_ec_profiles, store_ray_data, &
                              store_beam_shape, find_flux_surfaces, &
                              kinetic_profiles, ec_resonance, inputs_maps
-    use gray_errors,  only : is_critical, print_err_raytracing, print_err_ecrh_cd
+    use gray_errors,  only : gray_error, is_critical, has_new_errors, &
+                             print_err_raytracing, print_err_ecrh_cd
     use beams,        only : xgygcoeff, launchangles2n
     use beamdata,     only : pweight
     use magsurf_data, only : compute_flux_averages, dealloc_surfvec
@@ -34,13 +35,11 @@ contains
     class(abstract_plasma),   intent(in)    :: plasma
     type(contour),            intent(in)    :: limiter
     type(gray_results),       intent(out)   :: results
+    integer(kind=gray_error), intent(out)   :: error
 
     ! Predefined grid for the output profiles (optional)
     real(wp_), dimension(:), intent(in), optional :: rhout
 
-    ! Exit code
-    integer, intent(out) :: error
-
     ! local variables
     real(wp_), parameter :: taucr = 12._wp_, etaucr = exp(-taucr)
     character, dimension(2), parameter :: mode=['O','X']
@@ -53,11 +52,15 @@ contains
     real(wp_) :: rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx
     real(wp_) :: pabs_beam,icd_beam,cpl_beam1,cpl_beam2,cpl_cbeam1,cpl_cbeam2
 
-    integer :: iox,nharm,nhf,iokhawa,ierrn,ierrhcd, index_rt, parent_index_rt
+    integer :: iox,nharm,nhf,iokhawa,index_rt, parent_index_rt
     integer :: ip,ib,iopmin,ipar,child_index_rt
     integer :: igrad_b,istop_pass,nbeam_pass
     logical :: ins_pl,ins_wl,ent_pl,ext_pl,ent_wl,ext_wl,iboff
 
+    ! Error handing
+    integer(kind=gray_error) :: prev_error ! from previous step
+    integer(kind=gray_error) :: curr_error ! current step
+
     ! i: integration step, jk: global ray index
     integer :: i, jk
 
@@ -166,6 +169,9 @@ contains
     psipv(0) = params%antenna%psi
     chipv(0) = params%antenna%chi
 
+    ! Error value for the whole simulation
+    error = 0
+
     nbeam_pass=1                                                                           ! max n of beam per pass
     index_rt=0                                                                             ! global beam index:    1,O      2,X    1st pass
                                                                                            !                       | |      | |
@@ -265,7 +271,14 @@ contains
         ! ======= propagation loop  BEGIN =======
         call log_debug(' propagation loop start', mod='gray_core', proc='gray_main')
 
+        ! error value for the this step
+        curr_error = 0
+
         propagation_loop: do i=1,params%raytracing%nstep
+          ! rotate current and previous step errors
+          prev_error = curr_error
+          curr_error = 0
+
           ! advance one step with "frozen" grad(S_I)
           do jk=1,params%raytracing%nray
             if(iwait(jk)) cycle                                                            ! jk ray is waiting for next pass
@@ -274,9 +287,13 @@ contains
                         ypw(:,jk), gri(:,jk), ggri(:,:,jk), igrad_b)
           end do
           ! update position and grad
-          if(igrad_b == 1) call gradi_upd(params%raytracing, yw, ak0, xc, du1, gri, ggri)
+          if (igrad_b == 1) call gradi_upd(params%raytracing, yw, ak0, xc, &
+                                           du1, gri, ggri, curr_error)
+
+          if (has_new_errors(prev_error, curr_error)) then
+            call print_err_raytracing(curr_error, i, jk)
+          end if
 
-          error = 0
           iopmin = 10
 
           ! =========== rays loop BEGIN ===========
@@ -290,12 +307,7 @@ contains
                             xv, anv, gri(:,jk), ggri(:,:,jk), sox, bres,   &
                             xgcn,ypw(:,jk), psinv, dens, btot, bv, xg, yg, &
                             derxg, anpl, anpr, ddr, ddi, dersdst, derdnm,  &
-                            ierrn, igrad_b)
+                            igrad_b)
-            ! update global error code and print message
-            if(ierrn/=0) then
-              error = ior(error,ierrn)
-              call print_err_raytracing(ierrn, i, anpl)
-            end if
 
             ! check entrance/exit plasma/wall
             zzm = xv(3)*0.01_wp_
@@ -392,11 +404,7 @@ contains
                               xv, anv, gri(:,jk), ggri(:,:,jk), sox, bres, &
                               xgcn, ypw(:,jk), psinv, dens, btot, bv, xg,  &
                               yg, derxg, anpl, anpr, ddr, ddi, dersdst,    &
-                              derdnm, ierrn, igrad_b)                                      ! * update derivatives after reflection
+                              derdnm, igrad_b)                                             ! * update derivatives after reflection
-              if(ierrn/=0) then                                                            ! * update global error code and print message
-                error = ior(error,ierrn)
-                call print_err_raytracing(ierrn, i, anpl)
-              end if
 
               if(jk == 1 .and. ip == 1) then                                               ! * 1st pass, polarization angles at reflection for central ray
                 psipv(index_rt) = psipol
@@ -444,7 +452,7 @@ contains
             ! Compute ECRH&CD if (inside plasma & power available>0 & ray still active)
             ! Note: power check is τ + τ₀ + log(coupling O/X) < τ_critic
 
-            if(ierrn==0 .and. params%ecrh_cd%iwarm>0 .and. ins_pl .and. &
+            if (params%ecrh_cd%iwarm /= ABSORP_OFF .and. ins_pl .and. &
               (tau1(jk)+tau0(jk)+lgcpl1(jk))<=taucr .and. .not.iwait(jk)) then             ! H&CD computation check
               tekev = plasma%temp(psinv)
               block
@@ -452,12 +460,12 @@ contains
                 call alpha_effj(params%ecrh_cd, equil, plasma,            &
                                 psinv, xg, yg, dens, tekev, ak0, bres,    &
                                 derdnm, anpl, anpr, sox, Npr_warm, alpha, &
-                                didp, nharm, nhf, iokhawa, ierrhcd)
+                                didp, nharm, nhf, iokhawa, curr_error)
                 anprre = Npr_warm%re
                 anprim = Npr_warm%im
-                if(ierrhcd /= 0) then
+
-                  error = ior(error, ierrhcd)
+                if (has_new_errors(prev_error, curr_error)) then
-                  call print_err_ecrh_cd(ierrhcd, i, Npr_warm, alpha)
+                  call print_err_ecrh_cd(curr_error, i, jk, Npr_warm, alpha)
                 end if
               end block
             else
@@ -502,6 +510,9 @@ contains
                  nharm, nhf, iokhawa, index_rt, ddr, ddi, xg, yg, derxg)                   ! p0ray/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes)
             end if
 
+          ! Accumulate errors from the latest step
+          error = ior(error, curr_error)
+
           end do rays_loop
           ! ============ rays loop END ============
 
@@ -520,7 +531,7 @@ contains
           ! test whether further trajectory integration is unnecessary
           call vmaxmin(tau1+tau0+lgcpl1, taumn, taumx)                                     ! test on tau + coupling
 
-          if(is_critical(error)) then                                                      ! stop propagation for current beam
+          if (is_critical(curr_error)) then                                                ! stop propagation for current beam
             istop_pass = istop_pass +1                                                     ! * +1 non propagating beam
             if(ip < params%raytracing%ipass) call turnoffray(0,ip,npass,ib,iroff)          ! * de-activate derived beams
             exit
@@ -960,10 +971,9 @@ contains
   subroutine ywppla_upd(params, equil, plasma,                      &
                         xv, anv, dgr, ddgr, sox, bres, xgcn, dery,  &
                         psinv, dens, btot, bv, xg, yg, derxg, anpl, &
-                        anpr, ddr, ddi, dersdst, derdnm, error, igrad)
+                        anpr, ddr, ddi, dersdst, derdnm, igrad)
     ! Compute right-hand side terms of the ray equations (dery)
     ! used after full R-K step and grad(S_I) update
-    use gray_errors, only : raise_error, large_npl
     use gray_params, only : gray_parameters
     use gray_equil,  only : abstract_equil
     use gray_plasma, only : abstract_plasma
@@ -982,33 +992,24 @@ contains
     real(wp_), intent(out) :: psinv, dens, btot, xg, yg, anpl, anpr
     real(wp_), intent(out) :: ddr, ddi, dersdst, derdnm
     real(wp_), intent(out) :: bv(3)
-    integer,   intent(out) :: error
     real(wp_), intent(out) :: derxg(3)
     integer,   intent(in)  :: igrad
 
     ! local variables
-    real(wp_), dimension(3) :: deryg
+    real(wp_):: deryg(3)
-    real(wp_), dimension(3,3) :: derbv
+    real(wp_):: derbv(3,3)
-    real(wp_), parameter :: anplth1 = 0.99_wp_, anplth2 = 1.05_wp_
 
     call plas_deriv(equil, plasma,  xv, bres, xgcn, dens, btot, &
                     bv, derbv, xg, yg, derxg, deryg, psinv)
     call disp_deriv(params, anv, sox, xg, yg, derxg, deryg, bv, derbv, dgr, &
                     ddgr, igrad, dery, anpl, anpr, ddr, ddi, dersdst, derdnm)
-
-    if (abs(anpl) > anplth2) then
-      error = raise_error(error, large_npl, 1)
-    else if (abs(anpl) > anplth1) then
-      error = raise_error(error, large_npl, 0)
-    else
-      error = 0
-    end if
   end subroutine ywppla_upd
 
 
-  subroutine gradi_upd(params, ywrk,ak0,xc,du1,gri,ggri)
+  subroutine gradi_upd(params, ywrk, ak0, xc, du1, gri, ggri, error)
     use const_and_precisions, only : zero, half
     use gray_params,          only : raytracing_parameters
+    use gray_errors,          only : gray_error, unstable_beam, raise_error
 
     ! subroutine parameters
     type(raytracing_parameters),                        intent(in)    :: params
@@ -1017,6 +1018,7 @@ contains
     real(wp_), dimension(3,params%nrayth,params%nrayr), intent(inout) :: xc, du1
     real(wp_), dimension(3,params%nray),                intent(out)   :: gri
     real(wp_), dimension(3,3,params%nray),              intent(out)   :: ggri
+    integer(kind=gray_error),                           intent(inout) :: error
 
     ! local variables
     real(wp_), dimension(3,params%nrayth,params%nrayr) :: xco, du1o
@@ -1089,7 +1091,35 @@ contains
       dxv1 = xc(:,k ,j) -  xc(:,k ,jm)
       dxv2 = xc(:,kp,j) -  xc(:,km,j)
       dxv3 = xc(:,k ,j) - xco(:,k ,j)
-      call solg0(dxv1,dxv2,dxv3,dgu)
+
+      call solg0(dxv1, dxv2, dxv3, dgu)
+
+      ! Check for numerical instability in the computation of ∇S_I.
+      ! This can happen when the distances between the rays and the
+      ! integration step size differ wildly; in such cases the linear
+      ! system Δx⋅∇u = Δu becomes bad conditioned.
+      !
+      ! Source: https://discourse.julialang.org/t/what-is-the-best-way-to-solve-a-ill-conditioned-linear-problem/95894/10
+      block
+        use const_and_precisions, only : comp_eps
+        use utils, only : singvals
+        real(wp_) :: A(3,3), x(3), b(3), sigma(3)
+        real(wp_) :: R, normA, back_error, forw_error
+        A = reshape([dxv1, dxv2, dxv3], shape=[3,3], order=[2,1])
+        sigma = singvals(A)
+        b = [1, 0, 0]
+        x = dgu
+        R = maxval(sigma) / minval(sigma) ! condition number
+        normA = maxval(sigma)             ! L₂ norm of A
+        back_error = norm2(matmul(A, x) - b) / (normA * norm2(b))
+        forw_error = back_error * R * norm2(x)
+        ! TODO: figure out how to tie this threshold
+        ! to the derivatives of the dispersion relation
+        if (forw_error / comp_eps > 5000) then
+          error = raise_error(error, unstable_beam)
+        end if
+      end block
+
       du1(:,k,j) = dgu
       gri(:,jk) = dgu(:)*dffiu
     end do
@@ -1654,7 +1684,7 @@ contains
     use gray_plasma,          only : abstract_plasma
     use dispersion,           only : harmnumber, warmdisp
     use eccd,                 only : setcdcoeff, eccdeff, fjch0, fjch, fjncl
-    use gray_errors,          only : negative_absorption, raise_error
+    use gray_errors,          only : gray_error, negative_absorption, raise_error
     use magsurf_data,         only : fluxval
 
 
@@ -1691,7 +1721,7 @@ contains
     ! lowest/highest resonant harmonic numbers
     integer,   intent(out) :: nhmin, nhmax
     ! ECCD/overall error codes
-    integer,   intent(out) :: iokhawa, error
+    integer(kind=gray_error), intent(inout) :: iokhawa, error
 
     ! local variables
     real(wp_) :: rbavi, rrii, rhop
@@ -1708,7 +1738,6 @@ contains
     nhmin   = 0
     nhmax   = 0
     iokhawa = 0
-    error   = 0
 
     ! Absorption computation
 

src/gray_errors.f90

@@ -5,7 +5,7 @@
 ! function and combined simply with the intrinsic `ior` function.
 module gray_errors
 
-  use logger, only : log_error
+  use logger, only : log_warning
   use, intrinsic :: iso_fortran_env, only : int32
 
   implicit none
@@ -27,16 +27,13 @@ module gray_errors
 
   ! All GRAY errors
 
-  type(error_spec), parameter :: large_npl =        &
+  type(error_spec), parameter :: unstable_beam =  &
     error_spec(offset=0, subcases=2,              &
                mod='gray_core', proc='gray_main', &
-               msg=reshape([character(64) ::        &
+               msg=reshape(['beamtracing may be unstable'], [10], ['']))
-                 'N∥ is too large (>0.99)',         &
-                 'N∥ is too large (>1.05)'          &
-               ], [10], ['']))
 
   type(error_spec), parameter :: dielectric_tensor =            &
-    error_spec(offset=after(large_npl), subcases=2,             &
+    error_spec(offset=after(unstable_beam), subcases=2,         &
                mod='gray_core', proc='gray_main',               &
                msg=reshape([character(64) ::                    &
                  'ε tensor, overflow in `fsup`',                &
@@ -79,7 +76,7 @@ module gray_errors
                  ], [10], ['']))
 
   ! Errors occuring during raytracing
-  type(error_spec), parameter :: raytracing_errors(*) = [large_npl]
+  type(error_spec), parameter :: raytracing_errors(*) = [unstable_beam]
 
   ! Errors occuring during absorption and current drive computations
   type(error_spec), parameter :: ecrh_cd_errors(*) = &
@@ -115,6 +112,19 @@ contains
   end function has_error
 
 
+  pure function has_new_errors(error0, error1)
+    ! Checks whether the `error1` bitmask contains errors
+    ! not already present in `error0`
+
+    ! function arguments
+    integer(kind=gray_error), intent(in) :: error0, error1
+    logical                              :: has_new_errors
+
+    ! check if any bit flipped from 0→1
+    has_new_errors = iand(not(error0), error1) /= 0
+  end function has_new_errors
+
+
   pure function raise_error(error, spec, subcase)
     ! Raise the bits of error `spec` (with optional `subcase` number)
     ! in the `error` bitmask.
@@ -130,17 +140,14 @@ contains
   end function raise_error
 
 
-  subroutine print_err_raytracing(error, step, Npl)
+  subroutine print_err_raytracing(error, step, ray)
     ! Pretty prints raytracing errors
     !
-    ! The error and some context (integration step, N∥)
+    ! The error and some context (integration step, ray)
-    ! is logged to the stderr using the logger module.
+    ! is logged using the logger module.
-
-    use const_and_precisions, only : wp_
 
     ! subroutines arguments
-    integer,   intent(in) :: error, step
+    integer, intent(in) :: error, step, ray
-    real(wp_), intent(in) :: Npl
 
     ! local variables
     integer :: slice        ! a slice of the bitmask
@@ -150,7 +157,7 @@ contains
 
     ! format specifier of the log line
     character(*), parameter :: fmt = &
-      '(a,": ","error=",g0," N∥=",g0.3," step=",g0)'
+      '(a,": ","error=",g0," step=",g0," ray=",g0)'
 
     ! iterate on the known errors
     do i = 1, size(raytracing_errors)
@@ -161,23 +168,23 @@ contains
       ! iterate on the subcases
       do j = 0, spec%subcases - 1
         if (ibits(slice, j, 1) == 0) cycle
-        write(line, fmt) trim(spec%msg(j)), slice * 2**j, Npl, step
+        write(line, fmt) trim(spec%msg(j)), slice * 2**j, step, ray
-        call log_error(line, mod=spec%mod, proc=spec%proc)
+        call log_warning(line, mod=spec%mod, proc=spec%proc)
       end do
     end do
   end subroutine print_err_raytracing
 
 
-  subroutine print_err_ecrh_cd(error, step, Npr, alpha)
+  subroutine print_err_ecrh_cd(error, step, ray, Npr, alpha)
     ! Pretty prints ECRH & CD errors
     !
-    ! The error and some context (integration step, N⊥, α)
+    ! The error and some context (integration step, ray, N⊥, α)
-    ! is logged to the stderr using the logger module.
+    ! is logged using the logger module.
 
     use const_and_precisions, only : wp_
 
     ! subroutines arguments
-    integer,       intent(in) :: error, step
+    integer,       intent(in) :: error, step, ray
     complex(wp_),  intent(in) :: Npr
     real(wp_),     intent(in) :: alpha
 
@@ -189,7 +196,7 @@ contains
 
     ! format specifier of the log line
     character(*), parameter :: fmt = &
-      '(a,": ", "error=",g0, " N⊥=",f0.0,sp,f0.0,"i", " α=",g0, " step=",g0)'
+      '(a,": ", "error=",g0, " N⊥=",f0.0,sp,f0.0,"i", " α=",g0, " step=",g0," ray=",g0)'
 
     ! iterate on the known errors
     do i = 1, size(ecrh_cd_errors)
@@ -201,8 +208,8 @@ contains
       do j = 0, spec%subcases - 1
         if (ibits(slice, j, 1) == 0) cycle
         write(line, fmt) trim(spec%msg(j)), slice * 2**j, &
-                         Npr%re, Npr%im, alpha, step
+                         Npr%re, Npr%im, alpha, step, ray
-        call log_error(line, mod=spec%mod, proc=spec%proc)
+        call log_warning(line, mod=spec%mod, proc=spec%proc)
       end do
     end do
   end subroutine print_err_ecrh_cd

src/main.f90

@@ -6,6 +6,7 @@ program main
   use gray_cli,             only : cli_options, parse_cli_options, &
                                    parse_param_overrides
   use gray_core,            only : gray_main
+  use gray_errors,          only : gray_error
   use gray_equil,           only : abstract_equil, load_equil
   use gray_plasma,          only : abstract_plasma, load_plasma
   use gray_params,          only : gray_parameters, gray_results, &
@@ -13,6 +14,8 @@ program main
                                    make_gray_header
   implicit none
 
+  integer :: status ! exit status
+
   no_save: block
 
   ! CLI options
@@ -24,9 +27,10 @@ program main
   class(abstract_plasma), allocatable :: plasma     !
   type(contour)                       :: limiter    !
   type(gray_results)                  :: results    ! Outputs
-  integer                             :: err     ! Exit code
+  integer(kind=gray_error)            :: sim_error  !
 
   ! Store the original working directory
+  integer :: err
   character(len=256) :: cwd
   call getcwd(cwd)
 
@@ -62,17 +66,17 @@ program main
   ! Initialise antenna parameters before parsing the
   ! cli arguments, so antenna%pos can be overriden
   call init_antenna(params%antenna, err)
-  if (err /= 0) call exit(err)
+  if (err /= 0) call exit(1)
 
   ! Apply CLI overrides to the parameters
   call parse_param_overrides(params)
 
   ! Do some checks on the inputs
   associate (p => params%raytracing)
-    if (p%nrayr < 5) then
+    if (p%igrad == 1 .and. p%nrayr < 5) then
       p%igrad = 0
       call log_message(level=WARNING, mod='main', &
-                       msg='nrayr < 5 ⇒ optical case only')
+                       msg='igrad = 1 but nrayr < 5: disabling beamtracing')
     end if
     if (p%nrayr == 1) p%nrayth = 1
 
@@ -83,11 +87,11 @@ program main
 
   ! Read and initialise the equilibrium and limiter objects
   call load_equil(params%equilibrium, equil, limiter, err)
-  if (err /= 0) call exit(err)
+  if (err /= 0) call exit(1)
 
   ! Read and initialise the plasma state object
   call load_plasma(params, equil, plasma, err)
-  if (err /= 0) call exit(err)
+  if (err /= 0) call exit(1)
 
   ! Create a simple limiter, if necessary
   if (.not. allocated(limiter%R) .or. params%raytracing%ipass < 0) then
@@ -102,12 +106,14 @@ program main
   if (allocated(opts%sum_filelist)) then
     call log_message(level=INFO, mod='main', msg='summing profiles')
     call sum_profiles(params, equil, opts%sum_filelist, opts%output_dir, results)
+    status = 0
   else
     ! Activate the given output tables
     params%misc%active_tables = opts%tables
 
     ! Finally run the simulation
-    call gray_main(params, equil, plasma, limiter, results, err)
+    call gray_main(params, equil, plasma, limiter, results, sim_error)
+    status = 2 * sim_error
   end if
 
   print_res: block
@@ -138,6 +144,7 @@ program main
       if (err /= 0) then
         call log_message('failed to save '//trim(tbl%title)//' table', &
                          level=ERROR, mod='main')
+        status = 1
       end if
       end associate
     end do
@@ -145,6 +152,9 @@ program main
 
   end block no_save
 
+  ! Exit status is >1 if any error occurred
+  call exit(status)
+
 contains
 
   subroutine init_antenna(params, err)

src/utils.f90

@@ -218,4 +218,76 @@ contains
                  sin(phi),  cos(phi)], shape=[2,2], order=[2,1])
   end function rotate
 
+
+  pure function det(M)
+    ! Computes the determinant of a real 3x3 matrix M
+
+    ! function arguments
+    real(wp_), intent(in) :: M(3,3)
+    real(wp_)             :: det
+
+    associate(a => M(1,1), b => M(1,2), c => M(1,3), &
+              d => M(2,1), e => M(2,2), f => M(2,3), &
+              g => M(3,1), h => M(3,2), i => M(3,3))
+      det = (a*e*i + b*f*g + c*d*h) - (a*f*h + b*d*i + c*e*g)
+    end associate
+  end function det
+
+
+  pure function eigenvals(A)
+    ! Computes the eigenvalues of a real symmetric 3x3 matrix A
+    !
+    ! Source: https://doi.org/10.1145%2F355578.366316
+    use const_and_precisions, only : pi, one
+
+    ! function arguments
+    real(wp_), intent(in) :: A(3,3)
+    real(wp_)             :: eigenvals(3)
+
+    ! local variables
+    real(wp_) :: B(3,3)
+    real(wp_) :: p, q, r, p1, p2
+    real(wp_) :: phi
+    integer   :: i
+
+    p1 = A(1,2)**2 + A(1,3)**2 + A(2,3)**2
+
+    if (p1 == 0) then
+      ! A is diagonal
+      eigenvals = [A(1,1), A(2,2), A(3,3)]
+      return
+    end if
+
+    q  = sum([(A(i,i), i=1,3)])/3  ! trace(A)
+    p2 = (A(1,1) - q)**2 + (A(2,2) - q)**2 + (A(3,3) - q)**2 + 2*p1
+    p  = sqrt(p2 / 6)
+    B  = (A - q * diag([one, one, one])) / p
+    r = det(B) / 2
+
+    ! Ensure r ∈ [-1,1]
+    if (r <= -1) then
+       phi = pi / 3
+    else if (r >= 1) then
+       phi = 0
+    else
+       phi = acos(r) / 3
+    end if
+
+    eigenvals = q - 2*p*[-cos(phi), cos(phi + pi/3), sin(phi + pi/6)]
+  end function eigenvals
+
+
+  pure function singvals(A)
+    ! Computes the singular values of a real 3x3 matrix A
+
+    ! function arguments
+    real(wp_), intent(in) :: A(3,3)
+    real(wp_)             :: singvals(3)
+
+    ! By definition, the singular values of A are the
+    ! eigenvalues value of A†A. These are necessarily
+    ! real by the spectral theorem.
+    singvals = sqrt(eigenvals(matmul(A, transpose(A))))
+  end function singvals
+
 end module utils

tests/__init__.py

@@ -48,8 +48,9 @@ class GrayTest:
         # run gray to generate the candidate outputs
         proc = run_gray(cls.inputs, cls.candidate, params=cls.gray_params,
                         binary=options.binary)
-        assert proc.returncode == 0, \
+
-               f"gray failed with exit code {proc.returncode}"
+        # 0: all good, 1: input errors, >1: simulation errors
+        assert proc.returncode != 1, 'gray failed with exit code 1'
 
         # store the stderr for manual inspection
         with open(str(cls.candidate / 'log'), 'w') as log:
@@ -256,8 +257,11 @@ def run_gray(inputs: Path, outputs: Path,
     ] + list(itertools.chain(*params)) + options
     proc = subprocess.run(args, capture_output=True, text=True)
 
+    print()
     if proc.returncode != 0:
         # show the log on errors
+        print(f'Errors occurred (exit status {proc.returncode}), showing log:')
+        print(*proc.args)
         print(proc.stderr)
         print(proc.stdout)
     return proc