diff --git a/Makefile b/Makefile
index 241b4d4a9fed7282b5372db050ebf3547d7e4014..8901e87be09a9e8f6f7f48f4ef13585f52661842 100644
--- a/Makefile
+++ b/Makefile
@@ -91,7 +91,7 @@ $(OBJDIR)/utility_mod.o
$(OBJDIR)/closure_mod.o : src/closure_mod.F90 $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/closure_mod.F90 -o $@
- $(OBJDIR)/collision_mod.o : src/collision_mod.F90 $(OBJDIR)/initial_par_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
+ $(OBJDIR)/collision_mod.o : src/collision_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/collision_mod.F90 -o $@
$(OBJDIR)/compute_Sapj.o : src/compute_Sapj.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fourier_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o
@@ -133,7 +133,7 @@ $(OBJDIR)/utility_mod.o
$(OBJDIR)/main.o : src/main.F90 $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/main.F90 -o $@
- $(OBJDIR)/memory.o : src/memory.F90 $ $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/grid_mod.o
+ $(OBJDIR)/memory.o : src/memory.F90 $ $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/collision_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/grid_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/memory.F90 -o $@
$(OBJDIR)/model_mod.o : src/model_mod.F90 $(OBJDIR)/prec_const_mod.o
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index aedbeb58671263694fc373ef3582f0a81f0ae901..2bb4b5ef8f08bc0555c77bb66bb1d49ff977492a 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -41,8 +41,8 @@ while(CONTINUE)
W_HF = 0;%strcmp(h5readatt(filename,'/data/input','write_hf' ),'y');
W_PHI = strcmp(h5readatt(filename,'/data/input','write_phi' ),'y');
W_NA00 = strcmp(h5readatt(filename,'/data/input','write_Na00' ),'y');
- W_NAPJ = strcmp(h5readatt(filename,'/data/input','write_Napj' ),'y');
- W_SAPJ = strcmp(h5readatt(filename,'/data/input','write_Sapj' ),'y');
+ W_NAPJ = 0;%strcmp(h5readatt(filename,'/data/input','write_Napj' ),'y');
+ W_SAPJ = 0;%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');
@@ -135,6 +135,7 @@ while(CONTINUE)
TEMP_I_ = cat(4,TEMP_I_,TEMP_I); clear TEMP_I
end
+ Ts5D = [];
if W_NAPJ
[Nipj, Ts5D, ~] = load_5D_data(filename, 'moments_i');
Nipj_ = cat(6,Nipj_,Nipj); clear Nipj
diff --git a/matlab/load_params.m b/matlab/load_params.m
index 62ac961b49b1e113ee9f2744c6eb96c89ac01b04..1d44a63a5dac5bd3c4457f27dc7337dc110c6ac5 100644
--- a/matlab/load_params.m
+++ b/matlab/load_params.m
@@ -11,7 +11,8 @@ DATA.PMAXI = h5readatt(filename,'/data/input','pmaxi');
DATA.JMAXI = h5readatt(filename,'/data/input','jmaxi');
DATA.PMAXE = h5readatt(filename,'/data/input','pmaxe');
DATA.JMAXE = h5readatt(filename,'/data/input','jmaxe');
-DATA.NON_LIN = h5readatt(filename,'/data/input','NON_LIN');
+% DATA.LINEARITY = h5readatt(filename,'/data/input','NON_LIN');
+% DATA.LINEARITY = h5readatt(filename,'/data/input','LINEARITY');
DATA.NU = h5readatt(filename,'/data/input','nu');
DATA.Nx = h5readatt(filename,'/data/input','Nx');
DATA.Ny = h5readatt(filename,'/data/input','Ny');
@@ -26,23 +27,13 @@ DATA.W_NA00 = h5readatt(filename,'/data/input','write_Na00') == 'y';
DATA.W_NAPJ = h5readatt(filename,'/data/input','write_Napj') == 'y';
DATA.W_SAPJ = h5readatt(filename,'/data/input','write_Sapj') == 'y';
-if DATA.NON_LIN == 'y'
- DATA.NON_LIN = 1;
-else
- DATA.NON_LIN = 0;
-end
+% if DATA.LINEARITY == 'y'
+% DATA.LINEARITY = 1;
+% else
+% DATA.LINEARITY = 0;
+% end
-switch abs(DATA.CO)
- case 0; DATA.CONAME = 'LB';
- case 1; DATA.CONAME = 'DG';
- case 2; DATA.CONAME = 'SG';
- case 3; DATA.CONAME = 'PA';
- case 4; DATA.CONAME = 'FC';
- otherwise; DATA.CONAME ='UK';
-end
-if (DATA.CO <= 0); DATA.CONAME = [DATA.CONAME,'DK'];
-else; DATA.CONAME = [DATA.CONAME,'GK'];
-end
+DATA.CONAME = DATA.CO;
if (DATA.CLOS == 0); DATA.CLOSNAME = 'Trunc.';
elseif(DATA.CLOS == 1); DATA.CLOSNAME = 'Clos. 1';
@@ -57,7 +48,7 @@ end
degngrad = [degngrad,'_Kn_%1.1f_nu_%0.0e_',...
DATA.CONAME,'_CLOS_',num2str(DATA.CLOS),'_mu_%0.0e'];
degngrad = sprintf(degngrad,[DATA.K_N,DATA.NU,DATA.MU]);
-if ~DATA.NON_LIN; degngrad = ['lin_',degngrad]; end
+% if ~DATA.LINEARITY; degngrad = ['lin_',degngrad]; end
resolution = [num2str(DATA.Nx),'x',num2str(DATA.Ny),'_'];
gridname = ['L_',num2str(DATA.L),'_'];
DATA.PARAMS = [resolution,gridname,degngrad];
\ No newline at end of file
diff --git a/matlab/photomaton.m b/matlab/photomaton.m
index f1c8d13cbe470a73e9204c551c2b7e745000844b..ef43d8740afb61bafdce67a37b4e654f9f2153dc 100644
--- a/matlab/photomaton.m
+++ b/matlab/photomaton.m
@@ -16,10 +16,10 @@ FIGURE.fig = figure; set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 numel(FRAMES)
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
+% 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)]);
diff --git a/matlab/process_field.m b/matlab/process_field.m
index 68433ff5385e91b7e5b9b453c432bf3c470671dc..b2df3a5f5477454497307b3dd3270c62e46af163 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -34,12 +34,10 @@ switch OPTIONS.PLAN
XNAME = '$k_x$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.kx);
REALP = 0; COMPDIM = 2; POLARPLOT = 0; SCALE = 0;
- OPTIONS.COMP = 'max';
case 'kyz'
XNAME = '$k_y$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.ky);
REALP = 0; COMPDIM = 1; POLARPLOT = 0; SCALE = 0;
- OPTIONS.COMP = 'max';
end
Xmax_ = max(max(abs(X))); Ymax_ = max(max(abs(Y)));
Lmin_ = min([Xmax_,Ymax_]);
@@ -81,17 +79,25 @@ switch OPTIONS.COMP
case 1
i = OPTIONS.COMP;
compr = @(x) x(i,:,:);
- COMPNAME = sprintf([directions{COMPDIM},'=','%2.1f'],DATA.y(i));
+ if REALP
+ COMPNAME = sprintf(['x=','%2.1f'],DATA.x(i));
+ else
+ COMPNAME = sprintf(['k_x=','%2.1f'],DATA.kx(i));
+ end
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
case 2
i = OPTIONS.COMP;
compr = @(x) x(:,i,:);
- COMPNAME = sprintf([directions{COMPDIM},'=','%2.1f'],DATA.y(i));
+ if REALP
+ COMPNAME = sprintf(['y=','%2.1f'],DATA.y(i));
+ else
+ COMPNAME = sprintf(['k_y=','%2.1f'],DATA.ky(i));
+ end
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
case 3
i = OPTIONS.COMP;
compr = @(x) x(:,:,i);
- COMPNAME = sprintf([directions{COMPDIM},'=','%2.1f'],DATA.z(i)/pi);
+ COMPNAME = sprintf(['z=','%2.1f'],DATA.z(i));
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
end
end
@@ -181,7 +187,27 @@ switch OPTIONS.NAME
FIELD(:,:,it) = squeeze(compr(tmp));
end
end
- case '\Gamma_y'
+ case '\phi^{NZ}'
+ NAME = 'phiNZ';
+ if COMPDIM == 3
+ for it = FRAMES
+ tmp = squeeze(compr(DATA.PHI(:,:,:,it).*(KY~=0)));
+ FIELD(:,:,it) = squeeze(process(tmp));
+ end
+ else
+ if REALP
+ tmp = zeros(Nx,Ny,Nz);
+ else
+ tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ end
+ for it = FRAMES
+ for iz = 1:numel(DATA.z)
+ tmp(:,:,iz) = squeeze(process(DATA.PHI(:,:,iz,it).*(KY~=0)));
+ end
+ 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));
diff --git a/matlab/setup.m b/matlab/setup.m
index b972a0c2deb4ed0bfd53a160c8dd28137efaafc4..f696387e5cf5109889e5d32c83f20bb6fc288d57 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -15,10 +15,9 @@ GRID.shear = SHEAR; % shear
GRID.eps = EPS; % inverse aspect ratio
if SG; GRID.SG = '.true.'; else; GRID.SG = '.false.';end;
% Model parameters
-MODEL.CO = CO; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
MODEL.CLOS = CLOS;
MODEL.NL_CLOS = NL_CLOS;
-MODEL.NON_LIN = NON_LIN;
+MODEL.LINEARITY = LINEARITY;
MODEL.KIN_E = KIN_E;
if KIN_E; MODEL.KIN_E = '.true.'; else; MODEL.KIN_E = '.false.';end;
MODEL.mu = MU;
@@ -42,42 +41,39 @@ MODEL.K_E = K_E; % Electric
MODEL.GradB = GRADB; % Magnetic gradient
MODEL.CurvB = CURVB; % Magnetic curvature
MODEL.lambdaD = LAMBDAD;
-% if A0KH ~= 0; SIMID = [SIMID,'_Nz_',num2str(L/2/pi*KX0KH),'_A_',num2str(A0KH)]; end;
+% Collision parameters
+COLL.collision_model = CO;
+if (GKCO); COLL.gyrokin_CO = '.true.'; else; COLL.gyrokin_CO = '.false.';end;
+if (ABCO); COLL.interspecies = '.true.'; else; COLL.interspecies = '.false.';end;
+COLL.mat_file = '''null''';
+switch CO
+ case 'SG'
+ COLL.mat_file = '''../../../iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5''';
+ case 'LR'
+ COLL.mat_file = '''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5''';
+ case 'LD'
+ COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_75_kpm_6.0.h5''';
+end
% Time integration and intialization parameters
TIME_INTEGRATION.numerical_scheme = '''RK4''';
if (INIT_PHI); INITIAL.init_noisy_phi = '.true.'; else; INITIAL.init_noisy_phi = '.false.';end;
INITIAL.INIT_ZF = INIT_ZF;
INITIAL.wipe_turb = WIPE_TURB;
-INITIAL.wipe_zf = WIPE_ZF;
+INITIAL.ACT_ON_MODES = ACT_ON_MODES;
if (INIT_BLOB); INITIAL.init_blob = '.true.'; else; INITIAL.init_blob = '.false.';end;
INITIAL.init_background = (INIT_ZF>0)*ZF_AMP + BCKGD0;
INITIAL.init_noiselvl = NOISE0;
INITIAL.iseed = 42;
-INITIAL.mat_file = '''null''';
-if (abs(CO) == 2) %Sugama operator
- INITIAL.mat_file = ['''../../../iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'''];
-elseif (abs(CO) == 3) %pitch angle operator
- INITIAL.mat_file = ['''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'''];
-elseif (CO == 4) % Full Coulomb GK
-% INITIAL.mat_file = ['''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5'''];
- INITIAL.mat_file = ['''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_75_kpm_6.0.h5'''];
-% INITIAL.mat_file = ['''../../../iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5'''];
-% INITIAL.mat_file = ['''../../../iCa/gk_coulomb_NFLR_6_P_4_J_2_N_75_kpm_6.0.h5'''];
-elseif (CO == -1) % DGDK
- disp('Warning, DGDK not debugged')
-end
% Naming and creating input file
-switch abs(CO)
- case 0; CONAME = 'LB';
- case 1; CONAME = 'DG';
- case 2; CONAME = 'SG';
- case 3; CONAME = 'PA';
- case 4; CONAME = 'FC';
- otherwise; CONAME ='UK';
+CONAME = CO;
+if GKCO
+ CONAME = [CONAME,'GK'];
+else
+ CONAME = [CONAME,'DK'];
end
-if (CO <= 0); CONAME = [CONAME,'DK'];
-else; CONAME = [CONAME,'GK'];
+if ~ABCO
+ CONAME = [CONAME,'aa'];
end
if (CLOS == 0); CLOSNAME = 'Trunc.';
elseif(CLOS == 1); CLOSNAME = 'Clos. 1';
@@ -99,7 +95,7 @@ coll_ = ['_nu_%0.0e_',CONAME];
coll_ = sprintf(coll_,NU);
% nonlinear
lin_ = [];
-if ~NON_LIN; lin_ = '_lin'; end
+if ~LINEARITY; lin_ = '_lin'; end
adiabe_ = [];
if ~KIN_E; adiabe_ = '_adiabe'; end
% resolution and boxsize
@@ -161,7 +157,7 @@ if ~exist(BASIC.MISCDIR, 'dir')
mkdir(BASIC.MISCDIR)
end
%% Compile and WRITE input file
-INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC);
+INPUT = write_fort90(OUTPUTS,GRID,MODEL,COLL,INITIAL,TIME_INTEGRATION,BASIC);
nproc = 1;
MAKE = 'cd ..; make; cd wk';
% system(MAKE);
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 58092ab4055d785053046bb7c0635d86ecbe3dd6..e7901dd410458c0c333331425e12e0bba7b19e8a 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -1,4 +1,4 @@
-function [INPUT] = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC)
+function [INPUT] = write_fort90(OUTPUTS,GRID,MODEL,COLL,INITIAL,TIME_INTEGRATION,BASIC)
% Write the input script "fort.90" with desired parameters
INPUT = ['fort_',sprintf('%2.2d',OUTPUTS.job2load+1),'.90'];
fid = fopen(INPUT,'wt');
@@ -46,10 +46,9 @@ fprintf(fid,'/\n');
fprintf(fid,'&MODEL_PAR\n');
fprintf(fid,' ! Collisionality\n');
-fprintf(fid,[' CO = ', num2str(MODEL.CO),'\n']);
fprintf(fid,[' CLOS = ', num2str(MODEL.CLOS),'\n']);
fprintf(fid,[' NL_CLOS = ', num2str(MODEL.NL_CLOS),'\n']);
-fprintf(fid,[' NON_LIN = ', MODEL.NON_LIN,'\n']);
+fprintf(fid,[' LINEARITY = ', MODEL.LINEARITY,'\n']);
fprintf(fid,[' KIN_E = ', MODEL.KIN_E,'\n']);
fprintf(fid,[' mu = ', num2str(MODEL.mu),'\n']);
fprintf(fid,[' mu_p = ', num2str(MODEL.mu_p),'\n']);
@@ -69,16 +68,23 @@ fprintf(fid,[' CurvB = ', num2str(MODEL.CurvB),'\n']);
fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),'\n']);
fprintf(fid,'/\n');
+fprintf(fid,'&COLLISION_PAR\n');
+fprintf(fid,[' collision_model = ', COLL.collision_model,'\n']);
+fprintf(fid,[' gyrokin_CO = ', COLL.gyrokin_CO,'\n']);
+fprintf(fid,[' interspecies = ', COLL.interspecies,'\n']);
+fprintf(fid,[' mat_file = ', COLL.mat_file,'\n']);
+fprintf(fid,'/\n');
+
+
fprintf(fid,'&INITIAL_CON\n');
fprintf(fid,[' INIT_NOISY_PHI = ', INITIAL.init_noisy_phi,'\n']);
fprintf(fid,[' INIT_ZF = ', num2str(INITIAL.INIT_ZF),'\n']);
-fprintf(fid,[' WIPE_ZF = ', num2str(INITIAL.wipe_zf),'\n']);
+fprintf(fid,[' ACT_ON_MODES = ', INITIAL.ACT_ON_MODES,'\n']);
fprintf(fid,[' WIPE_TURB = ', num2str(INITIAL.wipe_turb),'\n']);
fprintf(fid,[' INIT_BLOB = ', INITIAL.init_blob,'\n']);
fprintf(fid,[' init_background = ', num2str(INITIAL.init_background),'\n']);
fprintf(fid,[' init_noiselvl = ', num2str(INITIAL.init_noiselvl),'\n']);
fprintf(fid,[' iseed = ', num2str(INITIAL.iseed),'\n']);
-fprintf(fid,[' mat_file = ', INITIAL.mat_file,'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&TIME_INTEGRATION_PAR\n');
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index 90def21c619f834eff31f2b698cd70edf23493b4..3978302c2f38791881c0680287c75454758adcec 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -2,13 +2,17 @@ addpath(genpath('../matlab')) % ... add
addpath(genpath('../matlab/plots')) % ... add
outfile ='';
%% Directory of the simulation
-if 0% Local results
+if 1% Local results
outfile ='';
-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';
+outfile ='';
+outfile ='';
+outfile ='';
+outfile ='';
+outfile ='';
+outfile ='';
+outfile ='';
+% outfile ='nonlin_FCGK/nadiab_op_continue';
+% outfile ='Cyclone/150x150x30_5x3_Lx_200_Ly_150_q0_1.4_e_0.18_kN_2.22_kT_6.9_nu_5e-02_SGGK_adiabe';
LOCALDIR = ['../results/',outfile,'/'];
MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
system(['mkdir -p ',MISCDIR]);
@@ -20,14 +24,14 @@ outfile ='';
outfile ='';
outfile ='';
outfile ='';
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x300_5x3_L_120_kN_1.6667_nu_1e-01_SGGK/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';
% BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
end
%% Load the results
% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 13;
+JOBNUMMIN = 00; JOBNUMMAX = 20;
data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
@@ -38,7 +42,7 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-TAVG_0 = 500; TAVG_1 = 1000; % Averaging times duration
+TAVG_0 = 11000; TAVG_1 = 11600; % Averaging times duration
fig = plot_radial_transport_and_shear(data,TAVG_0,TAVG_1);
save_figure(data,fig)
end
@@ -47,13 +51,13 @@ end
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
-options.INTERP = 0;
+options.INTERP = 1;
options.POLARPLOT = 0;
-% options.NAME = '\Gamma_x';
-options.NAME = 'n_i';
+options.NAME = '\phi^{NZ}';
+% options.NAME = 'n_i';
options.PLAN = 'kxky';
options.COMP = 1;
-options.TIME = 450:1:data.Ts3D(end);
+options.TIME = 1500:10:5500;
% options.TIME = 140:0.5:160;
data.a = data.EPS * 2000;
create_film(data,options,'.gif')
@@ -65,10 +69,11 @@ if 0
options.INTERP = 0;
options.POLARPLOT = 0;
% options.NAME = '\Gamma_x';
+% options.NAME = 'n_i^{NZ}|';
options.NAME = 'n_i';
-options.PLAN = 'kxky';
+options.PLAN = 'kxz';
options.COMP = 1;
-options.TIME = [50 100 1200];
+options.TIME = [50 70 200 500];
data.a = data.EPS * 1000;
fig = photomaton(data,options);
save_figure(data,fig)
diff --git a/wk/benchmark_linear_1D_entropy_mode.m b/wk/benchmark_linear_1D_entropy_mode.m
index 86528abce82fb8cbd6474d09e5d158b0d48e3761..f43f102d404acc430925e7825d47bb0cc324bbe5 100644
--- a/wk/benchmark_linear_1D_entropy_mode.m
+++ b/wk/benchmark_linear_1D_entropy_mode.m
@@ -6,7 +6,7 @@ cd ..
system('make');
cd wk
% Run HeLaZ in sequential
-system(['cd ',SIMDIR,';',' ./../../../bin/helaz 0; cd ../../../wk']);
+system(['cd ',SIMDIR,';',' ./../../../bin/helaz3.03 0; cd ../../../wk']);
% Run HeLaZ in parallel (discrepencies can occur at marginal growth rate)
% since the random seed will not be the same)
% system(['cd ',SIMDIR,';',' mpirun -np 6 ./../../../bin/helaz 1 6 0; cd ../../../wk'])
@@ -30,7 +30,7 @@ end
%% Plot
SCALE = 1;%sqrt(2);
-openfig([SIMDIR,'/benchmark_data.fig']); hold on;
+openfig([SIMDIR,'/benchmark_data_nadiab.fig']); hold on;
plot(kx,gamma_phi,'-x','DisplayName','new results');
hold on;
legend('show');
diff --git a/wk/cyclone_test_case.m b/wk/cyclone_test_case.m
index 3e58370e9ef01fce72c13df4d563bd4fed9e579e..0a6b241be4f7cb2c925f5707a68447d591b640fe 100644
--- a/wk/cyclone_test_case.m
+++ b/wk/cyclone_test_case.m
@@ -4,7 +4,7 @@ addpath(genpath('../matlab')) % ... add
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.01; % Collision frequency
+NU = 0.05; % Collision frequency
K_N = 2.22; % Density gradient drive
K_T = 6.9; % Temperature '''
K_E = 0.00; % Electrostat gradient
@@ -12,10 +12,10 @@ SIGMA_E = 0.05196; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
NU_HYP = 0.01;
KIN_E = 0; % Kinetic (1) or adiabatic (2) electron model
%% GRID PARAMETERS
-NX = 100; % Spatial radial resolution ( = 2x radial modes)
+NX = 150; % Spatial radial resolution ( = 2x radial modes)
LX = 200; % Radial window size
-NY = 100; % Spatial azimuthal resolution (= azim modes)
-LY = 100; % Azimuthal window size
+NY = 150; % Spatial azimuthal resolution (= azim modes)
+LY = 150; % Azimuthal window size
NZ = 30; % number of perpendicular planes (parallel grid)
P = 4;
J = 2;
@@ -28,7 +28,7 @@ CURVB = 1.0; % Magnetic curvature
SG = 1; % Staggered z grids option
%% TIME PARAMETERS
TMAX = 500; % Maximal time unit
-DT = 5e-2; % Time step
+DT = 2.5e-2; % Time step
SPS0D = 2; % Sampling per time unit for profiler
SPS2D = 1; % Sampling per time unit for 2D arrays
SPS3D = 1/2; % Sampling per time unit for 3D arrays
@@ -38,17 +38,17 @@ JOB2LOAD= -1;
%% OPTIONS AND NAMING
% Collision operator
% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; 4 : Coulomb; +/- for GK/DK)
-CO = 1;
+CO = 2;
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 = '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; % Non linear model (0: linear, 0.5: semi linear, 1: non linear)
+LINEARITY = 1; % Non linear model (0: linear, 0.5: semi linear, 1: non linear)
% INIT options
INIT_PHI= 0; % Start simulation with a noisy phi (0= noisy moments 00)
INIT_ZF = 0; ZF_AMP = 0.0;
-INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0;
+INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
%% OUTPUTS
W_DOUBLE = 1;
W_GAMMA = 1; W_HF = 1;
diff --git a/wk/linear_1D_damping.m b/wk/linear_1D_damping.m
index 95e0fca00b9fb8792266f52f3a4ed7ba36520e37..35e2a5313cf702586c0cf1280e67fd949c11492b 100644
--- a/wk/linear_1D_damping.m
+++ b/wk/linear_1D_damping.m
@@ -32,7 +32,7 @@ JOB2LOAD= -1;
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)
+LINEARITY = 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)
@@ -55,7 +55,7 @@ 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;
+INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 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;
diff --git a/wk/linear_1D_entropy_mode.m b/wk/linear_1D_entropy_mode.m
index 7f0d2c3ee0db84563dfc3b633db002e1c37997e0..5c23f2f39bcc3b54a50ba30b8d70dc96be6876df 100644
--- a/wk/linear_1D_entropy_mode.m
+++ b/wk/linear_1D_entropy_mode.m
@@ -33,11 +33,13 @@ SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
%% OPTIONS
SIMID = 'Linear_entropy_mode'; % Name of the simulation
-NON_LIN = 0; % activate non-linearity (is cancelled if KXEQ0 = 1)
+LINEARITY = 'linear'; % 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;
+% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+CO = 'LD';
+GKCO = 1; % gyrokinetic operator
+ABCO = 1; % interspecies collisions
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)
@@ -56,7 +58,7 @@ 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;
+INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 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;
@@ -90,7 +92,7 @@ for i = 1:Nparam
system(['rm fort*.90']);
% Run linear simulation
if RUN
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 1 6 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz3 1 6 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ./../../../bin/helaz 1 2 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz 0; cd ../../../wk'])
end
diff --git a/wk/local_run.m b/wk/local_run.m
index 3e58370e9ef01fce72c13df4d563bd4fed9e579e..c5d950fe21d2ed88de46e7f4262a6d8b123874ea 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -4,51 +4,49 @@ addpath(genpath('../matlab')) % ... add
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.01; % Collision frequency
-K_N = 2.22; % Density gradient drive
-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.01;
-KIN_E = 0; % Kinetic (1) or adiabatic (2) electron model
+NU = 0.1; % Collision frequency
+K_N = 1/0.6; % Density gradient drive
+K_T = 0.0; % Temperature '''
+K_E = 0.0; % Electrostat gradient
+SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+NU_HYP = 0.0;
+KIN_E = 1; % Kinetic (1) or adiabatic (2) electron model
%% GRID PARAMETERS
-NX = 100; % Spatial radial resolution ( = 2x radial modes)
-LX = 200; % Radial window size
+NX = 200; % Spatial radial resolution ( = 2x radial modes)
+LX = 120; % Radial window size
NY = 100; % Spatial azimuthal resolution (= azim modes)
-LY = 100; % Azimuthal window size
-NZ = 30; % number of perpendicular planes (parallel grid)
+LY = 120; % Azimuthal window size
+NZ = 1; % number of perpendicular planes (parallel grid)
P = 4;
J = 2;
%% GEOMETRY PARAMETERS
-Q0 = 1.4; % safety factor
+Q0 = 1.0; % safety factor
SHEAR = 0.0; % magnetic shear
-EPS = 0.18; % inverse aspect ratio
+EPS = 0.0; % inverse aspect ratio
GRADB = 1.0; % Magnetic gradient
CURVB = 1.0; % Magnetic curvature
-SG = 1; % Staggered z grids option
+SG = 0; % Staggered z grids option
%% TIME PARAMETERS
-TMAX = 500; % Maximal time unit
-DT = 5e-2; % Time step
-SPS0D = 2; % Sampling per time unit for profiler
+TMAX = 1000; % Maximal time unit
+DT = 1e-2; % Time step
+SPS0D = 1; % Sampling per time unit for profiler
SPS2D = 1; % Sampling per time unit for 2D 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;
+JOB2LOAD= 1;
%% OPTIONS AND NAMING
% Collision operator
% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; 4 : Coulomb; +/- for GK/DK)
-CO = 1;
+CO = 2;
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 = '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; % Non linear model (0: linear, 0.5: semi linear, 1: non linear)
+SIMID = 'semi_linear'; % Name of the simulation
+LINEARITY = 'nonlinear'; % (nonlinear, semilinear, linear)
% INIT options
-INIT_PHI= 0; % Start simulation with a noisy phi (0= noisy moments 00)
+INIT_PHI = 0; % Start simulation with a noisy phi (0= noisy moments 00)
INIT_ZF = 0; ZF_AMP = 0.0;
-INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0;
+INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 'donothing';
%% OUTPUTS
W_DOUBLE = 1;
W_GAMMA = 1; W_HF = 1;
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index e6ea6a5e084628d2ca50a3740e9cd8bc329b0927..f9836ed232d0504c7f18113acc5009525ff34986 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -2,8 +2,8 @@ clear all;
addpath(genpath('../matlab')) % ... add
SUBMIT = 1; % To submit the job automatically
CHAIN = 2; % To chain jobs (CHAIN = n will launch n jobs in chain)
-% EXECNAME = 'helaz_dbg';
- EXECNAME = 'helaz_3.5';
+% EXECNAME = 'helaz3_dbg';
+ EXECNAME = 'helaz3';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -50,10 +50,10 @@ NL_CLOS = 0; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=
% SIMID = 'test_chained_job'; % 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)
+LINEARITY = 1; % activate non-linearity (is cancelled if KXEQ0 = 1)
% INIT options
INIT_ZF = 0; ZF_AMP = 0.0;
-INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0;
+INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
%% OUTPUTS
W_DOUBLE = 1;
W_GAMMA = 1; W_HF = 1;
diff --git a/wk/mode_growth_meter.m b/wk/mode_growth_meter.m
new file mode 100644
index 0000000000000000000000000000000000000000..e00a08aab3eb88eda7a56665088929357f3993bb
--- /dev/null
+++ b/wk/mode_growth_meter.m
@@ -0,0 +1,87 @@
+MODES_SELECTOR = 1; %(1:Zonal, 2: NZonal, 3: ky=kx)
+NORMALIZED = 1;
+options.K2PLOT = 0.01:0.05:1.0;
+options.TIME =500:0.5:1000;
+DATA = data;
+OPTIONS = options;
+
+t = OPTIONS.TIME;
+iz = 1;
+
+if MODES_SELECTOR == 1
+ if NORMALIZED
+ plt = @(x,ik) movmean(abs(squeeze(x(ik,1,iz,FRAMES)))./max(abs(squeeze(x(ik,1,iz,FRAMES)))),1);
+ else
+ plt = @(x,ik) movmean(abs(squeeze(x(ik,1,iz,FRAMES))),1);
+ end
+ kstr = 'k_x';
+ k = DATA.kx;
+ MODESTR = 'Zonal modes';
+elseif MODES_SELECTOR == 2
+ if NORMALIZED
+ plt = @(x,ik) movmean(abs(squeeze(x(1,ik,iz,FRAMES)))./max(abs(squeeze(x(1,ik,iz,FRAMES)))),1);
+ else
+ plt = @(x,ik) movmean(abs(squeeze(x(1,ik,iz,FRAMES))),1);
+ end
+ kstr = 'k_y';
+ k = DATA.ky;
+ MODESTR = 'NZ modes';
+elseif MODES_SELECTOR == 3
+ if NORMALIZED
+ plt = @(x,ik) movmean(abs(squeeze(x(ik,ik,iz,FRAMES)))./max(abs(squeeze(x(ik,ik,iz,FRAMES)))),1);
+ else
+ plt = @(x,ik) movmean(abs(squeeze(x(ik,ik,iz,FRAMES))),1);
+ end
+ kstr = 'k_y=k_x';
+ k = DATA.ky;
+ MODESTR = 'Diag modes';
+end
+
+
+FRAMES = zeros(size(OPTIONS.TIME));
+for i = 1:numel(OPTIONS.TIME)
+ [~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts3D));
+end
+
+MODES = zeros(size(OPTIONS.K2PLOT));
+for i = 1:numel(OPTIONS.K2PLOT)
+ [~,MODES(i)] =min(abs(OPTIONS.K2PLOT(i)-k));
+end
+
+
+% plt = @(x,ik) abs(squeeze(x(1,ik,iz,FRAMES)))./max(abs(squeeze(x(1,ik,iz,FRAMES))));
+
+gamma = MODES;
+amp = MODES;
+i_=1;
+for ik = MODES
+ gr = polyfit(t',log(plt(DATA.PHI,ik)),1);
+ gamma(i_) = gr(1);
+ amp(i_) = gr(2);
+ i_=i_+1;
+end
+
+%plot
+figure; set(gcf, 'Position', [100 100 1200 350])
+subplot(1,3,1)
+ [YY,XX] = meshgrid(t,fftshift(k,1));
+ pclr = pcolor(XX,YY,abs(plt(fftshift(DATA.PHI,MODES_SELECTOR),1:numel(k))));set(pclr, 'edgecolor','none'); hold on;
+% xlim([t(1) t(end)]); %ylim([1e-5 1])
+ xlabel(['$',kstr,'\rho_s$']); ylabel('$t c_s /\rho_s$');
+ title(MODESTR)
+
+subplot(1,3,2)
+ mod2plot = round(linspace(2,numel(MODES),15));
+ for i_ = mod2plot
+ semilogy(t,plt(DATA.PHI,MODES(i_))); hold on;
+ semilogy(t,exp(gamma(i_).*t+amp(i_)),'--k')
+ end
+ xlim([t(1) t(end)]); %ylim([1e-5 1])
+ xlabel('$t c_s /\rho_s$'); ylabel(['$|\phi_{',kstr,'}|$']);
+ title('Measure time window')
+
+subplot(1,3,3)
+ plot(k(MODES),gamma); hold on;
+ plot(k(MODES(mod2plot)),gamma(mod2plot),'x')
+ xlabel(['$',kstr,'\rho_s$']); ylabel('$\gamma$');
+ title('Growth rates')
\ No newline at end of file
diff --git a/wk/shearless_linear_fluxtube.m b/wk/shearless_linear_fluxtube.m
index a86a52477f31b8baf1512e72b0e17dbf11162dbb..51ed3e45e7f3be15c5c4483dea4296e7ccfb38dc 100644
--- a/wk/shearless_linear_fluxtube.m
+++ b/wk/shearless_linear_fluxtube.m
@@ -58,12 +58,12 @@ MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
K_E = 0.00; % Electrostat '''
GRADB = 1.0;
CURVB = 1.0;
-INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0; INIT_PHI= 0;
+INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0; INIT_PHI= 0;
NOISE0 = 0.0e-4; % Init noise amplitude
BCKGD0 = 1.0e-4; % Init background
LAMBDAD = 0.0;
KXEQ0 = 0; % put kx = 0
-NON_LIN = 0; % activate non-linearity (is cancelled if KXEQ0 = 1)
+LINEARITY = 0; % activate non-linearity (is cancelled if KXEQ0 = 1)
%% PARAMETER SCANS
%% Parameter scan over PJ