TokaMaker: Add error handling to .save_eqdsk() (failed tracing)

 - Update error reporting in `.get_q()` to report flux coordinate where tracing failed

src/physics/grad_shaf.F90

@@ -4665,6 +4665,7 @@ subroutine gs_get_qprof(gseq,nr,psi_q,prof,dl,rbounds,zbounds,ravgs)
 real(8), parameter :: tol=1.d-10
 integer(4) :: i,j,cell
 type(gsinv_interp), target :: field
+CHARACTER(LEN=80) :: error_str
 !---
 raxis=gseq%o_point(1)
 zaxis=gseq%o_point(2)
@@ -4737,8 +4738,8 @@ subroutine gs_get_qprof(gseq,nr,psi_q,prof,dl,rbounds,zbounds,ravgs)
   pt_last=pt
   !---Skip point if trace fails
   if(active_tracer%status/=1)THEN
-    WRITE(*,*)j,pt
-    CALL oft_warn("gs_get_qprof: Trace did not complete")
+    WRITE(error_str,"(A,F10.4)")"gs_get_qprof: Trace did not complete at psi = ",1.d0-psi_q(j)
+    CALL oft_warn(error_str)
     CYCLE
   end if
   IF(j==1)THEN

src/physics/grad_shaf_util.F90

@@ -845,10 +845,11 @@ END SUBROUTINE gs_analyze
 !---------------------------------------------------------------------------
 !> Needs docs
 !---------------------------------------------------------------------------
-subroutine gs_save_decon(gseq,npsi,ntheta)
+subroutine gs_save_decon(gseq,npsi,ntheta,error_str)
 class(gs_eq), intent(inout) :: gseq
 integer(4), intent(in) :: npsi
 integer(4), intent(in) :: ntheta
+CHARACTER(LEN=80), OPTIONAL, INTENT(out) :: error_str
 type(gsinv_interp), target :: field
 type(oft_lag_brinterp) :: psi_int
 real(8) :: gop(3,3),psi_surf(1),pt_last(3)
@@ -859,6 +860,7 @@ subroutine gs_save_decon(gseq,npsi,ntheta)
 integer(4) :: j,k,cell,io_unit
 TYPE(spline_type) :: rz
 !---
+IF(PRESENT(error_str))error_str=""
 WRITE(*,'(2A)')oft_indent,'Saving DCON file'
 CALL oft_increase_indent
 !---
@@ -911,6 +913,9 @@ subroutine gs_save_decon(gseq,npsi,ntheta)
 ALLOCATE(ptout(3,active_tracer%maxsteps+1))
 !$omp do schedule(dynamic,1)
 do j=1,npsi-1
+  IF(PRESENT(error_str))THEN
+    IF(error_str/="")CYCLE
+  END IF
   !---------------------------------------------------------------------------
   ! Trace contour
   !---------------------------------------------------------------------------
@@ -928,7 +933,16 @@ subroutine gs_save_decon(gseq,npsi,ntheta)
   CALL tracinginv_fs(pt,ptout)
   pt_last=pt
   !---Exit if trace fails
-  if(active_tracer%status/=1)call oft_abort('Trace did not complete.','gs_save_decon',__FILE__)
+  IF(active_tracer%status/=1)THEN
+    IF(PRESENT(error_str))THEN
+      !$omp critical
+      WRITE(error_str,"(A,F10.4)")"Tracing failed at psi = ",1.d0-psi_surf
+      !$omp end critical
+      CYCLE
+    ELSE
+      call oft_abort('Trace did not complete.','gs_save_decon',__FILE__)
+    END IF
+  END IF
   !---------------------------------------------------------------------------
   ! Perform cubic spline interpolation
   !---------------------------------------------------------------------------
@@ -964,6 +978,12 @@ subroutine gs_save_decon(gseq,npsi,ntheta)
 CALL active_tracer%delete
 DEALLOCATE(ptout)
 !$omp end parallel
+IF(PRESENT(error_str))THEN
+  IF(error_str/="")THEN
+    DEALLOCATE(cout,rout,zout)
+    RETURN
+  END IF
+END IF
 !---Information for O-point
 rout(npsi,:)=raxis
 zout(npsi,:)=zaxis
