diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index e744f98923eb61341c9d0de39856a6ecc139f9e9..dcecd1eb9a0d89a4677d2fae410e79ce3364d219 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -62,7 +62,7 @@ while(CONTINUE)
% Loading from output file
% CPUTIME = h5readatt(filename,'/data/input','cpu_time');
% DT_SIM = h5readatt(filename,'/data/input','dt');
- [Pe, Je, Pi, Ji, kx, ky, z] = load_grid_data(filename);
+ [P, J, kx, ky, z] = load_grid_data(filename);
W_GAMMA = strcmp(h5readatt(filename,'/data/input','write_gamma'),'y');
W_HF = strcmp(h5readatt(filename,'/data/input','write_hf' ),'y');
W_PHI = strcmp(h5readatt(filename,'/data/input','write_phi' ),'y');
@@ -78,8 +78,8 @@ while(CONTINUE)
BETA = 0;
end
% Check polynomials degrees
- Pe_new= numel(Pe); Je_new= numel(Je);
- Pi_new= numel(Pi); Ji_new= numel(Ji);
+ Pe_new= numel(P); Je_new= numel(J);
+ Pi_new= numel(P); Ji_new= numel(J);
if(Pe_max < Pe_new); Pe_max = Pe_new; end;
if(Je_max < Je_new); Je_max = Je_new; end;
if(Pi_max < Pi_new); Pi_max = Pi_new; end;
@@ -305,13 +305,13 @@ else
DATA.PSI = PSI_;
DATA.KIN_E=KIN_E;
% grids
- DATA.Pe = Pe; DATA.Pi = Pi;
- DATA.Je = Je; DATA.Ji = Ji;
+ DATA.Pe = P; DATA.Pi = P;
+ DATA.Je = J; DATA.Ji = J;
DATA.kx = kx; DATA.ky = ky; DATA.z = z; DATA.Npol = -z(1)/pi;
DATA.x = x; DATA.y = y;
DATA.ikx0 = ikx0; DATA.iky0 = iky0;
DATA.Nx = Nx; DATA.Ny = Ny; DATA.Nz = Nz; DATA.Nkx = Nkx; DATA.Nky = Nky;
- DATA.Pmaxe = numel(Pe); DATA.Pmaxi = numel(Pi); DATA.Jmaxe = numel(Je); DATA.Jmaxi = numel(Ji);
+ DATA.Pmaxe = numel(P); DATA.Pmaxi = numel(P); DATA.Jmaxe = numel(J); DATA.Jmaxi = numel(J);
DATA.dir = DIRECTORY;
DATA.localdir = DIRECTORY;
DATA.param_title=['$\nu_{',DATA.CONAME,'}=$', num2str(DATA.NU), ...
diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index d449552ed0bc09b73e36b63f6241864e318e4abc..b97ab56be1247f12890edd7c71f7e0c6068f9fe2 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -160,7 +160,7 @@ for i = 1:2
lb_= '\gamma';
end
errorbar(x_,y_,e_,'-ok',...
- 'LineWidth',1.5,...
+ 'LineWidth',1.5,'MarkerSize',10,...
'DisplayName',...
['$\gamma$, (',num2str(DATA.inputs.PMAX),',',num2str(DATA.inputs.JMAX),')']);
hold on;
@@ -172,8 +172,8 @@ for i = 1:2
end
% ylabel('$\gamma$');
% yyaxis("right")
- errorbar(k(MODES),imag(omega(ik)),imag(err(ik)),'--ok',...
- 'LineWidth',1.5,...
+ errorbar(k(MODES),imag(omega(ik)),imag(err(ik)),'*k',...
+ 'LineWidth',1.5,'MarkerSize',10,...
'DisplayName',...
['$\omega$, (',num2str(DATA.inputs.PMAX),',',num2str(DATA.inputs.JMAX),')']);
hold on;
diff --git a/wk/analysis_demo.m b/wk/analysis_demo.m
new file mode 100644
index 0000000000000000000000000000000000000000..11a31c9a4c6471b76d4be68315d7e5dab1ec8a7b
--- /dev/null
+++ b/wk/analysis_demo.m
@@ -0,0 +1,90 @@
+gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory
+addpath(genpath([gyacomodir,'matlab'])) % ... add
+addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
+addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
+addpath(genpath([gyacomodir,'matlab/load'])) % ... add
+default_plots_options
+options.RESOLUTION = 256;
+% Partition of the computer where the data have to be searched
+
+J0 = 00; J1 = 10;
+
+% DATADIR = '/Users/ahoffmann/gyacomo/testcases/DIII-D_triangularity_fast_nonlinear/';
+% DATADIR = '/Users/ahoffmann/gyacomo/testcases/HEL_DIII-D_triangularity/';
+DATADIR = '/Users/ahoffmann/gyacomo/testcases/cyclone_example/';
+data = {};
+data = compile_results_low_mem(data,DATADIR,J0,J1);
+
+if 1
+%% Plot transport and phi radial profile
+[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
+options.TAVG_0 = 100;
+options.TAVG_1 = 1000;
+options.NCUT = 5; % Number of cuts for averaging and error estimation
+options.NMVA = 1; % Moving average for time traces
+options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+options.INTERP = 1;
+plot_radial_transport_and_spacetime(data,options);
+end
+
+if 0
+%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Options
+[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
+[data.Na00, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Na00');
+data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+options.INTERP = 1;
+options.POLARPLOT = 0;
+options.BWR = 1; % bluewhitered plot or gray
+options.CLIMAUTO = 0; % adjust the colormap auto
+options.NAME = '\phi';
+% options.NAME = 'N_i^{00}';
+options.PLAN = 'xy';
+options.COMP = 9;
+options.TIME = data.Ts3D(1:2:end);
+data.EPS = 0.1;
+data.a = data.EPS * 2000;
+create_film(data,options,'.gif')
+end
+
+if 0
+%% field snapshots
+% Options
+[data.Na00, data.Ts3D] = compile_results_3Da(data.folder,J0,J1,'Na00');
+[data.PHI, data.Ts3D] = compile_results_3D(data.folder,J0,J1,'phi');
+data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+
+options.INTERP = 0;
+options.POLARPLOT = 0;
+options.AXISEQUAL = 0;
+options.NORMALIZE = 0;
+options.LOGSCALE = 0;
+options.CLIMAUTO = 1;
+options.NAME = 'N_i^{00}';
+% options.NAME = 's_{Ey}';
+% options.NAME = '\phi';
+options.PLAN = 'yz';
+options.COMP = 'avg';
+options.TIME = [10 30];
+fig = photomaton(data,options);
+colormap(gray)
+clim('auto')
+% save_figure(data,fig)
+end
+if 0
+%% Performance profiler
+profiler(data)
+end
+
+if 0
+%% Hermite-Laguerre spectrum
+[data.Napjz, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Napjz');
+% [data.Napjz, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'Nipjz');
+options.ST = 0;
+options.NORMALIZED = 0;
+options.LOGSCALE = 1;
+options.FILTER = 0; %filter the 50% time-average of the spectrum from
+options.TAVG_2D = 0; %Show a 2D plot of the modes, 50% time averaged
+options.TAVG_2D_CTR= 0; %make it contour plot
+fig = show_moments_spectrum(data,options);
+end
diff --git a/wk/parameters/lin_DIIID_LM_rho90.m b/wk/parameters/lin_DIIID_LM_rho90.m
index 9995af00da479279d9fd14fab07431aa0c7780e0..c699e5e9aa6505fb71c4320acdd6f6d80fa341b5 100644
--- a/wk/parameters/lin_DIIID_LM_rho90.m
+++ b/wk/parameters/lin_DIIID_LM_rho90.m
@@ -1,7 +1,7 @@
%% Reference values
% See Neiser et al. 2019 Gyrokinetic GENE simulations of DIII-D near-edge L-mode plasmas
%% Set simulation parameters
-SIMID = 'lin_DIID_LM_rho90'; % Name of the simulation
+SIMID = 'lin_DIIID_LM_rho90'; % Name of the simulation
%% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
NU = 1.0; %(0.00235 in GENE)
@@ -47,9 +47,9 @@ PB_PHASE = 0;
%% TIME PARAMETERS
TMAX = 10; % Maximal time unit
DT = 5e-4; % Time step
-DTSAVE0D = 0.5; % Sampling time for 0D arrays
+DTSAVE0D = 0.1; % Sampling time for 0D arrays
DTSAVE2D = -1; % Sampling time for 2D arrays
-DTSAVE3D = 0.5; % Sampling time for 3D arrays
+DTSAVE3D = 0.1; % Sampling time for 3D arrays
DTSAVE5D = 100; % Sampling time for 5D arrays
JOB2LOAD = -1; % Start a new simulation serie
diff --git a/wk/parameters/lin_GASTD.m b/wk/parameters/lin_GASTD.m
new file mode 100644
index 0000000000000000000000000000000000000000..14d13da71466a0a8489758cf671db0b36e89e9a4
--- /dev/null
+++ b/wk/parameters/lin_GASTD.m
@@ -0,0 +1,99 @@
+%% Reference values
+% See Staebler et al. 2023
+%% Set simulation parameters
+SIMID = 'lin_GASTD'; % Name of the simulation
+%% Set up physical parameters
+CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
+NU = 0.05; %(0.00235 in GENE)
+TAU = 1; % i/e temperature ratio
+K_Ne = 1; % ele Density '''
+K_Te = 3; % ele Temperature '''
+K_Ni = 1; % ion Density gradient drive
+K_Ti = 3; % ion Temperature '''
+SIGMA_E = 0.0233380;%/sqrt(2); % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+NA = 2; % number of kinetic species
+ADIAB_E = (NA==1); % adiabatic electron model
+BETA = 0; % electron plasma beta in prct
+MHD_PD = 1;
+%% Set up grid parameters
+P = 2;
+J = P/2;%P/2;
+PMAX = P; % Hermite basis size
+JMAX = J; % Laguerre basis size
+NX = 4; % real space x-gridpoints
+NY = 2; % real space y-gridpoints
+LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
+LY = 2*pi/0.5; % Size of the squared frequency domain in y direction
+NZ = 16; % number of perpendicular planes (parallel grid)
+SG = 0; % Staggered z grids option
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+
+%% GEOMETRY
+% GEOMETRY= 's-alpha';
+GEOMETRY= 'miller';
+Q0 = 2; % safety factor
+SHEAR = 1; % magnetic shear
+EPS = 1/3; % inverse aspect ratio
+KAPPA = 1.0; % elongation
+S_KAPPA = 0;
+DELTA = 0; % triangularity
+S_DELTA = 0;
+ZETA = 0; % squareness
+S_ZETA = 0;
+PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
+SHIFT_Y = 0.0; % Shift in the periodic BC in z
+NPOL = 1; % Number of poloidal turns
+PB_PHASE = 0;
+%% TIME PARAMETERS
+TMAX = 20; % Maximal time unit
+DT = 1e-2; % Time step
+DTSAVE0D = 1.0; % Sampling time for 0D arrays
+DTSAVE2D = -1; % Sampling time for 2D arrays
+DTSAVE3D = 1.0; % Sampling time for 3D arrays
+DTSAVE5D = 100; % Sampling time for 5D arrays
+JOB2LOAD = -1; % Start a new simulation serie
+
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+GKCO = 1; % Gyrokinetic operator
+ABCO = 1; % INTERSPECIES collisions
+INIT_ZF = 0; % Initialize zero-field quantities
+HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+DMAX = -1;
+NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+NMAX = 0;
+KERN = 0; % Kernel model (0 : GK)
+INIT_OPT = 'phi'; % Start simulation with a noisy mom00/phi/allmom
+NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5)
+
+%% OUTPUTS
+W_DOUBLE = 1; % Output flag for double moments
+W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy)
+W_HF = 1; % Output flag for high-frequency potential energy
+W_PHI = 1; % Output flag for potential
+W_NA00 = 1; % Output flag for nalpha00 (density of species alpha)
+W_DENS = 1; % Output flag for total density
+W_TEMP = 1; % Output flag for temperature
+W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha)
+W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% UNUSED PARAMETERS
+% These parameters are usually not to play with in linear runs
+MU = 0.0; % Hyperdiffusivity coefficient
+MU_X = MU; % Hyperdiffusivity coefficient in x direction
+MU_Y = MU; % Hyperdiffusivity coefficient in y direction
+N_HD = 4; % Degree of spatial-hyperdiffusivity
+MU_Z = 5.0; % Hyperdiffusivity coefficient in z direction
+HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model
+MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
+MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
+LAMBDAD = 0.0; % Lambda Debye
+NOISE0 = 1.0e-5; % Initial noise amplitude
+BCKGD0 = 0.0e-5; % Initial background
+k_gB = 1.0; % Magnetic gradient strength
+k_cB = 1.0; % Magnetic curvature strength
+COLL_KCUT = 1; % Cutoff for collision operator
+ADIAB_I = 0; % adiabatic ion model
\ No newline at end of file