From f02cbe981ead5cc578276c1d517b11d57fb595e7 Mon Sep 17 00:00:00 2001
From: Antoine <antoine.hoffmann@epfl.ch>
Date: Sun, 14 Apr 2024 13:09:16 +0200
Subject: [PATCH] New tutorial to run the code for the first time - new_prob.sh
 creates a small environement with everything needed to run the first
 simulation - parameter files and minimal python analyzing scripts are now
 present in basic_script

---
 .gitignore                                    |   2 +
 bash scripts/fort_00.90                       | 107 -----------------
 bash scripts/submit_00.cmd                    |  13 ---
 basic_scripts/parameters_CBC.90               |  86 ++++++++++++++
 basic_scripts/parameters_Zpinch.90            |  90 +++++++++++++++
 basic_scripts/python_utilities/fourier.py     |  31 +++++
 basic_scripts/python_utilities/load_data.py   |  50 ++++++++
 .../python_utilities/minimal_analysis.py      |  82 +++++++++++++
 basic_scripts/python_utilities/tools.py       |  23 ++++
 basic_scripts/submit_marconi_example.cmd      |  15 +++
 basic_scripts/tutorial.md                     |  82 +++++++++++++
 fort.90                                       | 109 ------------------
 new_prob.sh                                   |  81 +++++++++++++
 13 files changed, 542 insertions(+), 229 deletions(-)
 delete mode 100644 bash scripts/fort_00.90
 delete mode 100644 bash scripts/submit_00.cmd
 create mode 100644 basic_scripts/parameters_CBC.90
 create mode 100644 basic_scripts/parameters_Zpinch.90
 create mode 100644 basic_scripts/python_utilities/fourier.py
 create mode 100644 basic_scripts/python_utilities/load_data.py
 create mode 100644 basic_scripts/python_utilities/minimal_analysis.py
 create mode 100644 basic_scripts/python_utilities/tools.py
 create mode 100644 basic_scripts/submit_marconi_example.cmd
 create mode 100644 basic_scripts/tutorial.md
 delete mode 100644 fort.90
 create mode 100644 new_prob.sh

diff --git a/.gitignore b/.gitignore
index 44012ae4..9a3873b1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -49,3 +49,5 @@ gyacomo_dbg
 gyacomo23*
 wk/parameters/profiles
 local/dirs.inc