@@ -1021,7 +1041,7 @@ end subroutine gs_save_decon
 !---------------------------------------------------------------------------
 !> Save equilibrium to General Atomics gEQDSK file
 !---------------------------------------------------------------------------
-subroutine gs_save_eqdsk(gseq,filename,nr,nz,rbounds,zbounds,run_info,limiter_file,psi_pad)
+subroutine gs_save_eqdsk(gseq,filename,nr,nz,rbounds,zbounds,run_info,limiter_file,psi_pad,error_str)
 class(gs_eq), intent(inout) :: gseq !< Equilibrium to save
 CHARACTER(LEN=OFT_PATH_SLEN), intent(in) :: filename 
 integer(4), intent(in) :: nr !< Number of radial points for flux/psi grid
@@ -1031,6 +1051,7 @@ subroutine gs_save_eqdsk(gseq,filename,nr,nz,rbounds,zbounds,run_info,limiter_fi
 CHARACTER(LEN=36), intent(in) :: run_info !< Run information string [36]
 CHARACTER(LEN=OFT_PATH_SLEN), intent(in) :: limiter_file !< Path to limiter file
 REAL(8), intent(in) :: psi_pad !< Padding at LCFS in normalized units
+CHARACTER(LEN=80), OPTIONAL, INTENT(out) :: error_str
 !
 real(8) :: psi_surf,rmax,x1,x2,raxis,zaxis,xr
 real(8) :: pt(3),pt_last(3),f(3),psi_tmp(1),gop(3,3)
@@ -1049,6 +1070,7 @@ subroutine gs_save_eqdsk(gseq,filename,nr,nz,rbounds,zbounds,run_info,limiter_fi
 REAL(8), ALLOCATABLE, DIMENSION(:) :: fpol,pres,ffprim,pprime,qpsi
 REAL(8), ALLOCATABLE, DIMENSION(:,:) :: psirz
 !---
+IF(PRESENT(error_str))error_str=""
 WRITE(*,'(2A)')oft_indent,'Saving EQDSK file'
 CALL oft_increase_indent
 !---
@@ -1104,6 +1126,9 @@ subroutine gs_save_eqdsk(gseq,filename,nr,nz,rbounds,zbounds,run_info,limiter_fi
 active_tracer%inv=.TRUE.
 !$omp do schedule(dynamic,1)
 do j=1,nr
+  IF(PRESENT(error_str))THEN
+    IF(error_str/="")CYCLE
+  END IF
   !---------------------------------------------------------------------------
   ! Trace contour
   !---------------------------------------------------------------------------
@@ -1127,7 +1152,16 @@ subroutine gs_save_eqdsk(gseq,filename,nr,nz,rbounds,zbounds,run_info,limiter_fi
     END IF
     pt_last=pt
     !---Exit if trace fails
-    if(active_tracer%status/=1)call oft_abort('Trace did not complete.','gs_save_eqdsk',__FILE__)
+    IF(active_tracer%status/=1)THEN
+      IF(PRESENT(error_str))THEN
+        !$omp critical
+        WRITE(error_str,"(A,F10.4)")"Tracing failed at psi = ",1.d0-(psi_surf-x1)/(x2-x1)
+        !$omp end critical
+        CYCLE
+      ELSE
+        call oft_abort('Trace did not complete.','gs_save_eqdsk',__FILE__)
+      END IF
+    END IF
   END IF
   IF(j==nr)THEN
     !---------------------------------------------------------------------------
@@ -1169,6 +1203,13 @@ subroutine gs_save_eqdsk(gseq,filename,nr,nz,rbounds,zbounds,run_info,limiter_fi
 end do
 CALL active_tracer%delete
 !$omp end parallel
+IF(PRESENT(error_str))THEN
+  IF(error_str/="")THEN
+    DEALLOCATE(rout,zout)
+    DEALLOCATE(fpol,pres,ffprim,pprime,qpsi,psirz)
+    RETURN
+  END IF
+END IF
 !---Extrapolate q on axis
 f(1) = x2
 f(2) = x2 - (x2-x1)*(1.d0/REAL(nr-1,8))
@@ -1249,7 +1290,7 @@ subroutine gs_save_eqdsk(gseq,filename,nr,nz,rbounds,zbounds,run_info,limiter_fi
 END IF
 CALL oft_decrease_indent
 !---
-DEALLOCATE(rout,zout)
+DEALLOCATE(rout,zout,rlim,zlim)
 DEALLOCATE(fpol,pres,ffprim,pprime,qpsi,psirz)
 end subroutine gs_save_eqdsk
 !---------------------------------------------------------------------------

src/python/OpenFUSIONToolkit/TokaMaker/_core.py

@@ -1273,7 +1273,10 @@ def save_eqdsk(self,filename,nr=65,nz=65,rbounds=None,zbounds=None,run_info='',l
             zbounds = numpy.r_[self.lim_contour[:,1].min(), self.lim_contour[:,1].max()]
             dr = zbounds[1]-zbounds[0]
             zbounds += numpy.r_[-1.0,1.0]*dr*0.05
-        tokamaker_save_eqdsk(cfilename,c_int(nr),c_int(nz),rbounds,zbounds,crun_info,c_double(lcfs_pad))
+        cstring = c_char_p(b""*200)
+        tokamaker_save_eqdsk(cfilename,c_int(nr),c_int(nz),rbounds,zbounds,crun_info,c_double(lcfs_pad),cstring)
+        if cstring.value != b'':
+            raise Exception(cstring.value)
 
     def eig_wall(self,neigs=4,pm=False):
         '''! Compute eigenvalues (1 / Tau_L/R) for conducting structures

src/python/OpenFUSIONToolkit/TokaMaker/_interface.py

@@ -199,7 +199,7 @@ class tokamaker_settings_struct(c_struct):
     [c_int_ptr,c_double_ptr_ptr])
 
 tokamaker_save_eqdsk = ctypes_subroutine(oftpy_lib.tokamaker_save_eqdsk, #(filename,nr,nz,rbounds,zbounds,run_info,psi_pad)
-    [c_char_p, c_int, c_int, ctypes_numpy_array(numpy.float64,1), ctypes_numpy_array(numpy.float64,1), c_char_p, c_double])
+    [c_char_p, c_int, c_int, ctypes_numpy_array(numpy.float64,1), ctypes_numpy_array(numpy.float64,1), c_char_p, c_double, c_char_p])
 ## @endcond
 
 

src/python/wrappers/tokamaker_f.F90

@@ -920,19 +920,21 @@ END SUBROUTINE tokamaker_set_coil_vsc
 !------------------------------------------------------------------------------
 !> Needs docs
 !------------------------------------------------------------------------------
-SUBROUTINE tokamaker_save_eqdsk(filename,nr,nz,rbounds,zbounds,run_info,psi_pad) BIND(C,NAME="tokamaker_save_eqdsk")
+SUBROUTINE tokamaker_save_eqdsk(filename,nr,nz,rbounds,zbounds,run_info,psi_pad,error_str) BIND(C,NAME="tokamaker_save_eqdsk")
 CHARACTER(KIND=c_char), INTENT(in) :: filename(80) !< Needs docs
 CHARACTER(KIND=c_char), INTENT(in) :: run_info(36) !< Needs docs
 INTEGER(c_int), VALUE, INTENT(in) :: nr !< Needs docs
 INTEGER(c_int), VALUE, INTENT(in) :: nz !< Needs docs
 REAL(c_double), INTENT(in) :: rbounds(2) !< Needs docs
 REAL(c_double), INTENT(in) :: zbounds(2) !< Needs docs
 REAL(c_double), VALUE, INTENT(in) :: psi_pad !< Needs docs
+CHARACTER(KIND=c_char), INTENT(out) :: error_str(80) !< Needs docs
 CHARACTER(LEN=36) :: run_info_f
-CHARACTER(LEN=80) :: filename_tmp,lim_file
+CHARACTER(LEN=80) :: filename_tmp,lim_file,error_flag
 CALL copy_string_rev(run_info,run_info_f)
 CALL copy_string_rev(filename,filename_tmp)
 lim_file='none'
-CALL gs_save_eqdsk(gs_global,filename_tmp,nr,nz,rbounds,zbounds,run_info_f,lim_file,psi_pad)
+CALL gs_save_eqdsk(gs_global,filename_tmp,nr,nz,rbounds,zbounds,run_info_f,lim_file,psi_pad,error_flag)
+CALL copy_string(TRIM(error_flag),error_str)
 END SUBROUTINE tokamaker_save_eqdsk
 END MODULE tokamaker_f