Skip to content
Snippets Groups Projects
Commit 8236a826 authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann :seedling:
Browse files

save latest versions of matlab scripts

parent 11c2edfa
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
matlab/compile_results.m
+ 1
− 1
View file @ 8236a826
......@@ -253,7 +253,7 @@ else
KPERP2 = KX.^2+KY.^2;
%% Add everything in output structure
% scaling
DATA.scale = (1/Nx/Ny)^2;
DATA.scale = 1;%(1/Nx/Ny)^2;
% Fields
DATA.GGAMMA_RI = GGAMMAI_; DATA.PGAMMA_RI = PGAMMAI_; DATA.HFLUX_X = HFLUXI_;
DATA.Nipj = Nipj_; DATA.Ni00 = Ni00_; DATA.DENS_I = DENS_I_; DATA.TEMP_I = TEMP_I_;
......
matlab/compute/compute_polarisation.m
+ 1
− 1
View file @ 8236a826
......@@ -14,7 +14,7 @@ KP = sqrt(KX.^2+KY.^2);
GAMMA2_ = 0.*KN_;
for in = 1:Ji+1
for in = 1:1
for iz = 1:data.Nz
GAMMA2_(:,:,iz) = GAMMA2_(:,:,iz) + kernel(in-1,KP*sqrt(2)).^2;
end
......
matlab/evolve_tracers.m
+ 8
− 8
View file @ 8236a826
% Options
SHOW_FILM = 0;
U_TIME = 50; % >0 for frozen velocity at a given time, -1 for evolving field
Evolve_U = 1; % 0 for frozen velocity at a given time, 1 for evolving field
Tfin = 50;
dt_ = 0.05;
Nstep = ceil(Tfin/dt_);
INIT = 'lin'; % lin (for a line)/ round (for a small round)/ gauss for random
U_TIME = 700; % >0 for frozen velocity at a given time, -1 for evolving field
Evolve_U = 1; % 0 for frozen velocity at a given time, 1 for evolving field
Tfin = 10;
dt_ = 0.0005;
Nstep = ceil(Tfin/dt_);
% Init tracers
Np = 100; %number of tracers
Np = 200; %number of tracers
% color = tcolors;
color = jet(Np);
tcolors = distinguishable_colors(Np); %Their colors
Na = 10000; %length of trace
Na = 1000000; %length of trace
Traj_x = zeros(Np,Nstep);
Traj_y = zeros(Np,Nstep);
......@@ -21,7 +22,6 @@ Disp_y = zeros(Np,Nstep);
xmax = max(data.x); xmin = min(data.x);
ymax = max(data.y); ymin = min(data.y);
INIT = 'round';
switch INIT
case 'lin'
% Evenly distributed initial positions
......
matlab/plot/spectrum_1D.m
+ 31
− 11
View file @ 8236a826
......@@ -5,6 +5,14 @@ options.COMP = 'avg';
options.INTERP = 0;
options.POLARPLOT = 0;
options.AXISEQUAL = 1;
switch options.COMPT
case 'avg'
[~,it0] = min(abs(data.Ts3D-options.TIME(1)));
[~,it1] = min(abs(data.Ts3D-options.TIME(end)));
options.TIME = data.Ts3D(it0:it1);
end
toplot = process_field(data,options);
t = data.Ts3D; frames = toplot.FRAMES;
......@@ -13,29 +21,40 @@ colors = jet(numel(frames));
switch options.NAME
case '\Gamma_x'
FIGURE.fig = figure; FIGURE.FIGNAME = ['transp_spectrum_',data.PARAMS];
yname = '$\sum_{k_y}|\Gamma_k|$'; fieldname = 'transport';
fieldname = 'transport';
case '\phi'
FIGURE.fig = figure; FIGURE.FIGNAME = ['phi_spectrum_',data.PARAMS];
yname = '$\sum_{k_y}|\phi|$'; fieldname = 'ES pot.';
fieldname = 'ES pot.';
case 'n_i'
FIGURE.fig = figure; FIGURE.FIGNAME = ['ni_spectrum_',data.PARAMS];
yname = '$\sum_{k_y}|n_i|$'; fieldname = 'ion dens.';
yname = 'n_i'; fieldname = 'ion dens.';
case 'n_e'
FIGURE.fig = figure; FIGURE.FIGNAME = ['ne_spectrum_',data.PARAMS];
yname = '$\sum_{k_y}|n_e|$'; fieldname = 'elec. dens.';
fieldname = 'elec. dens.';
end
PLOT2D = 0;
switch options.COMPXY
case 'avg'
compx = @(x) mean(x,2);
compy = @(x) mean(x,1);
compx = @(x) mean(x,2);
compy = @(x) mean(x,1);
ynamecx = ['$\sum_{k_x}',options.NAME,'$'];
ynamecy = ['$\sum_{k_y}',options.NAME,'$'];
case 'sum'
compx = @(x) sum(x,2);
compy = @(x) sum(x,1);
ynamecx = ['$\sum_{k_x}',options.NAME,'$'];
ynamecy = ['$\sum_{k_y}',options.NAME,'$'];
case 'max'
compx = @(x) max(x,2);
compy = @(x) max(x,1);
ynamecx = ['$\max_{k_x}',options.NAME,'$'];
ynamecy = ['$\max_{k_y}',options.NAME,'$'];
case 'zero'
compx = @(x) x(:,data.Nx/2);
compy = @(x) x(1,:);
ynamecx = ['$',options.NAME,'(k_x=',num2str(data.kx(1)),')$'];
ynamecy = ['$',options.NAME,'(k_y=',num2str(data.ky(1)),')$'];
otherwise
compx = @(x) x(:,:);
compy = @(x) x(:,:);
......@@ -85,8 +104,9 @@ if ~PLOT2D
end
end
grid on
title(['HeLaZ $k_x$ ',fieldname,' spectrum']); legend('show','Location','eastoutside')
xlabel(xname); ylabel(yname)
title(['GM $k_x$ ',fieldname]);
% legend('show','Location','eastoutside')
xlabel(xname); ylabel(ynamecx)
subplot(1,2,2)
......@@ -124,9 +144,9 @@ if ~PLOT2D
end
end
grid on
% title('HeLaZ $k_y$ transport spectrum');
legend('show','Location','eastoutside');
xlabel(xname); ylabel(yname)
title(['GM $k_y$ ',fieldname]);
% legend('show','Location','eastoutside');
xlabel(xname); ylabel(ynamecy)
else
% it0 = toplot.FRAMES(1); it1 = toplot.FRAMES(end);
Gk = mean(abs(toplot.FIELD(:,:,:)),3);
......
wk/analysis_HeLaZ.m
+ 14
− 14
View file @ 8236a826
......@@ -11,7 +11,7 @@ system(['mkdir -p ',LOCALDIR]);
CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
system(CMD);
% Load outputs from jobnummin up to jobnummax
JOBNUMMIN = 00; JOBNUMMAX = 00;
JOBNUMMIN = 00; JOBNUMMAX = 20;
data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
data.localdir = LOCALDIR;
data.FIGDIR = LOCALDIR;
......@@ -23,7 +23,7 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
options.TAVG_0 = 0.7*data.Ts3D(end); data.scale = (1/data.Nx/data.Ny)^2;
options.TAVG_0 = 0.7*data.Ts3D(end); data.scale = 1;%(1/data.Nx/data.Ny)^2;
options.TAVG_1 = data.Ts3D(end); % Averaging times duration
options.NMVA = 1; % Moving average for time traces
% options.ST_FIELD = '\Gamma_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
......@@ -45,18 +45,18 @@ if 0
options.INTERP = 1;
options.POLARPLOT = 0;
% options.NAME = '\phi';
options.NAME = 'N_i^{00}';
% options.NAME = 'v_y';
% options.NAME = 'N_i^{00}';
% options.NAME = 'v_x';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_e';
options.PLAN = 'kxky';
options.NAME = 'n_i';
options.PLAN = 'xy';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 1;
% options.TIME = data.Ts5D(end-30:end);
options.TIME = data.Ts3D(50:2:end);
% options.TIME = [350:600];
% options.TIME = data.Ts3D(20:2:end);
options.TIME = [700:1:750];
data.EPS = 0.1;
data.a = data.EPS * 2000;
create_film(data,options,'.gif')
......@@ -69,7 +69,7 @@ options.INTERP = 1;
options.POLARPLOT = 0;
options.AXISEQUAL = 1;
% options.NAME = '\phi';
options.NAME = 'n_e';
% options.NAME = 'n_e';
% options.NAME = 'N_i^{00}';
% options.NAME = 'T_i';
% options.NAME = '\Gamma_x';
......@@ -78,7 +78,7 @@ options.PLAN = 'xy';
% options.NAME 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
options.TIME = [600];
options.TIME = [200 800 1500];
data.a = data.EPS * 2e3;
fig = photomaton(data,options);
% save_figure(data,fig)
......@@ -131,16 +131,16 @@ end
if 0
%% Time averaged spectrum
options.TIME = 1100:10:1400;
options.TIME = [500 700];
options.NORM =1;
% options.NAME = '\phi';
% options.NAME = 'N_i^{00}';
options.NAME ='\Gamma_x';
options.NAME ='\Gamma_x';
options.PLAN = 'kxky';
options.COMPZ = 'avg';
options.OK = 0;
options.COMPXY = 'avg';
options.COMPT = 0;
options.COMPXY = 'avg'; % avg/sum/max/zero/ 2D plot otherwise
options.COMPT = 'avg';
options.PLOT = 'semilogy';
fig = spectrum_1D(data,options);
% save_figure(data,fig,'.png')
......
wk/analysis_gene.m
+ 12
− 12
View file @ 8236a826
......@@ -10,8 +10,8 @@
% folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/';
% folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/';
% folder = '/misc/gene_results/LD_zpinch_1.6/';
folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
% folder = '/misc/gene_results/ZP_HP_kn_2.5/';
% folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
folder = '/misc/gene_results/ZP_kn_2.5_large_box/';
gene_data = load_gene_data(folder);
gene_data = invert_kxky_to_kykx_gene_results(gene_data);
if 1
......@@ -39,15 +39,15 @@ options.INTERP = 1;
options.POLARPLOT = 0;
options.AXISEQUAL = 1;
% options.NAME = 'Q_x';
options.NAME = '\phi';
% options.NAME = 'n_i';
% options.NAME = '\phi';
options.NAME = 'n_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
options.PLAN = 'xy';
% options.NAME ='f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
options.TIME = [600];
options.TIME = [1:10];
gene_data.a = data.EPS * 2000;
fig = photomaton(gene_data,options);
save_figure(gene_data,fig,'.png')
......@@ -58,12 +58,12 @@ if 0
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
options.NAME = '\phi';
% options.NAME = '\phi';
% options.NAME = 'v_y';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
options.PLAN = 'xy';
options.NAME = 'n_i';
options.PLAN = 'kxky';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
......@@ -115,16 +115,16 @@ if 0
%% Time averaged spectrum
options.TIME = 300:600;
options.NORM =1;
options.NAME = '\phi';
% options.NAME = '\phi';
% options.NAME = 'n_i';
% options.NAME ='\Gamma_x';
options.NAME ='\Gamma_x';
options.PLAN = 'kxky';
options.COMPZ = 'avg';
options.OK = 0;
options.COMPXY = 'avg';
options.COMPXY = 'zero';
options.COMPT = 'avg';
options.PLOT = 'semilogy';
fig = spectrum_1D(data,options);
fig = spectrum_1D(gene_data,options);
% save_figure(data,fig)
end
......
wk/header_3D_results.m
+ 4
− 1
View file @ 8236a826
......@@ -28,7 +28,10 @@ helazdir = '/home/ahoffman/HeLaZ/';
% outfile = 'shearless_cyclone/64x32x16x5x3_CBC_CO/64x32x16x5x3_CBC_LRGK';
%% ZPINCH
% outfile ='Zpinch_rerun/Kn_2.5_200x48x5x3';
outfile ='Zpinch_rerun/Kn_2.0_200x48x9x5_large_box';
% outfile ='Zpinch_rerun/Kn_2.5_256x128x5x3_new_NL';
outfile ='Zpinch_rerun/Kn_2.5_256x64x5x3';
% outfile ='Zpinch_rerun/Kn_2.0_200x48x9x5_large_box';
% outfile ='Zpinch_rerun/Kn_2.0_256x64x9x5_Lx_240_Ly_120';
% outfile ='Zpinch_rerun/Kn_1.6_256x128x7x4';
% outfile ='Zpinch_rerun/Kn_1.6_200x48x11x6';
% outfile ='Zpinch_rerun/Kn_1.6_200x48x21x11';
......
wk/quick_run.m
+ 22
− 22
View file @ 8236a826
......@@ -13,48 +13,48 @@ EXECNAME = 'helaz3';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
NU = 0.0; % Collision frequency
NU = 0.5; % Collision frequency
TAU = 1.0; % e/i temperature ratio
K_N = 1.0;%2.0; % Density gradient drive
K_N = 6.0;%2.0; % Density gradient drive
K_T = 0;%0.25*K_N; % Temperature '''
K_E = 0.0; % Electrostat '''
% SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
%% GRID PARAMETERS
PMAXE = 30; % Hermite basis size of electrons
JMAXE = 15; % Laguerre "
PMAXI = 30; % " ions
JMAXI = 15; % "
NX = 1; % real space x-gridpoints
NY = 32; % '' y-gridpoints
PMAXE = 4; % Hermite basis size of electrons
JMAXE = 2; % Laguerre "
PMAXI = 4; % " ions
JMAXI = 2; % "
NX = 32; % real space x-gridpoints
NY = 1; % '' y-gridpoints
LX = 100; % Size of the squared frequency domain
LY = 30; % Size of the squared frequency domain
NZ = 1; % number of perpendicular planes (parallel grid)
LY = 2*pi/0.5; % Size of the squared frequency domain
NZ = 16; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
%% GEOMETRY
% GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
GEOMETRY= 's-alpha';
Q0 = 1.0; % safety factor
SHEAR = 0.0; % magnetic shear (Not implemented yet)
EPS = 0.0; % inverse aspect ratio
Q0 = 1.4; % safety factor
SHEAR = 0.796; % magnetic shear (Not implemented yet)
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
TMAX = 500; % Maximal time unit
DT = 5e-3; % Time step
TMAX = 50; % Maximal time unit
DT = 1e-3; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
SPS3D = 1; % Sampling per time unit for 2D arrays
SPS5D = 1/2; % Sampling per time unit for 5D arrays
SPS5D = 1/5; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
%% OPTIONS
SIMID = 'dbg'; % Name of the simulation
SIMID = 'shear_testcase_Pan22_linear'; % Name of the simulation
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
CO = 'DG';
GKCO = 1; % gyrokinetic operator
CO = 'SG';
GKCO = 0; % 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))s
......@@ -77,7 +77,7 @@ INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
MU_Y = MU; %
N_HD = 4;
MU_Z = 0.0; %
MU_Z = 0.05; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
......@@ -92,8 +92,8 @@ system(['rm fort*.90']);
% Run linear simulation
if RUN
% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 1 4 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',HELAZDIR,'bin/',EXECNAME,' 1 1 2 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
end
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment