diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index dc515acd2f7387742433b8d53e2cb9397aa673ea..aedbeb58671263694fc373ef3582f0a81f0ae901 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -45,8 +45,8 @@ while(CONTINUE)
W_SAPJ = strcmp(h5readatt(filename,'/data/input','write_Sapj' ),'y');
W_DENS = strcmp(h5readatt(filename,'/data/input','write_dens' ),'y');
W_TEMP = strcmp(h5readatt(filename,'/data/input','write_temp' ),'y');
-% KIN_E = strcmp(h5readatt(filename,'/data/input', 'KIN_E' ),'y');
- KIN_E = 1;
+ KIN_E = strcmp(h5readatt(filename,'/data/input', 'KIN_E' ),'y');
+% KIN_E = 1;
% Check polynomials degrees
Pe_new= numel(Pe); Je_new= numel(Je);
diff --git a/matlab/create_film.m b/matlab/create_film.m
index b48cf206c97d4afbb3e3549b13b6507829c3947d..5f9724726d2821d1ba205ca01925bc59a4e6d225 100644
--- a/matlab/create_film.m
+++ b/matlab/create_film.m
@@ -1,13 +1,13 @@
function create_film(DATA,OPTIONS,format)
%% Plot options
FPS = 30; DELAY = 1/FPS;
-INTERP = OPTIONS.INTERP; BWR = 1; NORMALIZED = 1;
+BWR = 1; NORMALIZED = 1;
T = DATA.Ts3D;
%% Processing
toplot = process_field(DATA,OPTIONS);
%%
-FILENAME = [toplot.FILENAME,format];
+FILENAME = [DATA.localdir,toplot.FILENAME,format];
XNAME = ['$',toplot.XNAME,'$'];
YNAME = ['$',toplot.YNAME,'$'];
FIELDNAME = ['$',toplot.FIELDNAME,'$'];
@@ -28,13 +28,13 @@ if hmax == hmin
disp('Warning : h = hmin = hmax = const')
else
% Setup figure frame
-fig = figure('Color','white','Position', toplot.DIMENSIONS);
+fig = figure('Color','white','Position', toplot.DIMENSIONS.*[1 1 1.2 1]);
pcolor(X,Y,FIELD(:,:,1)); % to set up
if BWR
colormap(bluewhitered)
end
axis tight manual % this ensures that getframe() returns a consistent size
- if INTERP
+ if toplot.INTERP
shading interp;
end
in = 1;
@@ -56,7 +56,7 @@ fig = figure('Color','white','Position', toplot.DIMENSIONS);
end
title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
,', scaling = ',sprintf('%.1e',scale)]);
- if INTERP
+ if toplot.INTERP
shading interp;
end
set(pclr, 'edgecolor','none'); pbaspect(toplot.ASPECT);
diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index 32ac436755659dbdcee89b79ab50d3d9f5e3612d..fa28a51fe520be2f73dfaa26d93940bf92844994 100644
--- a/matlab/extract_fig_data.m
+++ b/matlab/extract_fig_data.m
@@ -4,7 +4,7 @@ dataObjs = axObjs.Children;
X_ = dataObjs(1).XData;
Y_ = dataObjs(1).YData;
-
+n0 = 1500;
figure;
-plot(X_,Y_);
-plot(X_-X_(1),Y_);
\ No newline at end of file
+plot(X_(n0:end),Y_(n0:end));
+plot(X_(n0:end)-X_(n0),Y_(n0:end));
\ No newline at end of file
diff --git a/matlab/load_params.m b/matlab/load_params.m
index 7042b6c86376426fe37a03e7e76b3dd284220607..62ac961b49b1e113ee9f2744c6eb96c89ac01b04 100644
--- a/matlab/load_params.m
+++ b/matlab/load_params.m
@@ -1,9 +1,9 @@
DATA.CO = h5readatt(filename,'/data/input','CO');
-DATA.K_N = h5readatt(filename,'/data/input','eta_n');
-DATA.K_T = h5readatt(filename,'/data/input','eta_T');
-% DATA.K_N = h5readatt(filename,'/data/input','K_n');
-% DATA.K_T = h5readatt(filename,'/data/input','K_T');
-% DATA.K_E = h5readatt(filename,'/data/input','K_E');
+% DATA.K_N = h5readatt(filename,'/data/input','eta_n');
+% DATA.K_T = h5readatt(filename,'/data/input','eta_T');
+DATA.K_N = h5readatt(filename,'/data/input','K_n');
+DATA.K_T = h5readatt(filename,'/data/input','K_T');
+DATA.K_E = h5readatt(filename,'/data/input','K_E');
DATA.Q0 = h5readatt(filename,'/data/input','q0');
DATA.SHEAR = h5readatt(filename,'/data/input','shear');
DATA.EPS = h5readatt(filename,'/data/input','eps');
diff --git a/matlab/photomaton.m b/matlab/photomaton.m
index f054b94dcc7611d6ca3886c740515910a4d9fa7b..f1c8d13cbe470a73e9204c551c2b7e745000844b 100644
--- a/matlab/photomaton.m
+++ b/matlab/photomaton.m
@@ -10,10 +10,19 @@ TNAME = [];
FIGURE.fig = figure; set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 numel(FRAMES) 1])
for i_ = 1:numel(FRAMES)
subplot(1,numel(FRAMES),i_); TNAME = [TNAME,'_',sprintf('%.0f',DATA.Ts3D(FRAMES(i_)))];
- pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,FRAMES(i_))); set(pclr, 'edgecolor','none');pbaspect(toplot.ASPECT)
- colormap(bluewhitered);
- xlabel(toplot.XNAME); ylabel(toplot.YNAME);set(gca,'ytick',[]);
- title(sprintf('$t c_s/R=%.0f$',DATA.Ts3D(FRAMES(i_))));
+ if ~strcmp(OPTIONS.PLAN,'kxky')
+ scale = max(max(abs(toplot.FIELD(:,:,FRAMES(i_))))); % Scaling to normalize
+ else
+ scale = 1;
+ end
+ pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,FRAMES(i_))./scale); set(pclr, 'edgecolor','none');pbaspect(toplot.ASPECT)
+ if ~strcmp(OPTIONS.PLAN,'kxky')
+ caxis([-1,1]);
+ colormap(bluewhitered);
+ end
+ xlabel(toplot.XNAME); ylabel(toplot.YNAME);
+% if i_ > 1; set(gca,'ytick',[]); end;
+ title([sprintf('$t c_s/R=%.0f$',DATA.Ts3D(FRAMES(i_))),', max = ',sprintf('%.1e',scale)]);
if OPTIONS.INTERP
shading interp;
end
diff --git a/matlab/process_field.m b/matlab/process_field.m
index 563694e7dbd4b7194d191f643372ca80929e4fad..68433ff5385e91b7e5b9b453c432bf3c470671dc 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -98,11 +98,12 @@ end
switch REALP
case 1 % Real space plot
+ INTERP = OPTIONS.INTERP;
process = @(x) real(fftshift(ifft2(x,Nx,Ny)));
shift_x = @(x) x;
shift_y = @(x) x;
case 0 % Frequencies plot
- OPTIONS.INTERP = 0;
+ INTERP = 0;
switch COMPDIM
case 1
process = @(x) abs(fftshift(x,2));
@@ -180,6 +181,58 @@ switch OPTIONS.NAME
FIELD(:,:,it) = squeeze(compr(tmp));
end
end
+ case '\Gamma_y'
+ NAME = 'Gy';
+ [KY, ~] = meshgrid(DATA.ky, DATA.kx);
+ Gx = zeros(DATA.Nx,DATA.Ny,DATA.Nz,numel(FRAMES));
+ for it = FRAMES % Compute the 3D real transport for each frame
+ for iz = 1:DATA.Nz
+ Gx(:,:,iz,it) = real(ifft2(-1i*KY.*(DATA.PHI(:,:,iz,it)),DATA.Nx,DATA.Ny))...
+ .*real(ifft2(DATA.DENS_I(:,:,iz,it),DATA.Nx,DATA.Ny));
+ end
+ end
+ if REALP % plot in real space
+ for it = FRAMES
+ FIELD(:,:,it) = squeeze(compr(Gx(:,:,:,it)));
+ end
+ else % Plot spectrum
+ process = @(x) abs(fftshift(x,2));
+ shift_x = @(x) fftshift(x,2);
+ shift_y = @(x) fftshift(x,2);
+ tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ for it = FRAMES
+ for iz = 1:numel(DATA.z)
+ tmp(:,:,iz) = process(squeeze(fft2(Gx(:,:,iz,it),DATA.Nkx,DATA.Nky)));
+ end
+ FIELD(:,:,it) = squeeze(compr(tmp));
+ end
+ end
+ case '\Gamma_x'
+ NAME = 'Gx';
+ [KY, ~] = meshgrid(DATA.ky, DATA.kx);
+ Gx = zeros(DATA.Nx,DATA.Ny,DATA.Nz,numel(FRAMES));
+ for it = FRAMES % Compute the 3D real transport for each frame
+ for iz = 1:DATA.Nz
+ Gx(:,:,iz,it) = real((ifft2(-1i*KY.*(DATA.PHI(:,:,iz,it)),DATA.Nx,DATA.Ny)))...
+ .*real((ifft2(DATA.DENS_I(:,:,iz,it),DATA.Nx,DATA.Ny)));
+ end
+ end
+ if REALP % plot in real space
+ for it = FRAMES
+ FIELD(:,:,it) = squeeze(compr(Gx(:,:,:,it)));
+ end
+ else % Plot spectrum
+ process = @(x) abs(fftshift(x,2));
+ shift_x = @(x) fftshift(x,2);
+ shift_y = @(x) fftshift(x,2);
+ tmp = zeros(DATA.Nx,DATA.Ny,Nz);
+ for it = FRAMES
+ for iz = 1:numel(DATA.z)
+ tmp(:,:,iz) = process(squeeze(fft2(Gx(:,:,iz,it),DATA.Nx,DATA.Ny)));
+ end
+ FIELD(:,:,it) = squeeze(compr(tmp(1:DATA.Nkx,1:DATA.Nky,:)));
+ end
+ end
end
TOPLOT.FIELD = FIELD;
TOPLOT.X = shift_x(X);
@@ -191,6 +244,6 @@ TOPLOT.FILENAME = [NAME,'_',OPTIONS.PLAN,'_',COMPNAME,'_',POLARNAME];
TOPLOT.DIMENSIONS= DIMENSIONS;
TOPLOT.ASPECT = ASPECT;
TOPLOT.FRAMES = FRAMES;
-
+TOPLOT.INTERP = INTERP;
end
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index a73226576e53801425b2048d9c457dca85cb1d55..189db517b4c2b91d5090ace6a45cd3f94f959da3 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -4,9 +4,8 @@ outfile ='';
%% Directory of the simulation
if 1% Local results
outfile ='';
-outfile ='';
-outfile ='';
-outfile ='simulation_A/CO_damping_FCGK';
+outfile ='Cyclone/100x100x30_5x3_Lx_200_Ly_100_q0_1.4_e_0.18_kN_2.22_kT_6.9_nu_1e-02_DGGK_adiabe';
+% outfile ='simulation_A/CO_damping_FCGK';
% outfile ='fluxtube_salphaB_s0/100x100x30_5x3_Lx_200_Ly_100_q0_2.7_e_0.18_kN_2.22_kT_8_nu_1e-01_DGGK_adiabe';
% outfile ='fluxtube_salphaB_s0/50x100x20_5x3_L_300_q0_2.7_e_0.18_kN_2.22_kT_8_nu_1e-01_DGGK_adiabe';
% outfile ='fluxtube_salphaB_s0/50x100x20_5x3_L_300_q0_2.7_e_0.18_kN_2.22_kT_8_nu_1e-01_DGGK_adiabe_Sg';
@@ -16,15 +15,19 @@ outfile ='simulation_A/CO_damping_FCGK';
CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
system(CMD);
else% Marconi results
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x200_L_200_P_8_J_4_eta_0.6_nu_1e-01_PAGK_mu_0e+00/out.txt';
-% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x300_L_120_P_8_J_4_eta_0.6_nu_1e-01_PAGK_mu_0e+00/out.txt';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x200_L_200_P_8_J_4_eta_0.6_nu_1e-01_PAGK_mu_0e+00/out.txt';
+outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/500x500_L_120_P_4_J_2_eta_0.6_nu_1e-01_DGGK_mu_0e+00/out.txt';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/cw_FCGK_kp_3.0/out.txt';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/nonlin_FCGK/150x75_L_200_P_4_J_2_eta_0.6_nu_1e-01_FCGK_mu_0e+00/out.txt';
+
% BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
-MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
+MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'^{NZ}/'];
end
%% Load the results
% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 07; JOBNUMMAX = 07;
+JOBNUMMIN = 00; JOBNUMMAX = 13;
data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
@@ -35,7 +38,7 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-TAVG_0 = 1500; TAVG_1 = 5000; % Averaging times duration
+TAVG_0 = 100; TAVG_1 = 500; % Averaging times duration
fig = plot_radial_transport_and_shear(data,TAVG_0,TAVG_1);
save_figure(data,fig)
end
@@ -46,10 +49,11 @@ if 0
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
+% options.NAME = '\Gamma_x';
options.NAME = 'n_i^{NZ}';
options.PLAN = 'xy';
-options.COMP = 1;
-options.TIME = 00:20:6000;
+options.COMP = 16;
+options.TIME = 0:1:data.Ts3D(end);
% options.TIME = 140:0.5:160;
data.a = data.EPS * 2000;
create_film(data,options,'.gif')
@@ -60,10 +64,11 @@ if 0
% Options
options.INTERP = 0;
options.POLARPLOT = 0;
+% options.NAME = '\Gamma_x';
options.NAME = 'n_i^{NZ}';
-options.PLAN = 'xy';
+options.PLAN = 'kxky';
options.COMP = 1;
-options.TIME = [100 400 2500 5500];
+options.TIME = [50 800 1200];
data.a = data.EPS * 1000;
fig = photomaton(data,options);
save_figure(data,fig)
@@ -73,8 +78,8 @@ if 0
%% 3D plot on the geometry
options.INTERP = 1;
options.NAME = 'n_i';
-options.PLANES = 1;
-options.TIME = 5000;
+options.PLANES = 16;
+options.TIME = 50;
data.rho_o_R = 1e-3; % Sound larmor radius over Machine size ratio
FIGURE = show_geometry(data,options);
end
@@ -85,11 +90,6 @@ TAVG_0 = 1000; TAVG_1 = 5000; % Averaging times duration
ZF_fourier_analysis
end
-if 0
-%%
-plot_param_evol
-end
-
if 0
%% Radial shear profile (with moving average)
tf = 1000+[0:100:1000];
diff --git a/wk/linear_1D_damping.m b/wk/linear_1D_damping.m
new file mode 100644
index 0000000000000000000000000000000000000000..95e0fca00b9fb8792266f52f3a4ed7ba36520e37
--- /dev/null
+++ b/wk/linear_1D_damping.m
@@ -0,0 +1,156 @@
+RUN = 1; % To run or just to load
+addpath(genpath('../matlab')) % ... add
+default_plots_options
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Set Up parameters
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% PHYSICAL PARAMETERS
+NU = 0.1; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+%% GRID PARAMETERS
+NX = 150; % real space x-gridpoints
+NY = 1; % '' y-gridpoints
+LX = 200; % Size of the squared frequency domain
+LY = 1; % Size of the squared frequency domain
+NZ = 1; % number of perpendicular planes (parallel grid)
+Q0 = 1.0; % safety factor
+SHEAR = 0.0; % magnetic shear
+EPS = 0.0; % inverse aspect ratio
+SG = 0; % Staggered z grids option
+%% TIME PARMETERS
+TMAX = 20; % Maximal time unit
+DT = 1e-2; % Time step
+SPS0D = 1; % Sampling per time unit for 2D arrays
+SPS2D = 0; % Sampling per time unit for 2D arrays
+SPS3D = 5; % Sampling per time unit for 2D arrays
+SPS5D = 1/5; % Sampling per time unit for 5D arrays
+SPSCP = 0; % Sampling per time unit for checkpoints
+JOB2LOAD= -1;
+%% OPTIONS
+NOISE0 = 0.0; % Init noise amplitude
+BCKGD0 = 1.0; % Init background
+SIMID = 'Linear_damping'; % Name of the simulation
+NON_LIN = 0; % activate non-linearity (is cancelled if KXEQ0 = 1)
+KIN_E = 1;
+% Collision operator
+% (0:L.Bernstein, 1:Dougherty, 2:Sugama, 3:Pitch angle, 4:Full Couloumb ; +/- for GK/DK)
+CO = 4;
+INIT_ZF = 0; ZF_AMP = 0.0;
+CLOS = 0; % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))
+NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
+KERN = 0; % Kernel model (0 : GK)
+INIT_PHI= 0; % Start simulation with a noisy phi
+%% OUTPUTS
+W_DOUBLE = 1;
+W_GAMMA = 1; W_HF = 1;
+W_PHI = 1; W_NA00 = 1;
+W_DENS = 1; W_TEMP = 1;
+W_NAPJ = 1; W_SAPJ = 0;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% unused
+HD_CO = 0.5; % Hyper diffusivity cutoff ratio
+kmax = NX*pi/LX;% Highest fourier mode
+NU_HYP = 0.0; % Hyperdiffusivity coefficient
+MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
+INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0;
+MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
+MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
+LAMBDAD = 0.0;
+GRADB = 0.0;
+CURVB = 0.0;
+K_N = 0.0; % Density gradient drive
+K_T = 0.0; % Temperature '''
+K_E = 0.0; % Electrostat '''
+%% PARAMETER SCANS
+
+if 1
+%% Parameter scan over PJ
+% PA = [2 4];
+% JA = [1 2];
+PA = [2];
+JA = [1];
+DTA= DT*ones(size(JA));%./sqrt(JA);
+% DTA= DT;
+mup_ = MU_P;
+muj_ = MU_J;
+Nparam = numel(PA);
+param_name = 'PJ';
+gamma_Ni00 = zeros(Nparam,floor(NX/2)+1);
+gamma_Nipj = zeros(Nparam,floor(NX/2)+1);
+gamma_phi = zeros(Nparam,floor(NX/2)+1);
+for i = 1:Nparam
+ % Change scan parameter
+ PMAXE = PA(i); PMAXI = PA(i);
+ JMAXE = JA(i); JMAXI = JA(i);
+ DT = DTA(i);
+ setup
+ system(['rm fort*.90']);
+ % Run linear simulation
+ if RUN
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz3 1 6 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ./../../../bin/helaz3 1 2 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz3 0; cd ../../../wk'])
+ end
+% Load and process results
+ %%
+ filename = ['../results/',SIMID,'/',PARAMS,'/outputs_00.h5'];
+ load_results
+ for ikx = 1:NX/2+1
+ tend = 2;%max(Ts3D(abs(Ni00(ikx,1,1,:))~=0));
+ tstart = 0;
+ [~,itstart] = min(abs(Ts3D-tstart));
+ [~,itend] = min(abs(Ts3D-tend));
+ trange = itstart:itend;
+ % exp fit on moment 00
+ X_ = Ts3D(trange); Y_ = squeeze(abs(Ni00(ikx,1,1,trange)));
+ gamma_Ni00(i,ikx) = LinearFit_s(X_,Y_);
+ % exp fit on phi
+ X_ = Ts3D(trange); Y_ = squeeze(abs(PHI(ikx,1,1,trange)));
+ gamma_phi (i,ikx) = LinearFit_s(X_,Y_);
+ end
+ gamma_Ni00(i,:) = real(gamma_Ni00(i,:));% .* (gamma_Ni00(i,:)>=0.0));
+ gamma_Nipj(i,:) = real(gamma_Nipj(i,:));% .* (gamma_Nipj(i,:)>=0.0));
+ if 0
+ %% Fit verification
+ figure;
+ for i = 1:1:NX/2+1
+ X_ = Ts3D(:); Y_ = squeeze(abs(Ni00(i,1,1,:)));
+ semilogy(X_,Y_,'DisplayName',['k=',num2str(kx(i))]); hold on;
+ end
+end
+
+if 1
+%% Plot
+SCALE = 1;%sqrt(2);
+fig = figure; FIGNAME = 'linear_study';
+plt = @(x) x;
+% subplot(211)
+ for i = 1:Nparam
+ clr = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:);
+ linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1);
+ plot(plt(SCALE*kx),plt(gamma_Ni00(i,1:end)),...
+ 'Color',clr,...
+ 'LineStyle',linestyle{1},'Marker','^',...
+...% 'DisplayName',['$\kappa_N=',num2str(K_N),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, $P=',num2str(PA(i)),'$, $J=',num2str(JA(i)),'$']);
+ 'DisplayName',[CONAME,', $P,J=',num2str(PA(i)),',',num2str(JA(i)),'$']);
+ hold on;
+ end
+ grid on; xlabel('$k_y\rho_s^{R}$'); ylabel('$\gamma(\phi)L_\perp/c_s$'); xlim([0.0,max(kx)]);
+ title(['$\kappa_N=',num2str(K_N),'$, $\nu_{',CONAME,'}=',num2str(NU),'$'])
+ legend('show'); %xlim([0.01,10])
+saveas(fig,[SIMDIR,'/',PARAMS,'/gamma_vs_',param_name,'_',PARAMS,'.fig']);
+saveas(fig,[SIMDIR,'/',PARAMS,'/gamma_vs_',param_name,'_',PARAMS,'.png']);
+end
+end
+if 0
+%% Space time
+ [YT,XT] = meshgrid(Ts3D,kx);
+ figure;
+% pclr = surf(XT,YT,squeeze(abs(PHI_ST(1,:,:)))); set(pclr, 'edgecolor','none'); colorbar;
+% pclr = pcolor(XT,YT,squeeze(abs(Ni00_ST(1,:,:)))); set(pclr, 'edgecolor','none'); colorbar;
+ semilogy(Ts3D(1:TMAX/SPS3D),squeeze(abs(PHI_ST(1,50:5:100,:))));
+end
+end
\ No newline at end of file
diff --git a/wk/local_run.m b/wk/local_run.m
index 5490687b4d0922ff40419e082bd914aba4152513..7e7776b2ec10d5fd3ccab401dcb0a8e209899c7f 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -4,12 +4,12 @@ addpath(genpath('../matlab')) % ... add
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.1; % Collision frequency
+NU = 0.01; % Collision frequency
K_N = 2.22; % Density gradient drive
-K_T = 8.0; % Temperature '''
+K_T = 6.9; % Temperature '''
K_E = 0.00; % Electrostat gradient
SIGMA_E = 0.05196; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
-NU_HYP = 0.0;
+NU_HYP = 0.01;
KIN_E = 0; % Kinetic (1) or adiabatic (2) electron model
%% GRID PARAMETERS
NX = 100; % Spatial radial resolution ( = 2x radial modes)
@@ -20,18 +20,18 @@ NZ = 30; % number of perpendicular planes (parallel grid)
P = 4;
J = 2;
%% GEOMETRY PARAMETERS
-Q0 = 2.7; % safety factor
+Q0 = 1.4; % safety factor
SHEAR = 0.0; % magnetic shear
EPS = 0.18; % inverse aspect ratio
GRADB = 1.0; % Magnetic gradient
CURVB = 1.0; % Magnetic curvature
SG = 1; % Staggered z grids option
%% TIME PARAMETERS
-TMAX = 200; % Maximal time unit
+TMAX = 500; % Maximal time unit
DT = 5e-2; % Time step
-SPS0D = 1; % Sampling per time unit for profiler
+SPS0D = 2; % Sampling per time unit for profiler
SPS2D = 1; % Sampling per time unit for 2D arrays
-SPS3D = 5; % Sampling per time unit for 3D arrays
+SPS3D = 1/2; % Sampling per time unit for 3D arrays
SPS5D = 1/200; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints/10
JOB2LOAD= -1;
@@ -41,7 +41,7 @@ JOB2LOAD= -1;
CO = 1;
CLOS = 0; % Closure model (0: =0 truncation)
NL_CLOS = 0; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
-SIMID = 'fluxtube_salphaB_s0'; % Name of the simulation
+SIMID = 'Cyclone'; % Name of the simulation
% SIMID = 'simulation_A'; % Name of the simulation
% SIMID = ['v3.0_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation
NON_LIN = 1; % activate non-linearity (is cancelled if KXEQ0 = 1)