Skip to content
Snippets Groups Projects
Commit c259b712 authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann
Browse files

scripts update

parent ce192c00
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 240 additions and 128 deletions
Makefile
+ 2
− 2
View file @ c259b712
......@@ -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
......
matlab/compile_results.m
+ 3
− 2
View file @ c259b712
......@@ -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
......
matlab/load_params.m
+ 9
− 18
View file @ c259b712
......@@ -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
matlab/photomaton.m
+ 4
− 4
View file @ c259b712
......@@ -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)]);
......
matlab/process_field.m
+ 32
− 6
View file @ c259b712
......@@ -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));
......
matlab/setup.m
+ 24
− 28
View file @ c259b712
......@@ -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);
......
matlab/write_fort90.m
+ 11
− 5
View file @ c259b712
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');
......
wk/analysis_3D.m
+ 20
− 15
View file @ c259b712
......@@ -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)
......
wk/benchmark_linear_1D_entropy_mode.m
+ 2
− 2
View file @ c259b712
......@@ -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');
......
wk/cyclone_test_case.m
+ 8
− 8
View file @ c259b712
......@@ -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;
......
wk/linear_1D_damping.m
+ 2
− 2
View file @ c259b712
......@@ -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;
......
wk/linear_1D_entropy_mode.m
+ 7
− 5
View file @ c259b712
......@@ -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
......
wk/local_run.m
+ 23
− 25
View file @ c259b712
......@@ -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;
......
wk/marconi_run.m
+ 4
− 4
View file @ c259b712
......@@ -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;
......
wk/mode_growth_meter.m 0 → 100644
+ 87
− 0
View file @ c259b712
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
wk/shearless_linear_fluxtube.m
+ 2
− 2
View file @ c259b712
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment