scripts/gray/visual.py: degrade instead of failing on missing files

The poloidal flux matrix is written in fortran order!

default.nix

@@ -33,7 +33,7 @@ let
     src = builtins.filterSource sieve ./.;
 
     format = "other";
-    propagatedBuildInputs = with pkgs.python3Packages; [ numpy matplotlib ];
+    propagatedBuildInputs = with pkgs.python3Packages; [ numpy scipy matplotlib ];
     makeFlags = [ "PREFIX=$(out)" ];
 
     dontConfigure = true;

scripts/gray/__init__.py

@@ -58,7 +58,7 @@ def read_eqdsk(filepath: Path) -> Eqdsk:
     skip(file, nr * 4)
 
     # read flux and safety factor
-    psi = read(file, nr * nz).reshape(nr, nz)
+    psi = read(file, nr * nz).reshape((nr, nz), order='F').T
     q = read(file, nr)
 
     # normalise flux

scripts/gray/visual.py

@@ -39,6 +39,20 @@ def cli_args() -> argparse.Namespace:
     return parser.parse_args()
 
 
+def placeholder(ax: plt.Axes, text: str):
+    '''
+    Puts a placeholder text at the center of an axis
+    '''
+    left, width = .25, .5
+    bottom, height = .25, .5
+    right = left + width
+    top = bottom + height
+
+    ax.text(0.5*(left+right), 0.5*(top+bottom), text,
+            horizontalalignment='center',
+            verticalalignment='center')
+
+
 def align_yaxis(*axes: [plt.Axes]):
     '''
     Aligns the origins of two axes
@@ -116,30 +130,33 @@ def plot_poloidal(inputs: Path, outputs: Path, ax: plt.Axes):
     ax.set_ylabel('$z$ / m')
 
     # load flux surfaces
-    surfaces = gray.read_table(outputs / 'flux-surfaces.71.txt')
-    surfaces = surfaces.reshape(-1, int(surfaces['i'].max()))
+    try:
+        surfaces = gray.read_table(outputs / 'flux-surfaces.71.txt')
+        surfaces = surfaces.reshape(-1, int(surfaces['i'].max()))
 
-    # plot plasma boundary
-    bound = np.argmin(1 - surfaces['ψ_n'][:, 0])
-    ax.plot(surfaces[bound]['R'], surfaces[bound]['z'],
-            c='xkcd:slate grey', label='plasma boundary')
+        # plot plasma boundary
+        bound = np.argmin(1 - surfaces['ψ_n'][:, 0])
+        ax.plot(surfaces[bound]['R'], surfaces[bound]['z'],
+                c='xkcd:slate grey', label='plasma boundary')
 
-    # manual contourplot of ψ(R,z)
-    for i, surf in enumerate(surfaces[:bound]):
-        ax.plot(surf['R'], surf['z'], c='xkcd:grey',
-                label='flux surfaces' if i == 0 else '')
+        # manual contourplot of ψ(R,z)
+        for i, surf in enumerate(surfaces[:bound]):
+            ax.plot(surf['R'], surf['z'], c='xkcd:grey',
+                    label='flux surfaces' if i == 0 else '')
 
-    # label the rational surfaces
-    labels = ['q=2', 'q=3/2', 'q=1']
-    for label, surf in zip(labels, surfaces[1+bound:][::-1]):
-        line = np.vstack([surf['R'], surf['z']]).T
-        cont = ContourSet(ax, [1], [[line]], filled=False,
-                          colors='xkcd:ocean blue',
-                          linestyles='--')
-        ax.clabel(cont, [1], inline=True,
-                  fontsize=10, fmt={1: label})
-    ax.plot(np.nan, np.nan, '--', color='xkcd:ocean blue',
-            label='rational surfaces')
+        # label the rational surfaces
+        labels = ['q=2', 'q=3/2', 'q=1']
+        for label, surf in zip(labels, surfaces[1+bound:][::-1]):
+            line = np.vstack([surf['R'], surf['z']]).T
+            cont = ContourSet(ax, [1], [[line]], filled=False,
+                              colors='xkcd:ocean blue',
+                              linestyles='--')
+            ax.clabel(cont, [1], inline=True,
+                      fontsize=10, fmt={1: label})
+        ax.plot(np.nan, np.nan, '--', color='xkcd:ocean blue',
+                label='rational surfaces')
+    except FileNotFoundError:
+        pass
 
     # load limiter
     limiter = gray.get_limiter(conf, inputs)
@@ -186,12 +203,15 @@ def plot_toroidal(inputs: Path, outputs: Path, ax: plt.Axes):
     limiter = gray.get_limiter(conf, inputs)
 
     # plot plasma boundary
-    surfaces = gray.read_table(outputs / 'flux-surfaces.71.txt')
-    boundary = surfaces[np.isclose(surfaces['ψ_n'], 1, 1e-3)]
+    try:
+        surfaces = gray.read_table(outputs / 'flux-surfaces.71.txt')
+        boundary = surfaces[np.isclose(surfaces['ψ_n'], 1, 1e-3)]
 
-    # plot plasma boundary
-    plot_circle(ax, radius=boundary['R'].min(), color='xkcd:slate gray')
-    plot_circle(ax, radius=boundary['R'].max(), color='xkcd:slate gray')
+        # plot plasma boundary
+        plot_circle(ax, radius=boundary['R'].min(), color='xkcd:slate gray')
+        plot_circle(ax, radius=boundary['R'].max(), color='xkcd:slate gray')
+    except FileNotFoundError:
+        pass
 
     # plot limiter
     if limiter[0].size > 0:
@@ -374,36 +394,42 @@ def plot_inputs(inputs: Path, outputs: Path, axes: [plt.Axes]):
     '''
     from matplotlib.colors import TwoSlopeNorm
 
-    maps = gray.read_table(outputs / 'inputs-maps.72.txt')
+    conf = gray.read_conf(inputs / 'gray.ini')
 
-    # filter valid points
-    maps = maps[maps['ψ_n'] > 0]
-    flux = maps['R'], maps['z'], maps['ψ_n']
+    try:
+        maps = gray.read_table(outputs / 'inputs-maps.72.txt')
+        # filter valid points
+        maps = maps[maps['ψ_n'] > 0]
+        flux = maps['R'], maps['z'], maps['ψ_n']
+    except FileNotFoundError:
+        maps = None
+        placeholder(axes['B'], 'inputs-maps.72.txt not found')
 
     # contour levels
+    ψ_max = maps['ψ_n'].max() if maps is not None else 1.2
     inner_levels = np.linspace(0, 0.99, 8)
-    outer_levels = np.linspace(0.9991, maps['ψ_n'].max()*0.99, 6)
+    outer_levels = np.linspace(0.9991, ψ_max*0.99, 6)
     all_levels = np.concatenate([inner_levels, outer_levels])
 
     # contour style
-    norm = TwoSlopeNorm(vmin=0, vcenter=1, vmax=maps['ψ_n'].max())
+    norm = TwoSlopeNorm(vmin=0, vcenter=1, vmax=ψ_max)
     cmap = plt.cm.RdYlGn_r
     borders = dict(colors='xkcd:grey', linestyles='-', linewidths=0.8)
 
     # interpolated equilibrium
-    interp = axes['B']
-    interp.set_title('poloidal flux', loc='right')
-    interp.set_xlabel('$R$ / m')
-    interp.set_ylabel('$z$ / m')
-    interp.tricontourf(*flux, levels=all_levels, norm=norm, cmap=cmap)
-    interp.tricontour(*flux, levels=all_levels, **borders)
-    interp.tricontour(*flux, levels=[1], colors='xkcd:slate gray')
-    interp.plot(np.nan, np.nan, c='xkcd:slate gray', label='plasma boundary')
+    if maps is not None:
+        interp = axes['B']
+        interp.set_title('poloidal flux', loc='right')
+        interp.set_xlabel('$R$ / m')
+        interp.set_ylabel('$z$ / m')
+        interp.tricontourf(*flux, levels=all_levels, norm=norm, cmap=cmap)
+        interp.tricontour(*flux, levels=all_levels, **borders)
+        interp.tricontour(*flux, levels=[1], colors='xkcd:slate gray')
+        interp.plot(np.nan, np.nan, c='xkcd:slate gray', label='plasma boundary')
 
-    # add limiter
-    conf = gray.read_conf(inputs / 'gray.ini')
-    limiter = gray.get_limiter(conf, inputs)
-    interp.plot(*limiter, c='xkcd:black', label='limiter')
+        # add limiter
+        limiter = gray.get_limiter(conf, inputs)
+        interp.plot(*limiter, c='xkcd:black', label='limiter')
 
     # original MHD equilibrium
     orig = axes['A']
@@ -418,14 +444,18 @@ def plot_inputs(inputs: Path, outputs: Path, axes: [plt.Axes]):
         orig.plot(*eqdsk.limiter, c='xkcd:black')
         orig.plot(*eqdsk.boundary, c='xkcd:slate gray')
     else:
-        orig.axis('off')
+        placeholder(orig, 'not an EQDSK')
 
     # colorbar
     bar = plt.colorbar(plt.cm.ScalarMappable(norm, cmap), cax=axes['C'])
     bar.set_label(label='normalised poloidal flux', labelpad=10)
 
     # Plasma radial profiles
-    profiles = gray.read_table(outputs / 'kinetic-profiles.55.txt')
+    try:
+        profiles = gray.read_table(outputs / 'kinetic-profiles.55.txt')
+    except FileNotFoundError:
+        placeholder(axes['D'], 'kinetic-profiles.55.txt not found')
+        return
 
     axes['D'].set_title('plasma profiles', loc='right')
     add_alt_axis(axes['D'], profiles['ρ_t'], profiles['ψ_n']**0.5,

src/gray_params.f90

@@ -135,14 +135,14 @@ module gray_params
 
   ! Raytracing parameters
   type raytracing_parameters
-    real(wp_)                :: dst                   ! Integration step size
-    real(wp_)                :: rwmax = 1             ! Normalized maximum radius of beam power
-    integer                  :: nray, nrayr, nrayth   ! Total/radial/angular number of rays
-    logical                  :: igrad = .false.       ! Complex eikonal switch
-    integer                  :: nstep = 12000         ! Max number of integration steps
-    integer(kind(step_enum)) :: idst  = STEP_TIME     ! Choice of the integration variable
-    integer                  :: ipass = 1             ! Number of plasma passes
-    logical                  :: ipol  = .false.       ! Whether to compute wave polarisation (from chi, psi)
+    real(wp_)                :: dst                  ! Integration step size
+    real(wp_)                :: rwmax = 1            ! Normalized maximum radius of beam power
+    integer                  :: nray, nrayr, nrayth  ! Total/radial/angular number of rays
+    integer                  :: igrad = 0            ! Complex eikonal switch
+    integer                  :: nstep = 12000        ! Max number of integration steps
+    integer(kind(step_enum)) :: idst  = STEP_TIME    ! Choice of the integration variable
+    integer                  :: ipass = 1            ! Number of plasma passes
+    logical                  :: ipol  = .false.      ! Whether to compute wave polarisation (from chi, psi)
   end type
 
   ! EC resonant heating & Current Drive parameters
@@ -341,7 +341,7 @@ contains
     call append(header, line)
 
     ! code parameters
-    write(line, '("# COD igrad idst ipass ipol:",l4,2(1x,i4),1x,l4)') &
+    write(line, '("# COD igrad idst ipass ipol:",3(1x,i4),1x,l4)') &
       params%raytracing%igrad, params%raytracing%idst, &
       params%raytracing%ipass, params%raytracing%ipol
     call append(header, line)

src/gray_plasma.f90

@@ -21,7 +21,7 @@ module gray_plasma
   end type
 
   abstract interface
-    pure subroutine density_sub(self, psin, dens, ddens)
+    subroutine density_sub(self, psin, dens, ddens)
       ! Computes the density its first derivative as a function of
       ! normalised poloidal flux.
       !
@@ -32,7 +32,7 @@ module gray_plasma
       real(wp_),              intent(out) :: dens, ddens  ! density and first derivative
     end subroutine density_sub
 
-    pure function temp_fun(self, psin) result(temp)
+    function temp_fun(self, psin) result(temp)
       ! Computes the temperature as a function of the
       ! normalised poloidal flux.
       !
@@ -43,7 +43,7 @@ module gray_plasma
       real(wp_)                          :: temp
     end function temp_fun
 
-    pure function zeff_fun(self, psin) result(zeff)
+    function zeff_fun(self, psin) result(zeff)
       ! Computes the effective charge Z_eff as a
       ! function of the normalised poloidal flux.
       import :: abstract_plasma, wp_
@@ -97,7 +97,7 @@ module gray_plasma
 
 contains
 
-  pure subroutine analytic_density(self, psin, dens, ddens)
+  subroutine analytic_density(self, psin, dens, ddens)
     ! subroutine arguments
     class(analytic_plasma), intent(in)  :: self
     real(wp_),              intent(in)  :: psin         ! normalised poloidal flux
@@ -116,7 +116,7 @@ contains
   end subroutine analytic_density
 
 
-  pure function analytic_temp(self, psin) result(temp)
+  function analytic_temp(self, psin) result(temp)
     ! function arguments
     class(analytic_plasma), intent(in) :: self
     real(wp_),              intent(in) :: psin
@@ -139,7 +139,7 @@ contains
   end function analytic_temp
 
 
-  pure function analytic_zeff(self, psin) result(zeff)
+  function analytic_zeff(self, psin) result(zeff)
     ! function arguments
     class(analytic_plasma), intent(in) :: self
     real(wp_),              intent(in) :: psin
@@ -155,7 +155,7 @@ contains
   end function analytic_zeff
 
 
-  pure subroutine numeric_density(self, psin, dens, ddens)
+  subroutine numeric_density(self, psin, dens, ddens)
     use logger, only : log_error
 
     ! subroutine arguments
@@ -221,16 +221,16 @@ contains
       end block
     end if
 
-    ! if (dens < 0) then
-    !   write (msg, '("negative density:", 2(x,a,"=",g0.3))') &
-    !     'ne', dens, 'ψ', psin
-    !   call log_error(msg, mod='coreprofiles', proc='density')
-    ! end if
+    if (dens < 0) then
+      write (msg, '("negative density:", 2(x,a,"=",g0.3))') &
+        'ne', dens, 'ψ', psin
+      call log_error(msg, mod='coreprofiles', proc='density')
+    end if
 
   end subroutine numeric_density
 
 
-  pure function numeric_temp(self, psin) result(temp)
+  function numeric_temp(self, psin) result(temp)
     ! function arguments
     class(numeric_plasma), intent(in) :: self
     real(wp_),             intent(in) :: psin
@@ -246,7 +246,7 @@ contains
   end function numeric_temp
 
 
-  pure function numeric_zeff(self, psin) result(zeff)
+  function numeric_zeff(self, psin) result(zeff)
     ! function arguments
     class(numeric_plasma), intent(in) :: self
     real(wp_),             intent(in) :: psin

src/gray_tables.f90

@@ -121,7 +121,7 @@ contains
     real(wp_) :: rho_t, J_phi, dens, ddens
     real(wp_) :: psi_n, rho_p, drho_p
 
-    call tbl%init(title='kinetic-profiles', id=55, active=is_active(params, 25), &
+    call tbl%init(title='kinetic-profiles', id=55, active=is_active(params, 55), &
                   labels=[character(64) :: 'ψ_n', 'ρ_t', 'n_e', 'T_e', 'q', 'J_φ'])
 
     if (.not. tbl%active) return

src/main.f90

@@ -77,10 +77,10 @@ program main
 
   ! Do some checks on the inputs
   associate (p => params%raytracing)
-    if (p%igrad .and. p%nrayr < 5) then
-      p%igrad = .false.
+    if (p%igrad == 1 .and. p%nrayr < 5) then
+      p%igrad = 0
       call log_message(level=WARNING, mod='main', &
-                       msg='igrad = .true. but nrayr < 5: disabling beamtracing')
+                       msg='igrad = 1 but nrayr < 5: disabling beamtracing')
     end if
     if (p%nrayr == 1) p%nrayth = 1
 

src/main_convert.f90

@@ -156,8 +156,8 @@ contains
     write(u, fmt) "nrayth", params%raytracing%nrayth
     if (params%raytracing%rwmax /= defaults%raytracing%rwmax) &
       write(u, fmt) "rwmax", params%raytracing%rwmax
-    if (params%raytracing%igrad .neqv. defaults%raytracing%igrad) &
-      write(u, fmt) "igrad", format_bool(params%raytracing%igrad)
+    if (params%raytracing%igrad /= defaults%raytracing%igrad) &
+      write(u, fmt) "igrad", params%raytracing%igrad
     if (params%raytracing%ipass /= defaults%raytracing%ipass) &
       write(u, fmt) "ipass", params%raytracing%ipass
     if (params%raytracing%ipol .neqv. defaults%raytracing%ipol) &

tests/06-ITER-startup/inputs/gray.ini

@@ -1,7 +1,7 @@
 [raytracing]
 nrayr = 11
 nrayth = 16
-igrad = true
+igrad = 1
 dst = 0.1
 nstep = 8000