+simulations
+__pycache__
\ No newline at end of file
diff --git a/bash scripts/fort_00.90 b/bash scripts/fort_00.90
deleted file mode 100644
index 183029f9..00000000
--- a/bash scripts/fort_00.90	
+++ /dev/null
@@ -1,107 +0,0 @@
-&BASIC
-  nrun       = 100000000
-  dt         = 0.002
-  tmax       = 100
-  maxruntime = 356400
-  job2load   = -1
-/
-&GRID
-  pmax  = 2
-  jmax  = 1
-  Nx    = 192
-  Lx    = 90
-  Ny    = 48
-  Ly    = 200
-  Nz    = 32
-  SG    = .false.
-  Nexc  = 0
-/
-&GEOMETRY
-  geom     = 'miller'
-  q0       =-2.15
-  shear    = 3.62
-  eps      = 0.28
-  kappa    = 1.53
-  s_kappa  = 0.77
-  delta    = 0.23
-  s_delta  = 1.05
-  zeta     =-0.01
-  s_zeta   =-0.17
-  parallel_bc = 'dirichlet'
-  shift_y = 0
-  Npol    = 1
-  PB_PHASE= .false.
-/
-&OUTPUT_PAR
-  dtsave_0d = 0.5
-  dtsave_1d = -1
-  dtsave_2d = -1
-  dtsave_3d = 5.0
-  dtsave_5d = 10
-  write_doubleprecision = .false.
-  write_gamma = .true.
-  write_hf    = .true.
-  write_phi   = .true.
-  write_Na00  = .true.
-  write_Napj  = .true.
-  write_dens  = .true.
-  write_temp  = .true.
-/
-&MODEL_PAR
-LINEARITY = 'nonlinear'
-RM_LD_T_EQ= .false.
-  Na      = 2
-  mu_x    = 1.0
-  mu_y    = 1.0
-  N_HD    = 4
-  mu_z    = 2.0
-  HYP_V   = 'hypcoll'
-  mu_p    = 0
-  mu_j    = 0
-  nu      = 0.5
-  k_gB    = 1
-  k_cB    = 1
-  lambdaD = 0
-  beta    = 0.0034
-  ADIAB_E = .false.
-  ADIAB_I = .false.
-  MHD_PD  = .true.
-/
-&CLOSURE_PAR
-  hierarchy_closure='truncation'
-  dmax             =-1
-  nonlinear_closure='truncation'
-  nmax             =0
-/
-&SPECIES
-  name_  = 'ions' 
-  tau_   = 1
-  sigma_ = 1
-  q_     = 1
-  K_N_   = 1.33
-  K_T_   = 8.25
-/
-&SPECIES
-  name_  = 'electrons' 
-  tau_   = 1
-  sigma_ = 0.04!0.023338
-  q_     = -1
-  K_N_   = 1.33
-  K_T_   = 12
-/
-&COLLISION_PAR
-  collision_model = 'DG'
-  GK_CO      = .true.
-  INTERSPECIES    = .true.
-  mat_file        = '/Users/ahoffmann/gyacomo/iCa/null'
-  collision_kcut  = 1
-/
-&INITIAL_CON
-  INIT_OPT      = 'blob'
-  init_background  = 0
-  init_noiselvl = 1e-05
-  iseed         = 42
-/
-&TIME_INTEGRATION_PAR
-  numerical_scheme = 'RK4'
-/
\ No newline at end of file
diff --git a/bash scripts/submit_00.cmd b/bash scripts/submit_00.cmd
deleted file mode 100644
index 0d1f9436..00000000
--- a/bash scripts/submit_00.cmd	
+++ /dev/null
@@ -1,13 +0,0 @@
-#!/bin/bash
-#SBATCH --job-name=test_pipe
-#SBATCH --time=00:02:00
-#SBATCH --nodes=1
-#SBATCH --cpus-per-task=1
-#SBATCH --ntasks-per-node=48
-#SBATCH --mem=64GB
-#SBATCH --error=err_00.txt
-#SBATCH --output=out_00.txt
-#SBATCH --account=FUA36_TSVVT422
-#SBATCH --partition=skl_fua_dbg
-#SBATCH --qos=normal
-srun --cpu-bind=cores ./gyacomo 2 24 1 0
\ No newline at end of file
diff --git a/basic_scripts/parameters_CBC.90 b/basic_scripts/parameters_CBC.90
new file mode 100644
index 00000000..288b95aa
--- /dev/null
+++ b/basic_scripts/parameters_CBC.90
@@ -0,0 +1,86 @@
+&BASIC
+  nrun       = 99999999 ! number of time step to perform
+  dt         = 0.01     ! time step (not adaptive)
+  tmax       = 50       ! maximal time [c_s/R]
+  maxruntime = 72000    ! maximal wallclock runtime [sec]
+  job2load   = -1       ! index of the previous run to restart (-1:start a new run)
+/
+&GRID
+  pmax   = 2            ! maximal degree of the Hermite basis (parallel velocity)
+  jmax   = 1            ! maximal degree of the Laguerre basis (magnetic moment)
+  Nx     = 64           ! number of points in the radial direction
+  Lx     = 120          ! size of the box in the radial direction [rho_s]
+  Ny     = 48           ! bumber of points in the binormal direction
+  Ly     = 120          ! size of the box in the binormal direction [rho_s]
+  Nz     = 16           ! number of points in the magnetic field direction
+  SG     = .t.          ! use a staggered grid in z (experimental)
+  Nexc   = 0           ! to fullfill the sheared boundary condition (set -1 for automatic)
+/
+&GEOMETRY
+  geom   = 's-alpha'    ! geometry model (Z-pinch,s-alpha,miller,circular)
+  q0     = 1.4          ! safety factor
+  shear  = 0.8          ! magnetic shear
+  eps    = 0.18         ! inverse aspect ratio
+  kappa  = 1.0          ! elongation
+  s_kappa= 0.0          ! elongation derivative
+  delta  = 0.0          ! triangularity
+  s_delta= 0.0          ! triangularity derivative
+  zeta   = 0.0          ! squareness
+  s_zeta = 0.0          ! squareness derivative
+  parallel_bc = 'dirichlet' ! to change the type of parallel boundary condition (experimental)
+  shift_y= 0.0          ! to add a shift in the parallel BC (experimental)
+  Npol   = 1.0          ! set the length of the z domain (-pi N_pol < z < pi N_pol)
+/
+&DIAGNOSTICS
+  dtsave_0d = 1         ! period of 0D diagnostics (time traces)
+  dtsave_1d = -1        ! period of 1D diagnostics (nothing)
+  dtsave_2d = -1        ! period of 2D diagnostics (nothing)
+  dtsave_3d = 2         ! period of 3D diagnostics (phi, Aparallel, fluid moments etc.)
+  dtsave_5d = 20        ! period of 5D diagnostics (full set of GMs)
+/
+&MODEL
+  LINEARITY = 'nonlinear' ! set if we solve the linear or nonlinear problem
+  Na      = 1           ! number of species (this sets the number of species namelists to be read)
+  mu_x    = 0.0         ! numerical diffusion parameter in the radial direction
+  mu_y    = 0.0         ! numerical diffusion parameter in the binormal direction
+  N_HD    = 4           ! degree of the numerical diffusion
+  mu_z    = 0.0         ! numerical diffusion parameter in the parallel direction (fourth order only)
+  HYP_V   = 'hypcoll'   ! numerical diffusion scheme in the velocity space (experimental)
+  mu_p    = 0.0         ! numerical diffusion parameter in the parallel velocity (experimental)
+  mu_j    = 0.0         ! numerical diffusion parameter in the magnetic moment (experimental)
+  nu      = 0.05        ! collision rate, ~0.5 GENE parameter (better to use it instead of num. diff. in velocities)
+  beta    = 0.0         ! plasma beta (not in percent)
+  ADIAB_E = .t.         ! Use an adiabatic electron model (required if Na = 1)
+/
+&CLOSURE
+  hierarchy_closure='truncation' ! closure scheme
+  dmax = -1             ! set the maximal degree of moment to evolve (-1 evolves all th)
+  nonlinear_closure='anti_laguerre_aliasing' ! set the truncation of the Laguerre convolution (truncation,full_sum,anti_laguerre_aliasing)
+  nmax = -1             ! set the maximal degree in the truncation of the Laguerre convolution (experimental)
+/
+&SPECIES ! Defines the species a (out of Na species)
+ name_ = 'ions' ! name of the species a
+ tau_  = 1.0    ! temperature (Ta/Te)
+ sigma_= 1.0    ! sqrt mass ratio (sqrt(ma/mi))
+ q_    = 1.0    ! charge (qa/e)
+ k_N_  = 2.2    ! density background gradient, grad ln N (omn in GENE)
+ k_T_  = 7.0    ! temperature background gradient, grad ln T (omT in GENE)
+/
+&SPECIES ! defines electrons (if Na=1, not read)
+ name_ = 'electrons'
+ tau_  = 1.0
+ sigma_= 0.023338
+ q_    =-1.0
+ k_N_  = 2.0
+ k_T_  = 0.4
+/
+&COLLISION
+  collision_model = 'DG'   !LB/DG/SG/PA/LD (Lenhard-Bernstein, Dougherty, Sugama, pitch angle, Landau)
+  GK_CO           = .t.    ! gyrokinetic version of the collision operator (or longwavelength only)
+/
+&INITIAL
+  INIT_OPT         = 'blob' ! initilization of the system ('phi' put noise in electrostatic, 'blob' is like 'ppj' in GENE)
+/
+&TIME_INTEGRATION
+  numerical_scheme = 'RK4' ! numerical scheme for time integration (RK2,3,4 etc.)
+/
diff --git a/basic_scripts/parameters_Zpinch.90 b/basic_scripts/parameters_Zpinch.90
new file mode 100644
index 00000000..46501b32
--- /dev/null
+++ b/basic_scripts/parameters_Zpinch.90
@@ -0,0 +1,90 @@
+&BASIC
+  nrun       = 99999999 ! number of time step to perform
+  dt         = 0.01     ! time step (not adaptive)
+  tmax       = 100       ! maximal time [c_s/R]
+  maxruntime = 72000    ! maximal wallclock runtime [sec]
+  job2load   = -1       ! index of the previous run to restart (-1:start a new run)
+/
+&GRID
+  pmax   = 2            ! maximal degree of the Hermite basis (parallel velocity)
+  jmax   = 1            ! maximal degree of the Laguerre basis (magnetic moment)
+  Nx     = 64           ! number of points in the radial direction
+  Lx     = 200          ! size of the box in the radial direction [rho_s]
+  Ny     = 48           ! bumber of points in the binormal direction
+  Ly     = 60           ! size of the box in the binormal direction [rho_s]
+  Nz     = 1            ! number of points in the magnetic field direction
+  SG     = .f.          ! use a staggered grid in z (experimental)
+  Nexc   = -1           ! to fullfill the sheared boundary condition (set -1 for automatic)
+/
+&GEOMETRY
+  geom   = 'Z-pinch'    ! geometry model (Z-pinch,s-alpha,miller,circular)
+  q0     = 0.0          ! safety factor
+  shear  = 0.0          ! magnetic shear
+  eps    = 0.0          ! inverse aspect ratio
+  kappa  = 1.0          ! elongation
+  s_kappa= 0.0          ! elongation derivative
+  delta  = 0.0          ! triangularity
+  s_delta= 0.0          ! triangularity derivative
+  zeta   = 0.0          ! squareness
+  s_zeta = 0.0          ! squareness derivative
+  parallel_bc = 'dirichlet' ! to change the type of parallel boundary condition (experimental)
+  shift_y= 0.0          ! to add a shift in the parallel BC (experimental)
+  Npol   = 1.0          ! set the length of the z domain (-pi N_pol < z < pi N_pol)
+/
+&DIAGNOSTICS
+  dtsave_0d = 1         ! period of 0D diagnostics (time traces)
+  dtsave_1d = -1        ! period of 1D diagnostics (nothing)
+  dtsave_2d = -1        ! period of 2D diagnostics (nothing)
+  dtsave_3d = 1         ! period of 3D diagnostics (phi, Aparallel, fluid moments etc.)
+  dtsave_5d = 10        ! period of 5D diagnostics (full set of GMs)
+/
+&MODEL
+  LINEARITY = 'nonlinear' ! set if we solve the linear or nonlinear problem
+  Na      = 2           ! number of species (this sets the number of species namelists to be read)
+  mu_x    = 1.0         ! numerical diffusion parameter in the radial direction
+  mu_y    = 1.0         ! numerical diffusion parameter in the binormal direction
+  N_HD    = 4           ! degree of the numerical diffusion
+  mu_z    = 0.0         ! numerical diffusion parameter in the parallel direction (fourth order only)
+  HYP_V   = 'hypcoll'   ! numerical diffusion scheme in the velocity space (experimental)
+  mu_p    = 0.0         ! numerical diffusion parameter in the parallel velocity (experimental)
+  mu_j    = 0.0         ! numerical diffusion parameter in the magnetic moment (experimental)
+  nu      = 0.1         ! collision rate, ~0.5 GENE parameter (better to use it instead of num. diff. in velocities)
+  beta    = 0.0         ! plasma beta (not in percent)
+  ADIAB_E = .f.         ! Use an adiabatic electron model (required if Na = 1)
+  ADIAB_I = .f.         ! Use an adiabatic ion model (experimental)
+/
+&CLOSURE
+  hierarchy_closure='truncation' ! closure scheme
+  dmax = -1             ! set the maximal degree of moment to evolve (-1 evolves all th)
+  nonlinear_closure='anti_laguerre_aliasing' ! set the truncation of the Laguerre convolution (truncation,full_sum,anti_laguerre_aliasing)
+  nmax = -1             ! set the maximal degree in the truncation of the Laguerre convolution (experimental)
+/
+&SPECIES ! Defines the species a (out of Na species)
+ name_ = 'ions' ! name of the species a
+ tau_  = 1.0    ! temperature (Ta/Te)
+ sigma_= 1.0    ! sqrt mass ratio (sqrt(ma/mi))
+ q_    = 1.0    ! charge (qa/e)
+ k_N_  = 2.0    ! density background gradient, grad ln N (omn in GENE)
+ k_T_  = 0.4    ! temperature background gradient, grad ln T (omT in GENE)
+/
+&SPECIES ! defines electrons (if not, use ADIAB_E = true)
+ ! electrons
+ name_ = 'electrons'
+ tau_  = 1.0
+ sigma_= 0.023338
+ q_    =-1.0
+ k_N_  = 2.0
+ k_T_  = 0.4
+/
+&COLLISION
+  collision_model = 'DG'   !LB/DG/SG/PA/LD (Lenhard-Bernstein, Dougherty, Sugama, pitch angle, Landau)
+  GK_CO           = .t.    ! gyrokinetic version of the collision operator (or longwavelength only)
+/
+&INITIAL
+  INIT_OPT         = 'phi' ! initilization of the system ('phi' put noise in electrostatic, 'blob' is like 'ppj' in GENE)
+  init_noiselvl    = 0.005 ! amplitude of the noise init (for 'phi' init only)
+  iseed            = 42    ! random seed
+/
+&TIME_INTEGRATION
+  numerical_scheme = 'RK4' ! numerical scheme for time integration (RK2,3,4 etc.)
+/
diff --git a/basic_scripts/python_utilities/fourier.py b/basic_scripts/python_utilities/fourier.py
new file mode 100644
index 00000000..c44e51af
--- /dev/null
+++ b/basic_scripts/python_utilities/fourier.py
@@ -0,0 +1,31 @@
+import numpy as np
+def kx_to_x(var_kx, nx0, axis=-3):
+    """ Perform inverse FFT on kx spectral direction of variable
+
+    Note: The fft and ifft in python and GENE/IDL have the factor 1/N switched!
+    :param var_kx: Variable in kx space
+    :param nx0: Number of kx grid points
+    :param axis: Which axis of var_kx is the x direction, by default the third last one
+    :returns: variable in real x space
+    """
+    var_x = np.fft.ifft(nx0*var_kx, axis=axis)
+    var_x = np.real_if_close(var_x, tol=1e5)
+    return var_x
+
+
+def ky_to_y(var_ky, nky0, axis=-2):
+    """ Perform inverse FFT on ky spectral direction of variable
+
+    The GENE data only include the non-negative ky components, so we need to use the real
+    valued FFT routines
+    Note: The fft and ifft in python and GENE/IDL have the factor 1/N switched!
+    :param var_ky: Variable in kx space
+    :param nky0: Number of ky grid points
+    :param axis: Which axis of var_kx is the x direction, by default the third last one
+    :returns: variable in real y space
+    """
+    # The GENE data only include the non-negative ky components, so we need to use the real
+    # valued FFT routines
+    var_y = np.fft.irfft(2*nky0*var_ky, n=2*nky0, axis=axis)
+    var_y = np.real_if_close(var_y, tol=1e5)
+    return var_y
diff --git a/basic_scripts/python_utilities/load_data.py b/basic_scripts/python_utilities/load_data.py
new file mode 100644
index 00000000..07e1be6f
--- /dev/null
+++ b/basic_scripts/python_utilities/load_data.py
@@ -0,0 +1,50 @@
+import h5py
+import numpy as np
+import tools
+def load_data_0D(filename,dname):
+    with h5py.File(filename, 'r') as file:
+        # Load time data
+        time  = file['data/var0d/time']
+        time  = time[:]
+        var0D = file['data/var0d/'+dname]
+        var0D = var0D[:]
+    return time, var0D
+
+def load_data_3D_frame(filename,dname,tframe):
+    with h5py.File(filename, 'r') as file:
+        # load time
+        time  = file['data/var3d/time']
+        time  = time[:]
+        # find frame
+        iframe = tools.closest_index(time,tframe)
+        tf     = time[iframe]
+        # Load data
+        if dname == 'Ni00':
+            data = file[f'data/var3d/Na00/{iframe:06d}']
+            data = data[:,:,:,0]
+        else:
+            data = file[f'data/var3d/{dname}/{iframe:06d}']
+        data = data[:]
+        data = data['real']+1j*data['imaginary'] 
+        # Load the grids
+        kxgrid = file[f'data/grid/coordkx']
+        kygrid = file[f'data/grid/coordky']
+        zgrid  = file[f'data/grid/coordz']
+    return time,data,tf
+
+def load_grids(filename):
+    with h5py.File(filename, 'r') as file:
+        # Load the grids
+        kxgrid = file[f'data/grid/coordkx']
+        kygrid = file[f'data/grid/coordky']
+        zgrid  = file[f'data/grid/coordz']
+        pgrid  = file[f'data/grid/coordp']
+        jgrid  = file[f'data/grid/coordj']
+        Nx     = kxgrid.size
+        Nky    = kygrid.size
+        Ny     = 2*(Nky-1)
+        Lx     = 2*np.pi/kxgrid[1]
+        Ly     = 2*np.pi/kygrid[1]
+        xgrid  = np.linspace(-Lx/2,Lx/2,Nx)
+        ygrid  = np.linspace(-Ly/2,Ly/2,Ny)
+    return xgrid,kxgrid,ygrid,kygrid,zgrid,pgrid,jgrid
\ No newline at end of file
diff --git a/basic_scripts/python_utilities/minimal_analysis.py b/basic_scripts/python_utilities/minimal_analysis.py
new file mode 100644
index 00000000..5d361108
--- /dev/null
+++ b/basic_scripts/python_utilities/minimal_analysis.py
@@ -0,0 +1,82 @@
+"""
+This script loads data from the HDF5 output file of GYACOMO ('outputs_00.h5').
+It plots the time traces of the heat and particle fluxes in the radial direction by loading
+the following datasets:
+- 'data/var0d/time': Time values
+- 'data/var0d/hflux_x': Heat flux data
+- 'data/var0d/pflux_x': Particle flux data
+"""
+import h5py
+import load_data as loader
+import matplotlib.pyplot as plt
+import numpy as np
+import sys
+import tools
+
+# Open the HDF5 file
+if __name__ == "__main__":
+    # No command-line argument provided, use default filename
+   filename = "outputs_00.h5"
+elif len(sys.argv) == 2:
+    # Command-line argument provided
+    arg = sys.argv[1]
+    if arg.isdigit():
+        # Argument is an integer, format it as XX
+        filename = f"outputs_{int(arg):02d}.h5"
+    else:
+        # Argument is a string
+        filename = arg
+else:
+    print("Usage: python minimal_analysis.py [filename, int or string (opt)]")
+    sys.exit(1)
+
+t0D, hflux_x = loader.load_data_0D(filename,'hflux_x')
+t0D, pflux_x = loader.load_data_0D(filename,'pflux_x')
+t3D, phi, tf = loader.load_data_3D_frame(filename,'phi',10000)
+t3D, Ni00,tf = loader.load_data_3D_frame(filename,'Ni00',10000)
+x,kx,y,ky,z,p,j  = loader.load_grids(filename)
+
+phi = tools.zkxky_to_xy_const_z(phi,-1)
+Ni00= tools.zkxky_to_xy_const_z(Ni00,-1)
+
+# Plot hflux_x and pflux_x against time
+fig, axes = plt.subplots(2, 2, figsize=(10, 6))
+# Plot hflux_x
+nt0D = t0D.size
+if hflux_x.size > nt0D :
+    axes[0,0].plot(t0D, hflux_x[::2], label='ion heat flux')
+    axes[0,0].plot(t0D, hflux_x[1::2], label='electron heat flux')
+else:
+    axes[0,0].plot(t0D, hflux_x, label='ion heat flux')
+axes[0,0].set_title('radial heat flux')
+axes[0,0].set_xlabel('t c_s/R')
+axes[0,0].set_ylabel('Q_x')
+
+# Plot pflux_x
+if pflux_x.size > nt0D :
+    axes[1,0].plot(t0D, pflux_x[::2], label='ion heat flux')
+    axes[1,0].plot(t0D, pflux_x[1::2], label='electron heat flux')
+else:
+    axes[0,0].plot(t0D, pflux_x, label='ion heat flux')
+axes[1,0].set_title('radial particle flux')
+axes[1,0].set_xlabel('t c_s/R')
+axes[1,0].set_ylabel('P_x')
+
+# Plot slice of electrostatic potential at constant z
+axes[0,1].imshow(Ni00, extent=[x[0], x[-1], y[0], y[-1]],cmap='viridis',interpolation='quadric')
+#axes[0,1].pcolor(x,y,Ni00, cmap='viridis')
+axes[0,1].set_title(f'Ni00 (z = 0, t = {tf:3.1f})')
+axes[0,1].set_xlabel('x')
+axes[0,1].set_ylabel('y')
+#fig.colorbar(im, ax=axes[0,1], label='Ni00')
+
+# Plot slice of electrostatic potential at constant z
+axes[1,1].imshow(phi, extent=[x[0], x[-1], y[0], y[-1]],cmap='viridis',interpolation='quadric')
+#axes[1,1].pcolor(x,y,phi, cmap='viridis')
+axes[1,1].set_title(f'\phi (z = 0, t = {tf:3.1f})')
+axes[1,1].set_xlabel('x')
+axes[1,1].set_ylabel('y')
+#fig.colorbar(im, ax=axes[1,1], label='Potential')
+
+plt.tight_layout()
+plt.show()
diff --git a/basic_scripts/python_utilities/tools.py b/basic_scripts/python_utilities/tools.py
new file mode 100644
index 00000000..90c8d5b3
--- /dev/null
+++ b/basic_scripts/python_utilities/tools.py
@@ -0,0 +1,23 @@
+import fourier
+import numpy as np
+
+def zkxky_to_xy_const_z(array, iz):
+    # Get shape of the phi array
+    Nz, Nkx, Nky, = array.shape
+    if iz < 0: #outboard midplane for negative iz
+        iz = Nz // 2  # Using the middle value for z
+
+    array = array[iz,:,:]
+    array = fourier.kx_to_x(array,Nkx,-2)
+    array = fourier.ky_to_y(array,Nky-1,-1)
+    array = np.transpose(array)
+    return array
+    
+def closest_index(array, v):
+    # Compute absolute differences between each element of the array and v
+    absolute_diff = np.abs(array - v)
+    
+    # Find the index of the minimum difference
+    closest_index = np.argmin(absolute_diff)
+    
+    return closest_index
\ No newline at end of file
diff --git a/basic_scripts/submit_marconi_example.cmd b/basic_scripts/submit_marconi_example.cmd
new file mode 100644
index 00000000..054b8eef
--- /dev/null
+++ b/basic_scripts/submit_marconi_example.cmd
@@ -0,0 +1,15 @@
+#!/bin/bash
+#SBATCH --job-name=submit_example   # name of the job
+#SBATCH --time=00:15:00             # walltime
+#SBATCH --nodes=1                   # total number of nodes
+#SBATCH --cpus-per-task=1           # Number of threads per task (=1, no OpenMP feature yet)
+#SBATCH --ntasks-per-node=48        # MPI tasks (marconi, max 48)
+#SBATCH --mem=64GB                  # memory requirement
+#SBATCH --error=err.txt             # std-error file
+#SBATCH --output=out.txt            # std-output file
+#SBATCH --account=my_account        # account name
+#SBATCH --partition=skl_fua_dbg     # partition name
+#SBATCH --qos=normal                # quality of service
+
+ # run the code
+srun --cpu-bind=cores -l ./gyacomo.exe 2 24 1 0
\ No newline at end of file
diff --git a/basic_scripts/tutorial.md b/basic_scripts/tutorial.md
new file mode 100644
index 00000000..25c10d04
--- /dev/null
+++ b/basic_scripts/tutorial.md
@@ -0,0 +1,82 @@
+## Running Gyacomo Tutorial
+This tutorial provides guidance on running Gyacomo in a new installation. It is designed to be used in conjunction with the `new_prob.sh` script, providing a refresher on the basics of Gyacomo commands.
+It is meant to be followed once `sh new_prob.sh` has been run in `/gyacomo` which creates an example problem directory `/gyacomo/simulations/problem_01`. The following commands are expected to be run from the example problem directory.
+
+### How to Run the Code Locally After Executing `new_prob.sh` in `/gyacomo` Directory?
+- For a single-core run, execute:
+    ```bash
+    ./gyacomo.exe
+    ```
+- For a multi-core run, use:
+    ```bash
+    mpirun -np N ./gyacomo.exe np nky nz
+    ```
+    Here, `N` is the total number of processes, `np` is the number in the `p` direction (Hermite polynomials), `nky` is the number of poloidal wavenumbers, and `nz` is the number of parallel planes. The relation `N = np x nky x nz` must be satisfied.
+In both commands above, the input parameters are expected to be in a `fort.90` file present in the directory where the code is run. You can append an additional integer `X` at the end of both single and multi-core runs, which points to an input file `fort_XX.90` where `XX` is the two-digit version of the number `X`. This allows for clearer restarts by indexing each consecutive run.
+
+#### Examples
+- Single-core run with parameters located in `fort_00.90`:
+    ```bash
+    ./gyacomo.exe 0
+    ```
+- Multi-core run with 2 cores in Hermite, 4 in ky, 1 in z, and reading the file `fort_00.90`:
+    ```bash
+    mpirun -np 8 ./gyacomo.exe 2 4 1 0
+    ```
+- Same as above but redirecting the standard terminal output to the file `out_00.txt` and running in the background:
+    ```bash
+    mpirun -np 8 ./gyacomo.exe 2 4 1 0 > out_00.txt &
+    ```
+    You can monitor the simulation progress with:
+    ```bash
+    tail -f out_00.txt
+    ```
+
+### How to Stop the Simulation?
+Here are the stopping conditions of GYACOMO: 
+- A nan (overflow) is detected in any field, this can be caused by a too high `dt` (CFL condition) or an unsufficient resolution (`Nkx`, `Nky`, `Nz`)
+- The number of steps exceeds the number provided in `nrun`.
+- The simulation run time exceeds the time defined in `maxruntime`.
+- The simulation reaches the maximal physical time `tmax`.
+- A file named `mystop`is present in the current directory.
+Thus, you can smoothly stop your current simulation by creating an empty file named `mystop` in the directory where the code runs (here `gyacomo/simulations/problem_01`):
+```bash
+touch mystop
+```
+The code checks if such a file exists every 100 steps (this is set in the `gyacomo/src/tesend.F90` module). If it finds it, it will run the ending procedure and remove the file from the directory. 
+
+### Checking the Output
+Once the simulation is finished, the Python script `minimal_analysis.py` provides a basic example of result analysis. If you followed the tutorial, you should be able to obtain plots by typing:
+```bash
+python minimal_analysis.py
+```
+
+### Restarting a Simulation
+In the parameter file, the parameter `job2load` defines from which previous outputs the code has to continue the simulation. If `job2load=-1`, a new start is made, and the output is located in the `outputs_00.h5` file in the current directory.
+
+If the `outputs_XX` file exists, you can continue the run by setting `job2load=X` (where `XX` is the double-digit version of the integer `X`). It's recommended to link all runs with an input file `fort_XX.90`. In any case, you can find the input file used in the simulations in `outputs_XX.h5`, located in `outputs_XX.h5/files/STDIN.00`.
+
+You can use `minimal_analysis.py X` to analyze the restart. The script will here look for the `outputs_XX.h5` file.
+
+#### Example
+This is a minimal example for a restart procedure
+1. Copy the `fort_00.90` to a new `fort_01.90` input file.
+2. Adapt the parameters of `fort_01.90` (Tmax, grads, anything).
+3. **IMPORTANT:** Ensure that `fort_01.90` has `job2load = 0` (otherwise, it will restart a simulation from 0).
+4. Now you can run:
+    ```bash
+    mpirun -np 8 ./gyacomo.exe 2 4 1 1 > out_01.txt &
+    ```
+    Note the `1` input in the fourth position; the code reads `fort_01.90`. We also direct the std output to a new file out_01.txt
+5. You can monitor the run with:
+    ```bash
+    tail -f out_01.txt
+    ```
+5. Once the simulation is done (end or mystop), you can analyze the new data using:
+    ```bash
+    python minimal_analysis.py 1
+    ```
+
+### Running on Marconi
+
+You can adapt and use the provided template `submit_marconi.cmd`.
\ No newline at end of file
diff --git a/fort.90 b/fort.90
deleted file mode 100644
index 27a45669..00000000
--- a/fort.90
+++ /dev/null
@@ -1,109 +0,0 @@
-&BASIC
-  nrun       = 99999999
-  dt         = 0.01
-  tmax       = 50
-  maxruntime = 72000
-  job2load   = -1
-/
-&GRID
-  pmax   = 2
-  jmax   = 1
-  Nx     = 64
-  Lx     = 200
-  Ny     = 48
-  Ly     = 60
-  Nz     = 1
-  SG     = .f.
-  Nexc   = 1
-/
-&GEOMETRY
-  geom   = 'Z-pinch'
-  q0     = 0.0
-  shear  = 0.0
-  eps    = 0.0
-  kappa  = 1.0
-  s_kappa= 0.0
-  delta  = 0.0
-  s_delta= 0.0
-  zeta   = 0.0
-  s_zeta = 0.0
-  parallel_bc = 'shearless'
-  shift_y= 0.0
-  Npol   = 1.0
-/
-&DIAGNOSTICS
-  dtsave_0d = 1
-  dtsave_1d = -1
-  dtsave_2d = -1
-  dtsave_3d = 1
-  dtsave_5d = 10
-  write_doubleprecision = .f.
-  write_gamma = .t.
-  write_hf    = .t.
-  write_phi   = .t.
-  write_Na00  = .t.
-  write_Napj  = .t.
-  write_dens  = .t.
-  write_fvel  = .t.
-  write_temp  = .t.
-  !diag_mode   = 'txtonly'
-/
-&MODEL
-  LINEARITY = 'nonlinear'
-  Na      = 2 ! number of species
-  mu_x    = 1.0
-  mu_y    = 1.0
-  N_HD    = 4
-  mu_z    = 0.0
-  HYP_V   = 'hypcoll'
-  mu_p    = 0.0
-  mu_j    = 0.0
-  nu      = 0.1
-  beta    = 0.0
-  ADIAB_E = .f.
-  !KN_MODEL = 'taylor'
-  !ORDER = 1
-  !KN_MODEL = 'pade'
-  !ORDER_NUM = 2 
-  !ORDER_DEN = 2
-/
-&CLOSURE
-  hierarchy_closure='truncation'
-  dmax = -1
-  nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing)
-  nmax = 0
-/
-&SPECIES
- ! ions
- name_ = 'ions'
- tau_  = 1.0
- sigma_= 1.0
- q_    = 1.0
- k_N_  = 2.0
- k_T_  = 0.4
-/
-&SPECIES
- ! electrons
- name_ = 'electrons'
- tau_  = 1.0
- sigma_= 0.023338
- q_    =-1.0
- k_N_  = 2.0
- k_T_  = 0.4
-/
-&COLLISION
-  collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
-  GK_CO           = .t.
-  INTERSPECIES    = .true.
-  mat_file       = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
-/
-&INITIAL
-  INIT_OPT         = 'phi' !(phi,blob)
-  ACT_ON_MODES     = 'donothing'
-  init_background  = 0.0
-  init_noiselvl    = 0.005
-  iseed            = 42
-/
-&TIME_INTEGRATION
-  numerical_scheme = 'RK4'
-/
diff --git a/new_prob.sh b/new_prob.sh
new file mode 100644
index 00000000..91075ede
--- /dev/null
+++ b/new_prob.sh
@@ -0,0 +1,81 @@
+#!/bin/bash
+gyacomo_dir=$(pwd)
+echo "Setup of a new problem in the simulations directory..."
+
+# Check if an argument is provided
+if [ $# -ne 1 ]; then
+    echo "- using CBC parameter"
+    echo "  (you can ask for a Zpinch base case by adding "ZBC" after the call of the script)"
+    inputfile=$gyacomo_dir/basic_scripts/parameters_CBC.90
+else
+    # Check the value of the argument
+    if [ "$1" = "CBC" ]; then
+        echo "- using CBC parameter"
+        inputfile=$gyacomo_dir/basic_scripts/parameters_CBC.90
+    elif [ "$1" = "ZBC" ]; then
+        echo "- using Zpinch parameter"
+        inputfile=$gyacomo_dir/basic_scripts/parameters_Zpinch.90
+    else
+        echo "- using CBC parameter"
+        echo "(you can ask for a Zpinch test by adding "ZBC" after the call of the script)"
+        inputfile=$gyacomo_dir/basic_scripts/parameters_CBC.90
+    fi
+fi
+
+mkdir -p $gyacomo_dir/simulations
+cd $gyacomo_dir/simulations
+
+# Function to get the next available problem folder name
+get_next_problem_folder() {
+    local max_num=0
+    for folder in problem_*; do
+        num="${folder#problem_}"
+        if [[ $num =~ ^[0-99]+$ ]] && [ "$num" -gt "$max_num" ]; then
+            max_num="$num"
+        fi
+    done
+    printf "problem_%02d" "$((max_num + 1))"
+}
+
+# Function to display usage information
+usage() {
+    echo "- usage: $0 [base_case_type] [folder_name]"
+    echo "If folder_name is not provided, default is the next available 'problem_xx'"
+    exit 1
+}
+
+# Check if an argument is provided
+#if [ "$#" -gt 1 ]; then
+#    usage
+#fi
+
+# Set the folder name based on the input or default
+folder_name="${2:-$(get_next_problem_folder)}"
+
+# Check if the executable exists
+if [ ! -x $gyacomo_dir/bin/gyacomo23_dp ]; then
+    echo "Error: Executable 'bin/gyacomo23_dp' not found or not executable. Please compile the code first."
+    exit 1
+fi
+
+# Create the folder with the specified name if it doesn't exist
+mkdir -p "$folder_name"
+
+# Copy basic parameter file
+cp $inputfile "$folder_name/fort.90"
+
+# Copy the marconi job submission example
+cp $gyacomo_dir/basic_scripts/submit_marconi_example.cmd "$folder_name/submit_marconi.cmd"
+
+# Copy tutorial file
+cp $gyacomo_dir/basic_scripts/tutorial.md "$folder_name/."
+
+# Create a symbolic link to the executable
+ln -s $gyacomo_dir/bin/gyacomo23_dp "$folder_name/gyacomo.exe"
+
+# Create a symbolic link to the basic python analysis script
+ln -s $gyacomo_dir/basic_scripts/python_utilities/minimal_analysis.py "$folder_name/."
+
+echo "- folder 'simulations/$folder_name' created with symbolic link to the executable."
+echo "- check the tutorial file in there!"
+echo "...done"
\ No newline at end of file
-- 
GitLab