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

save scripts

parent fb72c0f3
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 1287 additions and 112 deletions
matlab/load/quick_gene_load.m
+ 9
− 4
View file @ f93045fb
...@@ -76,11 +76,16 @@ path = '/home/ahoffman/gene/linear_CBC_results/'; ...@@ -76,11 +76,16 @@ path = '/home/ahoffman/gene/linear_CBC_results/';
% fname = 'CBC_salpha_nz_24_nv_scan_nw_16_kine.txt'; % fname = 'CBC_salpha_nz_24_nv_scan_nw_16_kine.txt';
% fname = 'CBC_miller_nz_24_nv_scan_nw_16_kine.txt'; % fname = 'CBC_miller_nz_24_nv_scan_nw_16_kine.txt';
%---------- CBC %---------- CBC
% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
% fname = 'CBC_salpha_nx_8_nz_18_nv_18_nw_8_adiabe.txt';
% fname = 'CBC_salpha_nx_8_nz_18_nv_12_nw_8_adiabe.txt'; % fname = 'CBC_salpha_nx_8_nz_18_nv_12_nw_8_adiabe.txt';
% fname = 'CBC_salpha_nx_8_nz_24_nv_60_nw_30_adiabe.txt'; % fname = 'CBC_salpha_nx_8_nz_18_nv_18_nw_8_adiabe.txt';
fname = 'kT_5.3_salpha_nx_8_nz_24_nv_60_nw_30_adiabe.txt';
% fname = 'CBC_salpha_nx_8_nz_24_nv_8_nw_4_adiabe.txt';
% fname = 'CBC_salpha_nx_8_nz_24_nv_16_nw_8_adiabe.txt';
% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
fname = 'CBC_salpha_nx_8_nz_24_nv_48_nw_16_adiabe.txt';
% fname = 'kT_5.3_salpha_nx_8_nz_24_nv_60_nw_30_adiabe.txt';
% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_kine.txt'; % fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_kine.txt';
% fname = 'CBC_miller_nx_20_nz_32_nv_32_nw_12_adiabe.txt'; % fname = 'CBC_miller_nx_20_nz_32_nv_32_nw_12_adiabe.txt';
% fname = 'CBC_miller_nx_8_nz_24_nv_36_nw_16_adiabe.txt'; % fname = 'CBC_miller_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
......
matlab/plot/create_film.m
+ 2
− 2
View file @ f93045fb
...@@ -36,7 +36,7 @@ if hmax == hmin ...@@ -36,7 +36,7 @@ if hmax == hmin
disp('Warning : h = hmin = hmax = const') disp('Warning : h = hmin = hmax = const')
else else
% Setup figure frame % Setup figure frame
fig = figure('Color','white','Position', toplot.DIMENSIONS.*[0 0 1.2 1]); fig = figure('Color','white');%,'Position', toplot.DIMENSIONS.*[0.5 0.5 1.0 1]);
if ~strcmp(OPTIONS.PLAN,'sx') if ~strcmp(OPTIONS.PLAN,'sx')
pcolor(X,Y,FIELD(:,:,1)); % to set up pcolor(X,Y,FIELD(:,:,1)); % to set up
if BWR if BWR
...@@ -70,7 +70,7 @@ fig = figure('Color','white','Position', toplot.DIMENSIONS.*[0 0 1.2 1]); ...@@ -70,7 +70,7 @@ fig = figure('Color','white','Position', toplot.DIMENSIONS.*[0 0 1.2 1]);
if toplot.INTERP if toplot.INTERP
shading interp; shading interp;
end end
set(pclr, 'edgecolor','none'); pbaspect(toplot.ASPECT); set(pclr, 'edgecolor','none'); %pbaspect(toplot.ASPECT);
if BWR if BWR
colormap(bluewhitered) colormap(bluewhitered)
end end
......
matlab/plot/imagesc_custom.m
+ 2
− 2
View file @ f93045fb
...@@ -5,8 +5,8 @@ function [ fig ] = imagesc_custom( XX,YY,FF) ...@@ -5,8 +5,8 @@ function [ fig ] = imagesc_custom( XX,YY,FF)
fig = imagesc( XX(1,:), YY(:,1), FF); fig = imagesc( XX(1,:), YY(:,1), FF);
set(gca,'YDir','normal') set(gca,'YDir','normal')
xticks(XX(1,:)); xticks(XX(1,:));
% xticklabels(num2str(XX(1,:))); xticklabels(num2str(XX(1,:)));
yticks(YY(:,1)); yticks(YY(:,1));
% yticklabels(num2str(YY(:,1))); yticklabels(num2str(YY(:,1)));
end end
matlab/plot/show_moments_spectrum.m
+ 44
− 22
View file @ f93045fb
...@@ -70,32 +70,54 @@ for ia = 1:DATA.inputs.Na ...@@ -70,32 +70,54 @@ for ia = 1:DATA.inputs.Na
end end
% plots % plots
% scaling % scaling
if OPTIONS.NORMALIZED if OPTIONS.TAVG_2D
plt = @(x,i) reshape(x(i,:)./max(x(i,:)),numel(p2j),numel(Time_)); nt_half = ceil(numel(Time_)/2);
Napjz_tavg = mean(Napjz(:,:,nt_half:end),3);
x_ = DATA.grids.Parray;
y_ = DATA.grids.Jarray;
if OPTIONS.TAVG_2D_CTR
[YY,XX] = meshgrid(y_,x_);
surf(XX,YY,Napjz_tavg);
else
imagesc(x_,y_,Napjz_tavg');
set(gca,'YDir','normal')
xlabel('$p$'); ylabel('$j$');
end
clb = colorbar; colormap(jet);
clb.Label.String = '$\log\langle E(p,j) \rangle_t$';
clb.TickLabelInterpreter = 'latex';
clb.Label.Interpreter = 'latex';
clb.Label.FontSize= 18;
else else
plt = @(x,i) reshape(x(i,:),numel(p2j),numel(Time_)); if OPTIONS.NORMALIZED
end plt = @(x,i) reshape(x(i,:)./max(x(i,:)),numel(p2j),numel(Time_));
if OPTIONS.ST else
imagesc(Time_,p2j,plt(Na_ST,1:numel(p2j))); plt = @(x,i) reshape(x(i,:),numel(p2j),numel(Time_));
set(gca,'YDir','normal')
xlabel('$t$'); ylabel(yname);
if ~isempty(ticks_labels)
yticks(p2j);
yticklabels(ticks_labels)
end end
if OPTIONS.FILTER Nlabelmax = 15;
caxis([0 0.2]); nskip = floor(numel(p2j)/Nlabelmax);
title('Relative fluctuation from the average'); if OPTIONS.ST
imagesc(Time_,p2j,plt(Na_ST,1:numel(p2j)));
set(gca,'YDir','normal')
xlabel('$t$'); ylabel(yname);
if ~isempty(ticks_labels)
yticks(p2j(1:nskip:end));
yticklabels(ticks_labels(1:nskip:end))
end end
colorbar if OPTIONS.FILTER
else caxis([0 0.2]);
colors_ = jet(numel(p2j)); title('Relative fluctuation from the average');
for i = 1:numel(p2j) end
semilogy(Time_,squeeze(Na_ST(i,:)),... colorbar
'DisplayName',ticks_labels{i},... else
'color',colors_(i,:)); hold on; colors_ = jet(numel(p2j));
for i = 1:numel(p2j)
semilogy(Time_,squeeze(Na_ST(i,:)),...
'DisplayName',ticks_labels{i},...
'color',colors_(i,:)); hold on;
end
title(plotname)
end end
title(plotname)
end end
top_title(DATA.paramshort) top_title(DATA.paramshort)
......
wk/CBC_P_J_scan.m 0 → 100644
+ 224
− 0
View file @ f93045fb
gyacomodir = pwd;
gyacomodir = gyacomodir(1:end-2);
addpath(genpath([gyacomodir,'matlab'])) % ... add
addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
EXECNAME = 'gyacomo23_sp';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%
SIMID = 'p2_linear_new'; % Name of the simulation
RERUN = 0; % rerun if the data does not exist
RUN = 0;
NU = 1e-3;
K_T= 6.96;
P_a = 2:2:40;
J_a = 2:2:40;
% collision setting
CO = 'DG';
GKCO = 0; % gyrokinetic operator
COLL_KCUT = 1.75;
% model
KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 1e-4; % electron plasma beta
% background gradients setting
K_N = 2.22; % Density '''
% Geometry
% GEOMETRY= 'miller';
GEOMETRY= 's-alpha';
SHEAR = 0.8; % magnetic shear
% time and numerical grid
DT0 = 1e-2;
TMAX = 30;
kymin = 0.3;
NY = 2;
% arrays for the result
g_ky = zeros(numel(P_a),numel(J_a),NY/2+1);
g_avg= g_ky*0;
g_std= g_ky*0;
% Naming of the collision operator
if GKCO
CONAME = [CO,'GK'];
else
CONAME = [CO,'DK'];
end
j = 1;
for P = P_a
i = 1;
for J = J_a
%% PHYSICAL PARAMETERS
TAU = 1.0; % e/i temperature ratio
% 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)
K_Te = K_T; % ele Temperature '''
K_Ti = K_T; % ion Temperature '''
K_Ne = K_N; % ele Density '''
K_Ni = K_N; % ion Density gradient drive
%% GRID PARAMETERS
DT = DT0/sqrt(J);
PMAX = P; % Hermite basis size
JMAX = J; % Laguerre "
NX = 8; % real space x-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
LY = 2*pi/kymin; % Size of the squared frequency domain
NZ = 24; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
%% GEOMETRY
% GEOMETRY= 's-alpha';
EPS = 0.18; % inverse aspect ratio
Q0 = 1.4; % safety factor
KAPPA = 1.0; % elongation
DELTA = 0.0; % triangularity
ZETA = 0.0; % squareness
PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
SHIFT_Y = 0.0;
%% TIME PARMETERS
DTSAVE0D = 1; % Sampling per time unit for 0D arrays
DTSAVE2D = -1; % Sampling per time unit for 2D arrays
DTSAVE3D = 2; % Sampling per time unit for 3D arrays
DTSAVE5D = 100; % Sampling per time unit for 5D arrays
JOB2LOAD = -1; % Start a new simulation serie
%% OPTIONS
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator
ABCO = 1; % INTERSPECIES collisions
INIT_ZF = 0; ZF_AMP = 0.0;
CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
KERN = 0; % Kernel model (0 : GK)
INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom
NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
%% OUTPUTS
W_DOUBLE = 0;
W_GAMMA = 1; W_HF = 1;
W_PHI = 1; W_NA00 = 1;
W_DENS = 0; W_TEMP = 1;
W_NAPJ = 0; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
NA = 1;
ADIAB_E = (NA==1);
HD_CO = 0.0; % Hyper diffusivity cutoff ratio
MU = 0.0; % Hyperdiffusivity coefficient
INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
MU_Y = MU; %
N_HD = 4;
HYP_V = 'none';
HRCY_CLOS = 'truncation'; % Closure model for higher order moments
DMAX = -1;
NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
NMAX = 0;
MU_Z = 1.0; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
k_gB = 1.0;
k_cB = 1.0;
%% RUN
setup
% naming
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
% check if data exist to run if no data
data_ = {};
try
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
Ntime = numel(data_.Ts0D);
catch
data_.outfilenames = [];
end
if RUN && (RERUN || isempty(data_.outfilenames) || Ntime < 10)
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
end
try
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
[data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
if numel(data_.Ts3D)>10
if numel(data_.Ts3D)>10
% Load results after trying to run
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
[data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
% linear growth rate (adapted for 2D zpinch and fluxtube)
options.TRANGE = [0.5 1]*data_.Ts3D(end);
options.NPLOTS = 0; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3
[~,it1] = min(abs(data_.Ts3D-0.5*data_.Ts3D(end))); % start of the measurement time window
[~,it2] = min(abs(data_.Ts3D-1.0*data_.Ts3D(end))); % end of ...
field = 0;
field_t = 0;
for ik = 2:NY/2+1
field = squeeze(sum(abs(data_.PHI),3)); % take the sum over z
field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only
to_measure = log(field_t(it1:it2));
tw = double(data_.Ts3D(it1:it2));
% gr = polyfit(tw,to_measure,1);
gr = fit(tw,to_measure,'poly1');
err= confint(gr);
g_ky(i,j,ik) = gr.p1;
g_std(i,j,ik) = abs(err(2,1)-err(1,1))/2;
end
[gmax, ikmax] = max(g_ky(i,j,:));
msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
end
end
catch
g_ky(i,j,:) = 0;
g_std(i,j,:) = 0;
end
i = i + 1;
end
j = j + 1;
end
if 0
%% Check time evolution
figure;
to_measure = log(field_t);
plot(data_.Ts3D,to_measure); hold on
plot(data_.Ts3D(it1:it2),to_measure(it1:it2),'--');
end
%% take max growth rate among z coordinate
y_ = g_ky(:,:,2);
e_ = g_std(:,:,2);
%%
if(numel(J_a)>1 && numel(P_a)>1)
%% Save metadata
pmin = num2str(min(P_a)); pmax = num2str(max(P_a));
jmin = num2str(min(J_a)); jmax = num2str(max(J_a));
filename = [num2str(NX),'x',num2str(NZ),'_ky_',num2str(kymin),...
'_P_',pmin,'_',pmax,...
'_J_',jmin,'_',jmax,'_',CONAME,'_',num2str(NU),'_kT_',num2str(K_T),'.mat'];
metadata.name = filename;
metadata.kymin = kymin;
metadata.title = ['$\nu_{',CONAME,'}=$',num2str(NU),'$\kappa_T=$',num2str(K_T),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)];
metadata.par = [num2str(NX),'x1x',num2str(NZ)];
metadata.nscan = 2;
metadata.s2name = '$P$';
metadata.s2 = P_a;
metadata.s1name = '$J$';
metadata.s1 = J_a;
metadata.dname = '$\gamma c_s/R$';
metadata.data = y_;
metadata.err = e_;
save([SIMDIR,filename],'-struct','metadata');
disp(['saved in ',SIMDIR,filename]);
clear metadata tosave
end
wk/CBC_kT_PJ_scan.m
+ 14
− 11
View file @ f93045fb
...@@ -10,20 +10,20 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss ...@@ -10,20 +10,20 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
SIMID = 'p2_linear_new'; % Name of the simulation SIMID = 'p2_linear_new'; % Name of the simulation
RERUN = 0; % rerun if the data does not exist RERUN = 0; % rerun if the data does not exist
RUN = 1; RUN = 1;
KT_a = [3:0.5:6.5 6.96]; % KT_a = [3:0.5:6.5 6.96];
% KT_a = [5.5 6]; KT_a = [6.96];
% P_a = [25]; P_a = [2];
% P_a = [3:1:29]; % P_a = [3:1:29];
P_a = 2:2:10; % P_a = 2:2:10;
J_a = floor(P_a/2); J_a = floor(P_a/2);
% collision setting % collision setting
CO = 'LD'; CO = 'DG';
NU = 0.1; NU = 1e-3;
GKCO = 1; % gyrokinetic operator GKCO = 0; % gyrokinetic operator
COLL_KCUT = 1.75; COLL_KCUT = 1.75;
% model % model
KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 1e-4; % electron plasma beta BETA = 0; % electron plasma beta
% background gradients setting % background gradients setting
K_N = 2.22; % Density ''' K_N = 2.22; % Density '''
% Geometry % Geometry
...@@ -31,7 +31,7 @@ K_N = 2.22; % Density ''' ...@@ -31,7 +31,7 @@ K_N = 2.22; % Density '''
GEOMETRY= 's-alpha'; GEOMETRY= 's-alpha';
SHEAR = 0.8; % magnetic shear SHEAR = 0.8; % magnetic shear
% time and numerical grid % time and numerical grid
DT0 = 1e-2; DT0 = 5e-3;
TMAX = 30; TMAX = 30;
kymin = 0.3; kymin = 0.3;
NY = 2; NY = 2;
...@@ -134,13 +134,16 @@ for KT = KT_a ...@@ -134,13 +134,16 @@ for KT = KT_a
data_ = {}; data_ = {};
try try
data_ = compile_results_low_mem(data_,LOCALDIR,00,00); data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
if numel(data_.Ts3D) < 10
data_.outfilenames = [];
end
catch catch
data_.outfilenames = []; data_.outfilenames = [];
end end
if RUN && (RERUN || isempty(data_.outfilenames)) if RUN && (RERUN || isempty(data_.outfilenames))
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk']) system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
end end
data_ = compile_results_low_mem(data_,LOCALDIR,00,00); data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
[data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi'); [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
......
wk/CBC_kT_nu_scan.m
+ 12
− 11
View file @ f93045fb
...@@ -9,16 +9,17 @@ EXECNAME = 'gyacomo23_sp'; ...@@ -9,16 +9,17 @@ EXECNAME = 'gyacomo23_sp';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%% %%
SIMID = 'p2_linear_new'; % Name of the simulation SIMID = 'p2_linear_new'; % Name of the simulation
RERUN = 1; % rerun if the data does not exist RERUN = 0; % rerun if the data does not exist
RUN = 1; RUN = 1;
KT_a = [3.5 4.0 4.5 5.0 5.5 6.0 6.5 6.96]; KT_a = [2.5 3.0 3.5 4.0 4.5 5.0];
NU_a = [0 0.001 0.002 0.005 0.01 0.02 0.05 0.1 0.2 0.5]; % KT_a = [3.0 3.5 4.0 4.5];
% KT_a = 3.5; NU_a = [1e-3 2e-3 5e-3 1e-2 2e-2 5e-2 1e-1 2e-1 5e-1 1e+0];
% NU_a = 0.5; % KT_a = 4.5;
P = 8; % NU_a = 0.01;
J = 4; P = 16;
J = P/2;
% collision setting % collision setting
CO = 'LD'; CO = 'DG';
GKCO = 1; % gyrokinetic operator GKCO = 1; % gyrokinetic operator
COLL_KCUT = 1.75; COLL_KCUT = 1.75;
% model % model
...@@ -31,7 +32,7 @@ K_N = 2.22; % Density ''' ...@@ -31,7 +32,7 @@ K_N = 2.22; % Density '''
GEOMETRY= 's-alpha'; GEOMETRY= 's-alpha';
SHEAR = 0.8; % magnetic shear SHEAR = 0.8; % magnetic shear
% time and numerical grid % time and numerical grid
DT = 1e-2; DT = 5e-3;
TMAX = 30; TMAX = 30;
kymin = 0.3; kymin = 0.3;
NY = 2; NY = 2;
...@@ -136,8 +137,8 @@ for KT = KT_a ...@@ -136,8 +137,8 @@ for KT = KT_a
data_.outfilenames = []; data_.outfilenames = [];
end end
if RUN && (RERUN || isempty(data_.outfilenames)) if RUN && (RERUN || isempty(data_.outfilenames))
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk']) system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
end end
data_ = compile_results_low_mem(data_,LOCALDIR,00,00); data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
......
wk/CBC_ky_PJ_scan.m 0 → 100644
+ 222
− 0
View file @ f93045fb
gyacomodir = pwd;
gyacomodir = gyacomodir(1:end-2);
addpath(genpath([gyacomodir,'matlab'])) % ... add
addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
EXECNAME = 'gyacomo23_sp';
% EXECNAME = 'gyacomo23_dp';
% EXECNAME = 'gyacomo23_debug';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%
SIMID = 'p2_linear_new'; % Name of the simulation
RERUN = 0; % rerun if the data does not exist
RUN = 1;
K_T = 5.3;
P_a = [2 4 8 10];
% P_a = 10;
ky_a= [0.05:0.05:1.0];
% ky_a = 0.05;
% collision setting
CO = 'LD';
GKCO = 1; % gyrokinetic operator
COLL_KCUT = 1.75;
NU = 1e-2;
% model
KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 1e-4; % electron plasma beta
% background gradients setting
K_N = 2.22; % Density '''
% Geometry
% GEOMETRY= 'miller';
GEOMETRY= 's-alpha';
SHEAR = 0.8; % magnetic shear
% time and numerical grid
DT0 = 5e-3;
TMAX = 50;
% arrays for the result
g_ky = zeros(numel(ky_a),numel(P_a),2);
g_avg= g_ky*0;
g_std= g_ky*0;
% Naming of the collision operator
if GKCO
CONAME = [CO,'GK'];
else
CONAME = [CO,'DK'];
end
j = 1;
for P = P_a
J = P/2;
i = 1;
for ky = ky_a
%% PHYSICAL PARAMETERS
TAU = 1.0; % e/i temperature ratio
% 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)
K_Te = K_T; % ele Temperature '''
K_Ti = K_T; % ion Temperature '''
K_Ne = K_N; % ele Density '''
K_Ni = K_N; % ion Density gradient drive
%% GRID PARAMETERS
DT = DT0/sqrt(J);
PMAX = P; % Hermite basis size
JMAX = P/2; % Laguerre "
NX = 8; % real space x-gridpoints
NY = 2;
LX = 2*pi/0.8; % Size of the squared frequency domain
LY = 2*pi/ky; % Size of the squared frequency domain
NZ = 24; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
%% GEOMETRY
% GEOMETRY= 's-alpha';
EPS = 0.18; % inverse aspect ratio
Q0 = 1.4; % safety factor
KAPPA = 1.0; % elongation
DELTA = 0.0; % triangularity
ZETA = 0.0; % squareness
PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
SHIFT_Y = 0.0;
%% TIME PARMETERS
DTSAVE0D = 1; % Sampling per time unit for 0D arrays
DTSAVE2D = -1; % Sampling per time unit for 2D arrays
DTSAVE3D = 2; % Sampling per time unit for 3D arrays
DTSAVE5D = 100; % Sampling per time unit for 5D arrays
JOB2LOAD = -1; % Start a new simulation serie
%% OPTIONS
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator
ABCO = 1; % INTERSPECIES collisions
INIT_ZF = 0; ZF_AMP = 0.0;
CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
KERN = 0; % Kernel model (0 : GK)
INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom
NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
%% OUTPUTS
W_DOUBLE = 0;
W_GAMMA = 1; W_HF = 1;
W_PHI = 1; W_NA00 = 1;
W_DENS = 0; W_TEMP = 1;
W_NAPJ = 0; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
NA = 1;
ADIAB_E = (NA==1);
HD_CO = 0.0; % Hyper diffusivity cutoff ratio
MU = 0.0; % Hyperdiffusivity coefficient
INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
MU_Y = MU; %
N_HD = 4;
HYP_V = 'none';
HRCY_CLOS = 'truncation'; % Closure model for higher order moments
DMAX = -1;
NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
NMAX = 0;
MU_Z = 1.0; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
k_gB = 1.0;
k_cB = 1.0;
%% RUN
setup
% naming
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
% check if data exist to run if no data
data_ = {};
try
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
Ntime = numel(data_.Ts0D);
catch
data_.outfilenames = [];
end
if RUN && (RERUN || isempty(data_.outfilenames) || Ntime < 10)
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 2 2 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
end
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
[data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
if numel(data_.Ts3D)>10
if numel(data_.Ts3D)>10
% Load results after trying to run
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
[data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
% linear growth rate (adapted for 2D zpinch and fluxtube)
options.TRANGE = [0.5 1]*data_.Ts3D(end);
options.NPLOTS = 0; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3
[~,it1] = min(abs(data_.Ts3D-0.5*data_.Ts3D(end))); % start of the measurement time window
[~,it2] = min(abs(data_.Ts3D-1.0*data_.Ts3D(end))); % end of ...
field = 0;
field_t = 0;
for ik = 2:NY/2+1
field = squeeze(sum(abs(data_.PHI),3)); % take the sum over z
field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only
to_measure = log(field_t(it1:it2));
tw = double(data_.Ts3D(it1:it2));
% gr = polyfit(tw,to_measure,1);
gr = fit(tw,to_measure,'poly1');
err= confint(gr);
g_ky(i,j,ik) = gr.p1;
g_std(i,j,ik) = abs(err(2,1)-err(1,1))/2;
end
[gmax, ikmax] = max(g_ky(i,j,:));
msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
end
end
i = i + 1;
end
j = j + 1;
end
if 0
%% Check time evolution
figure;
to_measure = log(field_t);
plot(data_.Ts3D,to_measure); hold on
plot(data_.Ts3D(it1:it2),to_measure(it1:it2),'--');
end
%% take max growth rate among z coordinate
y_ = g_ky(:,:,2);
e_ = g_std(:,:,2);
%%
if(numel(ky_a)>1 && numel(P_a)>1)
%% Save metadata
pmin = num2str(min(P_a)); pmax = num2str(max(P_a));
kymin = num2str(min(ky_a)); kymax= num2str(max(ky_a));
filename = [num2str(NX),'x',num2str(NZ),'_ky_',kymin,'_',kymax,...
'_P_',pmin,'_',pmax,'_',CONAME,'_',num2str(NU),'_kT_',num2str(K_T),'.mat'];
metadata.name = filename;
metadata.kymin = ky;
metadata.title = ['$\nu_{',CONAME,'}=$',num2str(NU),'$\kappa_T=$',num2str(K_T),', $\kappa_N=$',num2str(K_N)];
metadata.par = [num2str(NX),'x1x',num2str(NZ)];
metadata.nscan = 2;
metadata.s2name = '$P$';
metadata.s2 = P_a;
metadata.s1name = '$ky$';
metadata.s1 = ky_a;
metadata.dname = '$\gamma c_s/R$';
metadata.data = y_;
metadata.err = e_;
save([SIMDIR,filename],'-struct','metadata');
disp(['saved in ',SIMDIR,filename]);
clear metadata tosave
end
wk/CBC_ky_PJ_scan_colless.m 0 → 100644
+ 221
− 0
View file @ f93045fb
gyacomodir = pwd;
gyacomodir = gyacomodir(1:end-2);
addpath(genpath([gyacomodir,'matlab'])) % ... add
addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
% EXECNAME = 'gyacomo23_sp';
EXECNAME = 'gyacomo23_dp';
% EXECNAME = 'gyacomo23_debug';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%
SIMID = 'p2_linear_new'; % Name of the simulation
RERUN = 0; % rerun if the data does not exist
RUN = 1;
K_T = 6.96;
% P_a = [2 4 8 16 32 48];
P_a = 64;
% ky_a= [0.05:0.05:1.0];
ky_a= [0.25:0.05:0.65];
% collision setting
CO = 'DG';
GKCO = 0; % gyrokinetic operator
COLL_KCUT = 1.75;
NU = 1e-3;
% model
KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 1e-4; % electron plasma beta
% background gradients setting
K_N = 2.22; % Density '''
% Geometry
% GEOMETRY= 'miller';
GEOMETRY= 's-alpha';
SHEAR = 0.8; % magnetic shear
% time and numerical grid
DT0 = 2e-3;
TMAX = 30;
% arrays for the result
g_ky = zeros(numel(ky_a),numel(P_a),NY/2+1);
g_avg= g_ky*0;
g_std= g_ky*0;
% Naming of the collision operator
if GKCO
CONAME = [CO,'GK'];
else
CONAME = [CO,'DK'];
end
j = 1;
for P = P_a
i = 1;
for ky = ky_a
%% PHYSICAL PARAMETERS
TAU = 1.0; % e/i temperature ratio
% 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)
K_Te = K_T; % ele Temperature '''
K_Ti = K_T; % ion Temperature '''
K_Ne = K_N; % ele Density '''
K_Ni = K_N; % ion Density gradient drive
%% GRID PARAMETERS
DT = DT0/sqrt(J);
PMAX = P; % Hermite basis size
JMAX = P/2; % Laguerre "
NX = 8; % real space x-gridpoints
NY = 2;
LX = 2*pi/0.8; % Size of the squared frequency domain
LY = 2*pi/ky; % Size of the squared frequency domain
NZ = 24; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
%% GEOMETRY
% GEOMETRY= 's-alpha';
EPS = 0.18; % inverse aspect ratio
Q0 = 1.4; % safety factor
KAPPA = 1.0; % elongation
DELTA = 0.0; % triangularity
ZETA = 0.0; % squareness
PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
SHIFT_Y = 0.0;
%% TIME PARMETERS
DTSAVE0D = 1; % Sampling per time unit for 0D arrays
DTSAVE2D = -1; % Sampling per time unit for 2D arrays
DTSAVE3D = 2; % Sampling per time unit for 3D arrays
DTSAVE5D = 100; % Sampling per time unit for 5D arrays
JOB2LOAD = -1; % Start a new simulation serie
%% OPTIONS
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator
ABCO = 1; % INTERSPECIES collisions
INIT_ZF = 0; ZF_AMP = 0.0;
CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
KERN = 0; % Kernel model (0 : GK)
INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom
NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
%% OUTPUTS
W_DOUBLE = 0;
W_GAMMA = 1; W_HF = 1;
W_PHI = 1; W_NA00 = 1;
W_DENS = 0; W_TEMP = 1;
W_NAPJ = 0; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
NA = 1;
ADIAB_E = (NA==1);
HD_CO = 0.0; % Hyper diffusivity cutoff ratio
MU = 0.0; % Hyperdiffusivity coefficient
INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
MU_Y = MU; %
N_HD = 4;
HYP_V = 'none';
HRCY_CLOS = 'truncation'; % Closure model for higher order moments
DMAX = -1;
NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
NMAX = 0;
MU_Z = 1.0; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
k_gB = 1.0;
k_cB = 1.0;
%% RUN
setup
% naming
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
% check if data exist to run if no data
data_ = {};
try
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
Ntime = numel(data_.Ts0D);
catch
data_.outfilenames = [];
end
if RUN && (RERUN || isempty(data_.outfilenames) || Ntime < 10)
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 2 2 1 0; cd ../../../wk'])
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
end
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
[data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
if numel(data_.Ts3D)>10
if numel(data_.Ts3D)>10
% Load results after trying to run
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
[data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
% linear growth rate (adapted for 2D zpinch and fluxtube)
options.TRANGE = [0.5 1]*data_.Ts3D(end);
options.NPLOTS = 0; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3
[~,it1] = min(abs(data_.Ts3D-0.5*data_.Ts3D(end))); % start of the measurement time window
[~,it2] = min(abs(data_.Ts3D-1.0*data_.Ts3D(end))); % end of ...
field = 0;
field_t = 0;
for ik = 2:NY/2+1
field = squeeze(sum(abs(data_.PHI),3)); % take the sum over z
field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only
to_measure = log(field_t(it1:it2));
tw = double(data_.Ts3D(it1:it2));
% gr = polyfit(tw,to_measure,1);
gr = fit(tw,to_measure,'poly1');
err= confint(gr);
g_ky(i,j,ik) = gr.p1;
g_std(i,j,ik) = abs(err(2,1)-err(1,1))/2;
end
[gmax, ikmax] = max(g_ky(i,j,:));
msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
end
end
i = i + 1;
end
j = j + 1;
end
if 0
%% Check time evolution
figure;
to_measure = log(field_t);
plot(data_.Ts3D,to_measure); hold on
plot(data_.Ts3D(it1:it2),to_measure(it1:it2),'--');
end
%% take max growth rate among z coordinate
y_ = g_ky(:,:,2);
e_ = g_std(:,:,2);
%%
if(numel(ky_a)>1 && numel(P_a)>1)
%% Save metadata
pmin = num2str(min(P_a)); pmax = num2str(max(P_a));
kymin = num2str(min(ky_a)); kymax= num2str(max(ky_a));
filename = [num2str(NX),'x',num2str(NZ),'_ky_',kymin,'_',kymax,...
'_P_',pmin,'_',pmax,'_',CONAME,'_',num2str(NU),'_kT_',num2str(K_T),'.mat'];
metadata.name = filename;
metadata.kymin = ky;
metadata.title = ['$\nu_{',CONAME,'}=$',num2str(NU),'$\kappa_T=$',num2str(K_T),', $\kappa_N=$',num2str(K_N)];
metadata.par = [num2str(NX),'x1x',num2str(NZ)];
metadata.nscan = 2;
metadata.s2name = '$P$';
metadata.s2 = P_a;
metadata.s1name = '$ky$';
metadata.s1 = ky_a;
metadata.dname = '$\gamma c_s/R$';
metadata.data = y_;
metadata.err = e_;
save([SIMDIR,filename],'-struct','metadata');
disp(['saved in ',SIMDIR,filename]);
clear metadata tosave
end
wk/CBC_nu_PJ_scan.m 0 → 100644
+ 223
− 0
View file @ f93045fb
gyacomodir = pwd;
gyacomodir = gyacomodir(1:end-2);
addpath(genpath([gyacomodir,'matlab'])) % ... add
addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
EXECNAME = 'gyacomo23_sp';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%
SIMID = 'p2_linear_new'; % Name of the simulation
RERUN = 0; % rerun if the data does not exist
RUN = 1;
nu_a = [1e-3,2e-3,5e-3,1e-2,2e-2,5e-2,1e-1];
% KT_a = [5.5 6];
% P_a = [25];
% P_a = [3:1:29];
P_a = 2:2:4;
J_a = floor(P_a/2);
% collision setting
CO = 'DG';
K_T = 6.96;
GKCO = 1; % gyrokinetic operator
COLL_KCUT = 1.75;
% model
KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 1e-4; % electron plasma beta
% background gradients setting
K_N = 2.22; % Density '''
% Geometry
% GEOMETRY= 'miller';
GEOMETRY= 's-alpha';
SHEAR = 0.8; % magnetic shear
% time and numerical grid
DT0 = 1e-2;
TMAX = 30;
kymin = 0.3;
NY = 2;
% arrays for the result
g_ky = zeros(numel(nu_a),numel(P_a),NY/2+1);
g_avg= g_ky*0;
g_std= g_ky*0;
% Naming of the collision operator
if GKCO
CONAME = [CO,'GK'];
else
CONAME = [CO,'DK'];
end
j = 1;
for P = P_a
i = 1;
for NU = nu_a
%% PHYSICAL PARAMETERS
TAU = 1.0; % e/i temperature ratio
% 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)
K_Te = K_T; % ele Temperature '''
K_Ti = K_T; % ion Temperature '''
K_Ne = K_N; % ele Density '''
K_Ni = K_N; % ion Density gradient drive
%% GRID PARAMETERS
J = floor(P/2);
DT = DT0/sqrt(J);
PMAX = P; % Hermite basis size
JMAX = J; % Laguerre "
NX = 8; % real space x-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
LY = 2*pi/kymin; % Size of the squared frequency domain
NZ = 24; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
%% GEOMETRY
% GEOMETRY= 's-alpha';
EPS = 0.18; % inverse aspect ratio
Q0 = 1.4; % safety factor
KAPPA = 1.0; % elongation
DELTA = 0.0; % triangularity
ZETA = 0.0; % squareness
PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
SHIFT_Y = 0.0;
%% TIME PARMETERS
DTSAVE0D = 1; % Sampling per time unit for 0D arrays
DTSAVE2D = -1; % Sampling per time unit for 2D arrays
DTSAVE3D = 2; % Sampling per time unit for 3D arrays
DTSAVE5D = 100; % Sampling per time unit for 5D arrays
JOB2LOAD = -1; % Start a new simulation serie
%% OPTIONS
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator
ABCO = 1; % INTERSPECIES collisions
INIT_ZF = 0; ZF_AMP = 0.0;
CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
KERN = 0; % Kernel model (0 : GK)
INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom
NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
%% OUTPUTS
W_DOUBLE = 0;
W_GAMMA = 1; W_HF = 1;
W_PHI = 1; W_NA00 = 1;
W_DENS = 0; W_TEMP = 1;
W_NAPJ = 0; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
ADIAB_E = (NA==1);
HD_CO = 0.0; % Hyper diffusivity cutoff ratio
MU = 0.0; % Hyperdiffusivity coefficient
INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
MU_Y = MU; %
N_HD = 4;
HYP_V = 'none';
HRCY_CLOS = 'truncation'; % Closure model for higher order moments
DMAX = -1;
NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
NMAX = 0;
MU_Z = 1.0; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
k_gB = 1.0;
k_cB = 1.0;
%% RUN
setup
% naming
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
% check if data exist to run if no data
data_ = {};
try
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
catch
data_.outfilenames = [];
end
if RUN && (RERUN || isempty(data_.outfilenames))
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
end
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
[data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
if numel(data_.Ts3D)>10
if numel(data_.Ts3D)>10
% Load results after trying to run
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
[data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
% linear growth rate (adapted for 2D zpinch and fluxtube)
options.TRANGE = [0.5 1]*data_.Ts3D(end);
options.NPLOTS = 0; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3
[~,it1] = min(abs(data_.Ts3D-0.5*data_.Ts3D(end))); % start of the measurement time window
[~,it2] = min(abs(data_.Ts3D-1.0*data_.Ts3D(end))); % end of ...
field = 0;
field_t = 0;
for ik = 2:NY/2+1
field = squeeze(sum(abs(data_.PHI),3)); % take the sum over z
field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only
to_measure = log(field_t(it1:it2));
tw = double(data_.Ts3D(it1:it2));
% gr = polyfit(tw,to_measure,1);
gr = fit(tw,to_measure,'poly1');
err= confint(gr);
g_ky(i,j,ik) = gr.p1;
g_std(i,j,ik) = abs(err(2,1)-err(1,1))/2;
end
[gmax, ikmax] = max(g_ky(i,j,:));
msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
end
end
i = i + 1;
end
j = j + 1;
end
if 0
%% Check time evolution
figure;
to_measure = log(field_t);
plot(data_.Ts3D,to_measure); hold on
plot(data_.Ts3D(it1:it2),to_measure(it1:it2),'--');
end
%% take max growth rate among z coordinate
y_ = g_ky(:,:,2);
e_ = g_std(:,:,2);
%%
if(numel(nu_a)>1 && numel(P_a)>1)
%% Save metadata
numin = num2str(min(nu_a)); numax = num2str(max(nu_a));
pmin = num2str(min(P_a)); pmax = num2str(max(P_a));
filename = [num2str(NX),'x',num2str(NZ),'_ky_',num2str(kymin),...
'_nu_',numin,'_',numax,...
'_P_',pmin,'_',pmax,'_',CONAME,'_kT_',num2str(K_T),'.mat'];
metadata.name = filename;
metadata.kymin = kymin;
metadata.title = ['$\kappa_T=$',num2str(K_T),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)];
metadata.par = [num2str(NX),'x1x',num2str(NZ)];
metadata.nscan = 2;
metadata.s1name = ['$\nu_{',CONAME,'}=$'];
metadata.s1 = nu_a;
metadata.s2name = '$P$, $J=\lfloor P/2 \rfloor$';
metadata.s2 = P_a;
metadata.dname = '$\gamma c_s/R$';
metadata.data = y_;
metadata.err = e_;
metadata.date = date;
save([SIMDIR,filename],'-struct','metadata');
disp(['saved in ',SIMDIR,filename]);
clear metadata tosave
end
wk/analysis_gene.m
+ 3
− 2
View file @ f93045fb
...@@ -53,13 +53,14 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add ...@@ -53,13 +53,14 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add
% folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x32x16_Nexc_5_00/'; % folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x32x16_Nexc_5_00/';
% folder = '/misc/gene_results/CBC/KT_4.5_192x96x24x30x16_00/'; % folder = '/misc/gene_results/CBC/KT_4.5_192x96x24x30x16_00/';
folder = '/misc/gene_results/CBC/KT_5.0_192x96x24x30x16_00/'; % folder = '/misc/gene_results/CBC/KT_5.0_192x96x24x30x16_00/';
% folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x16x8_Nexc_5/'; % folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x16x8_Nexc_5/';
% folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x8x4_Nexc_5/'; % folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x8x4_Nexc_5/';
% folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x8x4_Nexc_5_smallvbox/'; % folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x8x4_Nexc_5_smallvbox/';
% folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x16x8_Nexc_5_largexbox/'; % folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x16x8_Nexc_5_largexbox/';
% folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x16x8_Muz_0.02/'; % folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x16x8_Muz_0.02/';
% folder = '/misc/gene_results/kT_scan_nu0/30x16x128x64x24/kT_4.0/';
folder = '/misc/gene_results/kT_scan_nu0/42x24x128x64x24/kT_6.5/';
% Miller CBC % Miller CBC
% folder = '/misc/gene_results/CBC/Miller_KT_6.96_128x64x24x16x8/'; % folder = '/misc/gene_results/CBC/Miller_KT_6.96_128x64x24x16x8/';
% folder = '/misc/gene_results/CBC/Miller_KT_4.1_128x64x24x16x8/'; % folder = '/misc/gene_results/CBC/Miller_KT_4.1_128x64x24x16x8/';
......
wk/benchmark_and_scan_scripts/Dimits_2000_fig3_dashed_low.txt 0 → 100644
+ 24
− 0
View file @ f93045fb
6.074498567335244, 0.0664424234873966
6.275071633237822, 0.6469830987085263
6.50429799426934, 1.1655306198722162
6.762177650429798, 1.6634928544732954
7.048710601719197, 2.2443894712488195
7.335243553008594, 2.7631742867821103
7.679083094555873, 3.4064199991694704
8.13753581661891, 3.9880284990537884
8.538681948424069, 4.569399704568509
9.083094555873924, 5.254883814744286
9.570200573065902, 5.712387359328931
9.999999999999998, 6.148949675796567
10.716332378223493, 6.752329934091489
11.20343839541547, 7.16842561118131
11.518624641833808, 7.376770067688219
12.177650429799423, 7.814281561634245
13.15186246418338, 8.377321777097535
13.810888252148997, 8.752721469801324
14.84240687679083, 9.191775377149739
15.702005730659023, 9.54730166639971
16.446991404011456, 9.798833698173425
17.335243553008596, 10.133774700860785
18.28080229226361, 10.406841196675504
19.16905444126074, 10.617558596878393
wk/benchmark_and_scan_scripts/Dimits_2000_fig3_full_no_GF.txt
+ 11
− 21
View file @ f93045fb
5.297267537770027 -0.051044540400214444 5.3, 0.03534795644626776
4.754512410510305 0.004436559014473929 5.5, 0.4494240176739783
5.089866441088001 0.2660247587319642 6.0, 0.0578349376676659
5.451261144141954 0.43432466047664775 6.0, 1.2154805112829408
5.709658557756645 0.4345657778143899 6.5, 1.568171058860658
6.046499478583758 0.005642145703188106 6.96, 2.0654884014517894
5.2693381128147845 0.9193804088144457 6.96, 0.9906106990689576
5.914809225953066 1.1625954673959136 6.96, 2.4583399084740396
6.948639997749572 1.0515850450989888 9.0, 4.489742780495503
6.327280618385599 1.6108808217278856 9.0, 5.027378885908156
6.791552052209941 2.003226953703461 9.0, 2.216032823102413
6.894710086541032 2.0966358103451412
6.9714657723892515 2.4512953024314665
9.039247874651133 2.1732870120136702
9.036515211490045 3.4423357840231628
8.905066076197096 4.487314214067991
8.955620344677234 5.00991193189239
12.023959026127072 8.073428266587367
13.934050389504964 9.026999113893783
16.026145156655943 9.458189248979773
17.988719727215916 10.038558680927418
wk/benchmark_and_scan_scripts/Lin_1999_fig2.txt
+ 1
− 1
View file @ f93045fb
0, 0.022113821138211365 0.001, 0.022113821138211365
0.04541040549724383, 0.17300813008130078 0.04541040549724383, 0.17300813008130078
0.08992977422034283, 0.2809756097560976 0.08992977422034283, 0.2809756097560976
0.1803843539983388, 0.45788617886178873 0.1803843539983388, 0.45788617886178873
......
wk/fast_analysis.m
+ 32
− 8
View file @ f93045fb
...@@ -48,10 +48,19 @@ PARTITION = '/misc/gyacomo23_outputs/'; ...@@ -48,10 +48,19 @@ PARTITION = '/misc/gyacomo23_outputs/';
% resdir = 'paper_2_GYAC23/collision_study/nuLDGK_scan_CBC/9x2x64x48x16/nu_0.1'; % resdir = 'paper_2_GYAC23/collision_study/nuLDGK_scan_CBC/9x2x64x48x16/nu_0.1';
% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/9x5x128x64x24/nu_0.5'; % resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/9x5x128x64x24/nu_0.5';
% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/5x3x128x64x24/nu_0.5'; % resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/5x3x128x64x24/nu_0.5';
resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/9x5x128x64x24/nu_0.5'; % resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/9x5x128x64x24/nu_0.5';
% resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/9x5x128x64x24_Lx200/nu_0.5'; % resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/9x5x128x64x24_Lx200/nu_0.5';
% resdir = 'paper_2_GYAC23/collision_study/nuLDGK_scan_kT_5.3/9x2x128x64x24/nu_0.01'; % resdir = 'paper_2_GYAC23/collision_study/nuLDGK_scan_kT_5.3/9x5x128x64x24/nu_0.005';
% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/11x6x192x96x24/nu_0.1';
% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/9x5x192x96x24/nu_0.2';
% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/17x9x192x96x24/nu_0.2';
% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/9x5x192x96x24/nu_0.02';
resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/17x9x192x96x24/nu_0.02';
% resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/9x5x192x96x24/nu_0.2';
% resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/17x9x192x96x24/nu_0.02';
%% kT eff study %% kT eff study
% resdir = 'paper_2_GYAC23/kT_eff_study/1x3x128x64x24_kT_3.0/Lx120'; % resdir = 'paper_2_GYAC23/kT_eff_study/1x3x128x64x24_kT_3.0/Lx120';
% resdir = 'paper_2_GYAC23/kT_eff_study/1x3x128x64x24_kT_3.0/Lx240'; % resdir = 'paper_2_GYAC23/kT_eff_study/1x3x128x64x24_kT_3.0/Lx240';
...@@ -61,6 +70,16 @@ resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/9x5x128x64x24/nu_0.5 ...@@ -61,6 +70,16 @@ resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/9x5x128x64x24/nu_0.5
% resdir = 'paper_2_GYAC23/kT_eff_study/7x3x128x64x24_kT_3.6/truncation'; % resdir = 'paper_2_GYAC23/kT_eff_study/7x3x128x64x24_kT_3.6/truncation';
% resdir = 'paper_2_GYAC23/kT_eff_study/7x3x128x64x24_kT_3.6/dmax'; % resdir = 'paper_2_GYAC23/kT_eff_study/7x3x128x64x24_kT_3.6/dmax';
%% dev
% PARTITION='';
% resdir = '/home/ahoffman/gyacomo/testcases/zpinch_3D';
% resdir = '/home/ahoffman/gyacomo';
%% CBC benchmark
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/3x2x128x64x24/kT_7.0';
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/5x3x128x64x24/kT_7.0';
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/9x5x128x64x24/kT_7.0';
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/17x9x128x64x24/kT_7.0';
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/31x16x128x64x24/kT_7.0';
%% %%
J0 = 00; J1 = 10; J0 = 00; J1 = 10;
...@@ -87,13 +106,15 @@ end ...@@ -87,13 +106,15 @@ end
if 0 if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options % Options
[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi'); % [data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
% [data.Na00, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Na00');
% data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
options.INTERP = 1; options.INTERP = 1;
options.POLARPLOT = 0; options.POLARPLOT = 0;
options.NAME = '\phi'; % options.NAME = '\phi';
% options.NAME = '\phi^{NZ}'; % options.NAME = '\phi^{NZ}';
% options.NAME = '\omega_z'; % options.NAME = '\omega_z';
% options.NAME = 'N_i^{00}'; options.NAME = 'N_i^{00}';
% options.NAME = 's_{Ey}'; % options.NAME = 's_{Ey}';
% options.NAME = 'n_i^{NZ}'; % options.NAME = 'n_i^{NZ}';
% options.NAME = 'Q_x'; % options.NAME = 'Q_x';
...@@ -104,7 +125,7 @@ options.PLAN = 'xy'; ...@@ -104,7 +125,7 @@ options.PLAN = 'xy';
% options.PLAN = 'sx'; % options.PLAN = 'sx';
options.COMP = 'avg'; options.COMP = 'avg';
% options.TIME = data.Ts5D(end-30:end); % options.TIME = data.Ts5D(end-30:end);
options.TIME = data.Ts3D; options.TIME = data.Ts3D(1:1:end);
% options.TIME = [0:1500]; % options.TIME = [0:1500];
data.EPS = 0.1; data.EPS = 0.1;
data.a = data.EPS * 2000; data.a = data.EPS * 2000;
...@@ -116,6 +137,7 @@ if 0 ...@@ -116,6 +137,7 @@ if 0
%% fields snapshots %% fields snapshots
% Options % Options
[data.Na00, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Na00'); [data.Na00, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Na00');
[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
options.INTERP = 0; options.INTERP = 0;
...@@ -126,7 +148,7 @@ options.NAME = 'N_i^{00}'; ...@@ -126,7 +148,7 @@ options.NAME = 'N_i^{00}';
% options.NAME = '\phi'; % options.NAME = '\phi';
options.PLAN = 'xy'; options.PLAN = 'xy';
options.COMP = 'avg'; options.COMP = 'avg';
options.TIME = [10 50 80]; options.TIME = [100 200 300];
options.RESOLUTION = 256; options.RESOLUTION = 256;
fig = photomaton(data,options); fig = photomaton(data,options);
% save_figure(data,fig) % save_figure(data,fig)
...@@ -144,10 +166,12 @@ options.ST = 1; ...@@ -144,10 +166,12 @@ options.ST = 1;
options.NORMALIZED = 0; options.NORMALIZED = 0;
options.LOGSCALE = 1; options.LOGSCALE = 1;
options.FILTER = 0; %filter the 50% time-average of the spectrum from options.FILTER = 0; %filter the 50% time-average of the spectrum from
options.TAVG_2D = 1; %Show a 2D plot of the modes, 50% time averaged
options.TAVG_2D_CTR= 1; %make it contour plot
fig = show_moments_spectrum(data,options); fig = show_moments_spectrum(data,options);
end end
if 1 if 0
%% Mode evolution %% Mode evolution
[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi'); [data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
......
wk/heat_flux_convergence_analysis.m
+ 6
− 3
View file @ f93045fb
figure
hold on
if 0 if 0
%% GYAC Local low res %% GYAC Local low res
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/local_runs/nu=0.05/3x2x64x48x16/CBC/'; % resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/local_runs/nu=0.05/3x2x64x48x16/CBC/';
...@@ -22,7 +24,7 @@ if 0 ...@@ -22,7 +24,7 @@ if 0
legend('show') legend('show')
disp(['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']); disp(['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
end end
if 1 if 0
%% Manually gathered data %% Manually gathered data
QvsPJ = [... QvsPJ = [...
%vp mu Qxav Qxer kT %vp mu Qxav Qxer kT
...@@ -36,7 +38,6 @@ QvsPJ = [... ...@@ -36,7 +38,6 @@ QvsPJ = [...
11 06 27.47 05.48 6.96;... 11 06 27.47 05.48 6.96;...
13 05 16.45 04.63 6.96;... 13 05 16.45 04.63 6.96;...
]; ];
figure
errorbar(QvsPJ(:,1).*QvsPJ(:,2),QvsPJ(:,3),QvsPJ(:,4),'--s',... errorbar(QvsPJ(:,1).*QvsPJ(:,2),QvsPJ(:,3),QvsPJ(:,4),'--s',...
'DisplayName','GYAC 64x48x16 ($\nu_{DGDK}=0.001$)'); hold on 'DisplayName','GYAC 64x48x16 ($\nu_{DGDK}=0.001$)'); hold on
end end
...@@ -62,7 +63,9 @@ if 0 ...@@ -62,7 +63,9 @@ if 0
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/CBC/7x4x128x64x24/'; % resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/CBC/7x4x128x64x24/';
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/local_runs/7x4x128x64x24/CBC_L120/'; % resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/local_runs/7x4x128x64x24/CBC_L120/';
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/local_runs/7x4x128x64x24/CBC_L180/'; % resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/local_runs/7x4x128x64x24/CBC_L180/';
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/CBC/21x6x192x96x24/'; % resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/CBC/21x6x192x96x24/';
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/CBC/21x6x192x96x24/';
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/CBC/21x6x192x96x24/';
data = {}; data = {};
data = compile_results_low_mem(data,resdir,00,10); data = compile_results_low_mem(data,resdir,00,10);
% fast heat flux analysis % fast heat flux analysis
......
wk/lin_ITG_salpha.m
+ 2
− 2
View file @ f93045fb
...@@ -31,8 +31,8 @@ NA = 1; % number of kinetic species ...@@ -31,8 +31,8 @@ NA = 1; % number of kinetic species
ADIAB_E = (NA==1); % adiabatic electron model ADIAB_E = (NA==1); % adiabatic electron model
BETA = 0.0; % electron plasma beta BETA = 0.0; % electron plasma beta
%% Set up grid parameters %% Set up grid parameters
P = 12; P = 2;
J = 6;%P/2; J = 1;%P/2;
PMAX = P; % Hermite basis size PMAX = P; % Hermite basis size
JMAX = J; % Laguerre basis size JMAX = J; % Laguerre basis size
NX = 8; % real space x-gridpoints NX = 8; % real space x-gridpoints
......
wk/load_metadata_scan.m
+ 47
− 17
View file @ f93045fb
...@@ -49,10 +49,27 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0' ...@@ -49,10 +49,27 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0'
% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_10_LDGK_0.05.mat'; % datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_10_LDGK_0.05.mat';
% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_20_DGDK_0.1.mat'; % datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_20_DGDK_0.1.mat';
datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_20_DGGK_0.1.mat'; % datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_20_DGGK_0.1.mat';
% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_20_SGGK_0.1.mat'; % datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_20_SGGK_0.1.mat';
% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_10_LDGK_0.1.mat'; % datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_10_LDGK_0.1.mat';
% datafname = 'p2_linear_new/8x24_ky_0.3_P_2_J_1_kT_3_5_nu_0.001_0.1_DGGK.mat';
% datafname = 'p2_linear_new/8x24_ky_0.3_P_4_J_2_kT_2.5_4.5_nu_0.001_0.1_DGGK.mat';
% datafname = 'p2_linear_new/8x24_ky_0.3_P_6_J_3_kT_2.5_4.5_nu_0.001_0.1_DGGK.mat';
% datafname = 'p2_linear_new/8x24_ky_0.3_P_8_J_4_kT_2.5_4.5_nu_0.001_0.1_DGGK.mat';
% datafname = 'p2_linear_new/8x24_ky_0.3_P_16_J_8_kT_2.5_4.5_nu_0.001_0.1_DGGK.mat';
% datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_64_DGDK_0.001_kT_6.96.mat';
%% ky pj scan
% datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_64_DGDK_0.001_kT_6.96.mat';
% datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_48_DGDK_0.001_kT_6.96.mat';
% datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_16_DGGK_0.01_kT_5.3.mat';
% datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_16_SGGK_0.01_kT_5.3.mat';
datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_10_LDGK_0.01_kT_5.3.mat';
%% Data for the paper :
% Collisionless kT threshold
% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_30_DGDK_0.001.mat';
%% Chose if we filter gamma>0.05 %% Chose if we filter gamma>0.05
FILTERGAMMA = 1; FILTERGAMMA = 1;
...@@ -63,18 +80,20 @@ if FILTERGAMMA ...@@ -63,18 +80,20 @@ if FILTERGAMMA
d.data = d.data.*(d.data>0.025); d.data = d.data.*(d.data>0.025);
d.err = d.err.*(d.data>0.025); d.err = d.err.*(d.data>0.025);
end end
if 1 if 0
%% Pcolor of the peak %% Pcolor of the peak
figure; figure;
% [XX_,YY_] = meshgrid(d.s1,d.s2); % [XX_,YY_] = meshgrid(d.s1,d.s2);
[XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2)); [XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2));
pclr=imagesc_custom(XX_,YY_,d.data'.*(d.data>0)'); pclr=imagesc_custom(XX_,YY_,d.data'.*(d.data>0)');
% pclr=contourf(1:numel(d.s1),1:numel(d.s2),d.data'.*(d.data>0)');
% pclr=surf(1:numel(d.s1),1:numel(d.s2),d.data'.*(d.data>0)');
title(d.title); title(d.title);
xlabel(d.s1name); ylabel(d.s2name); xlabel(d.s1name); ylabel(d.s2name);
set(gca,'XTicklabel',d.s1) set(gca,'XTick',1:numel(d.s1),'XTicklabel',d.s1)
set(gca,'YTicklabel',d.s2) set(gca,'YTick',1:numel(d.s2),'YTicklabel',d.s2)
colormap(jet) colormap(jet)
% colormap(bluewhitered) colormap(bluewhitered)
clb=colorbar; clb=colorbar;
clb.Label.String = '$\gamma c_s/R$'; clb.Label.String = '$\gamma c_s/R$';
clb.Label.Interpreter = 'latex'; clb.Label.Interpreter = 'latex';
...@@ -99,9 +118,11 @@ xlabel(d.s1name); ylabel(d.dname);title(d.title); ...@@ -99,9 +118,11 @@ xlabel(d.s1name); ylabel(d.dname);title(d.title);
xlim([d.s1(1) d.s1(end)]); xlim([d.s1(1) d.s1(end)]);
colormap(colors_); colormap(colors_);
clb = colorbar; clb = colorbar;
caxis([d.s2(1)-0.5,d.s2(end)+0.5]); % caxis([d.s2(1)-0.5,d.s2(end)+0.5]);
clim([1 numel(d.s2)+1]);
clb.Ticks=linspace(d.s2(1),d.s2(end),numel(d.s2)); clb.Ticks=linspace(d.s2(1),d.s2(end),numel(d.s2));
clb.YTick=d.s2; clb.Ticks =1.5:numel(d.s2)+1.5;
clb.TickLabels=d.s2;
clb.Label.String = d.s2name; clb.Label.String = d.s2name;
clb.Label.Interpreter = 'latex'; clb.Label.Interpreter = 'latex';
clb.Label.FontSize= 18; clb.Label.FontSize= 18;
...@@ -169,29 +190,38 @@ if 0 ...@@ -169,29 +190,38 @@ if 0
figure; i_ = 0; figure; i_ = 0;
% target_ = 2.72724991618068013377e-01; % Value for nuDGDK = 1.0, kT=6.96, (40,20), Nkx=8 % target_ = 2.72724991618068013377e-01; % Value for nuDGDK = 1.0, kT=6.96, (40,20), Nkx=8
% target_ = 2.79666916212537142172e-01; % Value for nuDGDK = 0.05, kT=6.96, (40,20), Nkx=8 % target_ = 2.79666916212537142172e-01; % Value for nuDGDK = 0.05, kT=6.96, (40,20), Nkx=8
target_ = 2.73048910051283844069e-01; % Value for nuDGDK = 0.0, kT=6.96, (40,20), Nkx=8 % target_ = 2.71455e-01; % Value for nuDGDK = 0.001, kT=6.96, (40,20), Nkx=8
% target_ = 1.39427e-01; % Value for nuDGDK = 0.001, kT=5.3, (30,16), Nkx=8
target_ = 2.72510405826983714839e-01 % Value for nuDGDK = 0.001, kT=6.96, (50,25), Nkx=8
% target_ = 2.73048910051283844069e-01; % Value for nuDGDK = 0.0, kT=6.96, (40,20), Nkx=8
% target_ = 0.25*(d.data(end,end)+d.data(end-i_,end)+d.data(end,end-i_)+d.data(end-i_,end-i_)); % target_ = 0.25*(d.data(end,end)+d.data(end-i_,end)+d.data(end,end-i_)+d.data(end-i_,end-i_));
% eps_ = log(abs(target_ - d.data)/abs(target_)); % eps_ = log(abs(target_ - d.data)/abs(target_));
eps_ = max(-10,log(abs(target_ - d.data)/abs(target_))); eps_ = max(-10,log(abs(target_ - d.data)/abs(target_)));
sign_ = 1;%sign(d.data - target_); sign_ = 1;%sign(d.data - target_);
eps_ = d.data; eps_ = d.data;
for i = 1:numel(d.s1) for i = 1:numel(d.s1)
target_ = d.data(i,end); for j = 1:numel(d.s2)
eps_(i,:) = log(abs(target_ - d.data(i,1:end))); % target_ = d.data(i,end);
if target_ > 0 % target_ = d.data(i,end);
% eps_(i,:) = max(-12,log(abs(target_ - d.data(i,1:end))/abs(target_))); % target_ = d.data(end,j);
% eps_(i,:) = log(abs(target_ - d.data(i,1:end))/abs(target_)); eps_(i,j) = log(abs(target_ - d.data(i,j))/target);
% eps_(i,:) = min(100,100*abs(target_ - d.data(i,1:end))/abs(target_)); if target_ > 0
else % eps_(i,:) = max(-12,log(abs(target_ - d.data(i,1:end))/abs(target_)));
% eps_(i,:) = log(abs(target_ - d.data(i,1:end))/abs(target_));
% eps_(i,:) = min(100,100*abs(target_ - d.data(i,1:end))/abs(target_));
else
end
end end
end end
[XX_,YY_] = meshgrid(d.s1,d.s2); [XX_,YY_] = meshgrid(d.s1,d.s2);
[XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2)); [XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2));
pclr=imagesc_custom(XX_,YY_,eps_'.*(d.data>0)'.*sign_'); pclr=imagesc_custom(XX_,YY_,eps_'.*(d.data>0)'.*sign_');
% pclr=contourf(1:numel(d.s1),1:numel(d.s2),eps_'.*(d.data>0)'.*sign_',5);
% pclr=surf(1:numel(d.s1),1:numel(d.s2),eps_'.*(d.data>0)'.*sign_');
title(d.title); title(d.title);
xlabel(d.s1name); ylabel(d.s2name); xlabel(d.s1name); ylabel(d.s2name);
set(gca,'XTicklabel',d.s1) set(gca,'XTick',1:numel(d.s1),'XTicklabel',d.s1)
set(gca,'YTicklabel',d.s2) set(gca,'YTick',1:numel(d.s2),'YTicklabel',d.s2)
% colormap(jet) % colormap(jet)
colormap(bluewhitered) colormap(bluewhitered)
% caxis([-10, 0]); % caxis([-10, 0]);
......
wk/local_run.m
+ 7
− 6
View file @ f93045fb
...@@ -13,14 +13,14 @@ addpath(genpath([gyacomodir,'matlab/load'])) % Add load module ...@@ -13,14 +13,14 @@ addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
%% Set simulation parameters %% Set simulation parameters
SIMID = 'dbg'; % Name of the simulation SIMID = 'dbg'; % Name of the simulation
RUN = 1; % To run or just to load RUN = 0; % To run or just to load
default_plots_options default_plots_options
EXECNAME = 'gyacomo23_sp'; % single precision EXECNAME = 'gyacomo23_sp'; % single precision
% EXECNAME = 'gyacomo23_dp'; % double precision % EXECNAME = 'gyacomo23_dp'; % double precision
%% Set up physical parameters %% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
NU = 0.05; % Collision frequency NU = 0.001; % Collision frequency
TAU = 1.0; % e/i temperature ratio TAU = 1.0; % e/i temperature ratio
K_Ne = 0*2.22; % ele Density K_Ne = 0*2.22; % ele Density
K_Te = 0*6.96; % ele Temperature K_Te = 0*6.96; % ele Temperature
...@@ -31,13 +31,13 @@ NA = 1; % number of kinetic species ...@@ -31,13 +31,13 @@ NA = 1; % number of kinetic species
ADIAB_E = (NA==1); % adiabatic electron model ADIAB_E = (NA==1); % adiabatic electron model
BETA = 0.0; % electron plasma beta BETA = 0.0; % electron plasma beta
%% Set up grid parameters %% Set up grid parameters
P = 16; P = 60;
J = P/2; J = P/2;
DT = 1e-2; % Time step DT = 1e-2; % Time step
PMAX = P; % Hermite basis size PMAX = P; % Hermite basis size
JMAX = J; % Laguerre basis size JMAX = J; % Laguerre basis size
NX = 8; % real space x-gridpoints NX = 8; % real space x-gridpoints
NY = 12; % real space y-gridpoints NY = 2; % real space y-gridpoints
LX = 2*pi/0.1; % Size of the squared frequency domain in x direction LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
LY = 2*pi/0.1; % Size of the squared frequency domain in y direction LY = 2*pi/0.1; % Size of the squared frequency domain in y direction
NZ = 24; % number of perpendicular planes (parallel grid) NZ = 24; % number of perpendicular planes (parallel grid)
...@@ -71,7 +71,8 @@ CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:S ...@@ -71,7 +71,8 @@ CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:S
GKCO = 1; % Gyrokinetic operator GKCO = 1; % Gyrokinetic operator
ABCO = 1; % INTERSPECIES collisions ABCO = 1; % INTERSPECIES collisions
INIT_ZF = 0; % Initialize zero-field quantities INIT_ZF = 0; % Initialize zero-field quantities
HRCY_CLOS = 'truncation'; % Closure model for higher order moments % HRCY_CLOS = 'truncation'; % Closure model for higher order moments
HRCY_CLOS = 'monomial'; % Closure model for higher order moments
DMAX = -1; DMAX = -1;
NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
NMAX = 0; NMAX = 0;
......
wk/new_CBC_convergence_GENE_GYAC_comparison.m 0 → 100644
+ 181
− 0
View file @ f93045fb
figure
hold on
TRTW = [0.5 1];
ERRBAR = 1;
%%%%%%%%%%%%%%%%%%%%%%%%
%% GYAC Marconi standard res
resdirs = {...
'/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/3x2x128x64x24/kT_7.0/';
'/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/5x3x128x64x24/kT_7.0/';
'/misc/gyacomo23_outputs/paper_2_GYAC23/save/kT_scan_1e-3/7x4x128x64x24/kT_7.0/';
'/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/9x5x128x64x24/kT_7.0/';
'/misc/gyacomo23_outputs/paper_2_GYAC23/save/kT_scan_1e-3/11x6x128x64x24/kT_7.0/';
'/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/17x9x128x64x24/kT_7.0/';
'/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/31x16x128x64x24/kT_7.0/'
};
QvsPJ = zeros(numel(resdirs),4);
N = numel(resdirs);
for i = 1:N
data = {};
data = compile_results_low_mem(data,resdirs{i},00,10);
% fast heat flux analysis
T = data.Ts0D;
Qx = data.HFLUX_X;
Trange = T(end)*TRTW;
[~,it0] = min(abs(Trange(1)-T));
[~,it1] = min(abs(Trange(2)-T));
Qavg = mean(Qx(it0:it1));
Qstd = std(Qx(it0:it1));
disp(['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
QvsPJ(i,1) = data.grids.Np;
QvsPJ(i,2) = data.grids.Nj;
QvsPJ(i,3) = Qavg;
QvsPJ(i,4) = Qstd;
end
if ERRBAR
errorbar(QvsPJ(:,1).*QvsPJ(:,2),QvsPJ(:,3),QvsPJ(:,4),'LineStyle','-',...
'DisplayName','GM, $\nu_{DGDK}=0.001$','Marker','o',...
'MarkerSize',7,'LineWidth',2);
else
plot((QvsPJ(:,1).*QvsPJ(:,2)),QvsPJ(:,3),'-o','DisplayName','GM, $\nu_{DGDK}=0.001$')
end
%% GYAC Marconi nu = 0.01
resdirs = {...
'/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-2/5x3x128x64x24/kT_7.0/';
'/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-2/7x4x128x64x24/kT_7.0/';
'/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-2/9x5x128x64x24/kT_7.0/';
'/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-2/11x6x128x64x24/kT_7.0/';
'/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-2/13x7x128x64x24/kT_7.0/'
};
QvsPJ = zeros(numel(resdirs),4);
N = numel(resdirs);
for i = 1:N
data = {};
data = compile_results_low_mem(data,resdirs{i},00,10);
% fast heat flux analysis
T = data.Ts0D;
Qx = data.HFLUX_X;
Trange = T(end)*TRTW;
[~,it0] = min(abs(Trange(1)-T));
[~,it1] = min(abs(Trange(2)-T));
Qavg = mean(Qx(it0:it1));
Qstd = std(Qx(it0:it1));
disp(['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
QvsPJ(i,1) = data.grids.Np;
QvsPJ(i,2) = data.grids.Nj;
QvsPJ(i,3) = Qavg;
QvsPJ(i,4) = Qstd;
end
if ERRBAR
errorbar(QvsPJ(:,1).*QvsPJ(:,2),QvsPJ(:,3),QvsPJ(:,4),'LineStyle','none',...
'DisplayName','GM, $\nu_{DGDK}=0.01$','Marker','o',...
'MarkerSize',7,'LineWidth',2);
else
plot((QvsPJ(:,1).*QvsPJ(:,2)),QvsPJ(:,3),'-o','DisplayName','GM, $\nu_{DGDK}=0.01$')
end
%% GYAC Marconi nu = 0.05
resdirs = {...
'/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_5e-2/3x2x128x64x24/kT_7.0/';
'/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_5e-2/5x3x128x64x24/kT_7.0/';
'/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_5e-2/7x4x128x64x24/kT_7.0/';
'/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_5e-2/9x5x128x64x24/kT_7.0/';
'/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_5e-2/11x6x128x64x24/kT_7.0/';
'/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_5e-2/13x7x128x64x24/kT_7.0/'
'/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_5e-2/15x8x128x64x24/kT_7.0/'
};
QvsPJ = zeros(numel(resdirs),4);
N = numel(resdirs);
for i = 1:N
data = {};
data = compile_results_low_mem(data,resdirs{i},00,10);
% fast heat flux analysis
T = data.Ts0D;
Qx = data.HFLUX_X;
Trange = T(end)*TRTW;
[~,it0] = min(abs(Trange(1)-T));
[~,it1] = min(abs(Trange(2)-T));
Qavg = mean(Qx(it0:it1));
Qstd = std(Qx(it0:it1));
disp(['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
QvsPJ(i,1) = data.grids.Np;
QvsPJ(i,2) = data.grids.Nj;
QvsPJ(i,3) = Qavg;
QvsPJ(i,4) = Qstd;
end
if ERRBAR
errorbar(QvsPJ(:,1).*QvsPJ(:,2),QvsPJ(:,3),QvsPJ(:,4),'LineStyle','none',...
'DisplayName','GM, $\nu_{DGDK}=0.05$','Marker','o',...
'MarkerSize',7,'LineWidth',2);
else
plot((QvsPJ(:,1).*QvsPJ(:,2)),QvsPJ(:,3),'-o','DisplayName','GM, $\nu_{DGDK}=0.05$')
end
%% GENE normal box
resdirs = {...
'/misc/gene_results/kT_scan_nu0/8x4x128x64x24/kT_7.0/';
'/misc/gene_results/kT_scan_nu0/16x8x128x64x24/kT_7.0/';
'/misc/gene_results/kT_scan_nu0/30x16x128x64x24/kT_7.0/';
'/misc/gene_results/kT_scan_nu0/42x24x128x64x24/kT_7.0/'
};
QvsNv = zeros(numel(resdirs),4);
N = numel(resdirs);
for i = 1:N
% fast heat flux analysis
nrgfile = 'nrg.dat.h5';
% nrgfile = 'nrg_1.h5';
T = h5read([resdirs{i},nrgfile],'/nrgions/time');
Qx = h5read([resdirs{i},nrgfile],'/nrgions/Q_es');
Trange = T(end)*TRTW;
[~,it0] = min(abs(Trange(1)-T));
[~,it1] = min(abs(Trange(2)-T));
Qavg = mean(Qx(it0:end));
Qstd = std(Qx(it0:end));
disp(['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
namelist = read_namelist([resdirs{i},'parameters']);
QvsNv(i,1) = namelist.box.nv0;
QvsNv(i,2) = namelist.box.nw0;
QvsNv(i,3) = Qavg;
QvsNv(i,4) = Qstd;
end
if ERRBAR
errorbar((QvsNv(:,1).*QvsNv(:,2)),QvsNv(:,3),QvsNv(:,4),'k','LineStyle','-',...
'DisplayName','GENE','Marker','s',...
'MarkerSize',10,'LineWidth',2);
else
plot((QvsNv(:,1).*QvsNv(:,2)),QvsNv(:,3),'sk','DisplayName','GENE')
end
%% GENE half vbox
resdirs = {...
'/misc/gene_results/kT_scan_nu0/8x4x128x64x24_half_vbox/kT_7.0/';
'/misc/gene_results/kT_scan_nu0/16x8x128x64x24_half_vbox/kT_7.0/';
};
QvsNv = zeros(numel(resdirs),4);
N = numel(resdirs);
for i = 1:N
% fast heat flux analysis
nrgfile = 'nrg.dat.h5';
% nrgfile = 'nrg_1.h5';
T = h5read([resdirs{i},nrgfile],'/nrgions/time');
Qx = h5read([resdirs{i},nrgfile],'/nrgions/Q_es');
[~,it0] = min(abs(0.25*T(end)-T));
Qavg = mean(Qx(it0:end));
Qstd = std(Qx(it0:end));
disp(['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
namelist = read_namelist([resdirs{i},'parameters']);
QvsNv(i,1) = namelist.box.nv0;
QvsNv(i,2) = namelist.box.nw0;
QvsNv(i,3) = Qavg;
QvsNv(i,4) = Qstd;
end
if ERRBAR
errorbar((QvsNv(:,1).*QvsNv(:,2)),QvsNv(:,3),QvsNv(:,4),'k','LineStyle','none',...
'DisplayName','GENE, $L_{v_\parallel}/2,L_\mu/2$','Marker','x',...
'MarkerSize',7,'LineWidth',2);
else
plot((QvsNv(:,1).*QvsNv(:,2)),QvsNv(:,3),'xk','DisplayName','GENE, $L_{v_\parallel}/2,L_\mu/2$')
end
%% plot formatting
set(gca,'xscale','log');
xlabel('$N_{v\parallel}\times N_{\mu}$'); ylabel('$Q_x$');
grid on
\ No newline at end of file
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