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

scripts updates

parent 65a77462
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 199 additions and 161 deletions
matlab/COSOlver_matrices_analysis.m deleted 100644 → 0
+ 0
− 72
View file @ 65a77462
addpath(genpath('../matlab')) % ... add
%% Grid configuration
N = 10; % Frequency gridpoints (Nkx = N/2)
L = 120; % Size of the squared frequency domain
dk = 2*pi/L;
kmax = N/2*dk;
kx = dk*(0:N/2);
ky = dk*(0:N/2);
[ky, kx]= meshgrid(ky,kx);
KPERP = sqrt(kx.^2 + ky.^2);
kperp = reshape(KPERP,[1,numel(kx)^2]);
kperp = uniquetol(kperp,1e-14);
Nperp = numel(kperp);
COSOlver_path = '../../Documents/MoliSolver/COSOlver/';
COSOLVER.pmaxe = 10;
COSOLVER.jmaxe = 5;
COSOLVER.pmaxi = 10;
COSOLVER.jmaxi = 5;
COSOLVER.neFLR = max(5,ceil(COSOLVER.kperp^2)); % rule of thumb for sum truncation
COSOLVER.niFLR = max(5,ceil(COSOLVER.kperp^2));
COSOLVER.neFLRs = 0; % ... only for GK abel
COSOLVER.npeFLR = 0; % ... only for GK abel
COSOLVER.niFLRs = 0; % ... only for GK abel
COSOLVER.npiFLR = 0; % ... only for GK abel
COSOLVER.gk = 1;
COSOLVER.co = 3;
if 0
%% plot the kperp distribution
figure
plot(kperp)
end
%% Load the matrices
C_self_i = zeros((COSOLVER.pmaxi+1)*(COSOLVER.jmaxi+1),(COSOLVER.pmaxi+1)*(COSOLVER.jmaxi+1),Nperp);
for n_ = 1:Nperp
COSOLVER.kperp = kperp(n_);
k_string = sprintf('%0.4f',kperp(n_));
mat_file_name = ['../iCa/self_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),...
'_ESELF_',num2str(COSOLVER.co),'_ISELF_',num2str(COSOLVER.co),...
'_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),...
'_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),...
'_JE_12',...
'_NFLR_',num2str(COSOLVER.neFLR),...
'_kperp_',k_string,'.h5'];
tmp = h5read(mat_file_name,'/Caapj/Ceepj');
C_self_i(:,:,n_) = tmp;
end
%% Post processing
dC_self_i = diff(C_self_i,1,3);
gvar_dC = squeeze(sum(sum(abs(dC_self_i),2),1));
%% Plots
%% all coeffs evolution
figure
for ip_ = 1:COSOLVER.pmaxi+1
for ij_ = 1:COSOLVER.jmaxi+1
plot(kperp,squeeze(C_self_i(ip_,ij_,:)),'o'); hold on;
end
end
%% global matrix variation
figure;
kperpmid = 0.5*(kperp(2:end)+kperp(1:end-1));
plot(kperpmid,gvar_dC./diff(kperp)');
matlab/assemble_cosolver_matrices.m 0 → 100644
+ 71
− 0
View file @ f225dc1f
codir = '/home/ahoffman/cosolver/';
matname= 'gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0';
matdir = [codir,matname,'/'];
kperp = load([matdir,'kperp.log']);
Nkp = numel(kperp);
outdir = '/home/ahoffman/HeLaZ/iCa/';
outfile= [outdir,matname,'.h5'];
if(exist(outfile)>0)
system(['rm ',outdir,matname,'.h5']);
end
Nmax = numel(kperp);
for n_=1:Nmax
n_string = sprintf('%5.5d',n_-1); disp(n_string);
scandir = ['scanfiles_',n_string,'/'];
if(exist([matdir,scandir, 'ei.h5'])>0)
% Load matrices and write them in one h5 file
olddname = '/Ceipj/CeipjF'; infile = [matdir,scandir, 'ei.h5'];
newdname = ['/',n_string,olddname];
load_write_mat(infile,olddname,outfile,newdname);
olddname = '/Ceipj/CeipjT'; infile = [matdir,scandir, 'ei.h5'];
newdname = ['/',n_string,olddname];
load_write_mat(infile,olddname,outfile,newdname);
olddname = '/Ciepj/CiepjT'; infile = [matdir,scandir, 'ie.h5'];
newdname = ['/',n_string,olddname];
load_write_mat(infile,olddname,outfile,newdname);
olddname = '/Ciepj/CiepjF'; infile = [matdir,scandir, 'ie.h5'];
newdname = ['/',n_string,olddname];
load_write_mat(infile,olddname,outfile,newdname);
olddname = '/Caapj/Ceepj'; infile = [matdir,scandir,'self.h5'];
newdname = ['/',n_string,olddname];
mat_ee = load_write_mat(infile,olddname,outfile,newdname);
olddname = '/Caapj/Ciipj'; infile = [matdir,scandir,'self.h5'];
newdname = ['/',n_string,olddname];
mat_ii = load_write_mat(infile,olddname,outfile,newdname);
% to verify symmetry
verif = max(abs(imag(eig(mat_ee)))) + max(abs(imag(eig(mat_ii))));
if(verif>0) disp(['Warning, non sym matrix at ', n_string]); end;
% to verify negative EV
verif = max(real(eig(mat_ee))) + max(real(eig(mat_ii)));
if(verif>0) disp(['Warning, positive EV at ', n_string]); end;
else
break
end
end
olddname = '/Ceipj/CeipjF'; infile = [matdir,scandir, 'ei.h5'];
Pmaxe = h5readatt(infile,olddname,'Pmaxe');
Jmaxe = h5readatt(infile,olddname,'Jmaxe');
Pmaxi = h5readatt(infile,olddname,'Pmaxi');
Jmaxi = h5readatt(infile,olddname,'Jmaxi');
dims_e = [0 0];
dims_e(1)= Pmaxe; dims_e(2) = Jmaxe;
dims_i = [Pmaxi; Jmaxi];
h5create(outfile,'/dims_e',numel(dims_e));
h5write (outfile,'/dims_e',dims_e);
h5create(outfile,'/dims_i',numel(dims_i));
h5write (outfile,'/dims_i',dims_i);
%
h5create(outfile,'/coordkperp',numel(kperp(1:Nmax)));
h5write (outfile,'/coordkperp',kperp(1:Nmax)');
fid = H5F.open(outfile);
H5F.close(fid);
\ No newline at end of file
matlab/extract_fig_data.m
+ 9
− 9
View file @ f225dc1f
......@@ -7,20 +7,20 @@ for i = 1:numel(dataObjs)
X_ = [X_ dataObjs(i).XData];
Y_ = [Y_ dataObjs(i).YData];
end
n0 = 100;
n0 = 1;
figure;
mvm = @(x) movmean(x,1);
shift = X_(n0);
% shift = 0;
% plot(X_(n0:end),Y_(n0:end));
plot(mvm(X_(n0:end)),mvm(Y_(n0:end))./max(Y_(200:end))); hold on;
plot(mvm(X_(n0:end)-shift),mvm(Y_(n0:end))); hold on;
n1 = n0+1; n2 = min(n1 + 50000,numel(Y_));
avg_ = mean(Y_(n1:n2));
std_ = std(Y_(n1:n2));
title(['avg = ',num2str(avg_),' std = ', num2str(std_)])
plot([X_(n1),X_(n2)],[1 1]*avg_,'--k')
%
% n1 = n0+1; n2 = min(n1 + 50000,numel(Y_));
% avg_ = mean(Y_(n1:n2));
% std_ = std(Y_(n1:n2));
% title(['avg = ',num2str(avg_),' std = ', num2str(std_)])
% plot([X_(n1),X_(n2)],[1 1]*avg_,'--k')
% ylim([0,avg_*3]);
sum(Y_(2:end)./X_(2:end).^(3/2))
......@@ -35,4 +35,4 @@ if 0
end
xlim([-3,3]); ylim([0,4.5]);
\ No newline at end of file
% xlim([-3,3]); ylim([0,4.5]);
\ No newline at end of file
matlab/load/load_params.m
+ 2
− 2
View file @ f225dc1f
......@@ -20,8 +20,8 @@ DATA.L = h5readatt(filename,'/data/input','Lx');
DATA.CLOS = h5readatt(filename,'/data/input','CLOS');
DATA.DT_SIM = h5readatt(filename,'/data/input','dt');
DATA.MU = h5readatt(filename,'/data/input','mu');
DATA.MUx = -99;%h5readatt(filename,'/data/input','mu_x');
DATA.MUy = -99;%h5readatt(filename,'/data/input','mu_y');
DATA.MUx = h5readatt(filename,'/data/input','mu_x');
DATA.MUy = h5readatt(filename,'/data/input','mu_y');
% MU = str2num(filename(end-18:end-14)); %bad...
DATA.W_GAMMA = h5readatt(filename,'/data/input','write_gamma') == 'y';
DATA.W_PHI = h5readatt(filename,'/data/input','write_phi') == 'y';
......
matlab/load_write_mat.m 0 → 100644
+ 15
− 0
View file @ f225dc1f
function mat_ = load_write_mat(infile,olddname,outfile,newdname)
mat_ = h5read (infile,olddname);
h5create (outfile,newdname,size(mat_));
h5write (outfile,newdname,mat_);
if(~strcmp(olddname(end-3:end-2),'ii'))
neFLR = h5readatt(infile,olddname,'neFLR');
h5writeatt(outfile,newdname,'neFLR',neFLR);
end
if(~strcmp(olddname(end-3:end-2),'ee'))
niFLR = h5readatt(infile,olddname,'niFLR');
h5writeatt(outfile,newdname,'niFLR',niFLR);
end
end
matlab/plot/plot_cosol_mat.m
+ 11
− 10
View file @ f225dc1f
......@@ -116,11 +116,12 @@ if 0
%%
mfns = {'/home/ahoffman/HeLaZ/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',...
'/home/ahoffman/cosolver/gk.coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
'/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
'/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_75_kpm_6.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0.h5',...
};
CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'FC 4 2 NFLR 16'};
CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'Hacked SG A', 'Hacked SG B'};
figure
for j_ = 1:numel(mfns)
mat_file_name = mfns{j_}; disp(mat_file_name);
......@@ -150,14 +151,14 @@ end
%% Van Kampen plot
if 0
%%
kperp= 3.9;
kperp= 0.1;
mfns = {'/home/ahoffman/HeLaZ/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',...
'/home/ahoffman/cosolver/gk.coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
'/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
'/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_75_kpm_6.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
'/home/ahoffman/HeLaZ/wk/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5',...
};
CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'FC 4 2 NFLR 12 extended'};
CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'Hacked SG'};
grow = {};
puls = {};
for j_ = 1:numel(mfns)
......@@ -177,4 +178,4 @@ for j_ = 1:numel(mfns)
end
legend('show'); grid on;
xlabel('$\omega$ from Eig(iCa)'); ylabel('$\gamma$ from Eig(iCa)')
end
end
\ No newline at end of file
matlab/plot/plot_param_evol.m
+ 8
− 0
View file @ f225dc1f
......@@ -9,12 +9,20 @@ MUy_EVOL = data.MUy_EVOL;
K_N_EVOL = data.K_N_EVOL;
L_EVOL = data.L_EVOL;
CO_LIST = {};
for i_ = 1:numel(CO_EVOL)/6
CO_LIST{i_} = CO_EVOL(6*(i_-1)+1:6*i_);
end
fig = figure;
hold on; set(gcf, 'Position', [100, 100, 400, 1500])
subplot(4,1,1);
title('Parameter evolution'); hold on;
yyaxis left
plot(TJOB_SE,NU_EVOL,'DisplayName','$\nu$');
for i_ = 1:numel(CO_EVOL)/12
k_ = 2*(i_-1)+1;
text(double(TJOB_SE(k_)),NU_EVOL(k_)*2.0,CO_LIST{k_});
end
yyaxis right
plot(TJOB_SE,DT_EVOL,'DisplayName','dt');
xlim([TJOB_SE(1) TJOB_SE(end)]);legend('show');grid on;
......
matlab/plot/plot_radial_transport_and_spacetime.m
+ 17
− 15
View file @ f225dc1f
......@@ -45,21 +45,6 @@ mvm = @(x) movmean(x,Nmvm);
grid on; set(gca,'xticklabel',[]); ylabel('$\Gamma_x$')
ylim([0,5*abs(gamma_infty_avg)]); xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
title([DATA.param_title,', $\Gamma^{\infty} = $',num2str(gamma_infty_avg),'$\pm$',num2str(gamma_infty_std)]);
% plot
subplot(312)
it0 = 1; itend = Ns3D;
trange = it0:itend;
plt1 = @(x) x;%-x(1);
plt2 = @(x) x./max((x(its3D:ite3D)));
% plty = @(x) x(500:end)./max(squeeze(x(500:end)));
yyaxis left
plot(plt1(DATA.Ts3D(trange)),plt2(E_Zmode_SK(trange)),'DisplayName','$k_{zf}^2|\phi_{kzf}|^2$');
ylim([-0.1, 1.5]); ylabel('$E_{Z}$')
yyaxis right
plot(plt1(DATA.Ts3D(trange)),plt2(E_NZmode_SK(trange)),'DisplayName','$(1+k^2)|\phi^k|^2$');
xlim([DATA.Ts3D(it0), DATA.Ts3D(itend)]);
ylim([-0.1, 1.5]); ylabel('$E_{NZ}$')
xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]);
%% radial shear radial profile
% computation
Ns3D = numel(DATA.Ts3D);
......@@ -108,4 +93,21 @@ mvm = @(x) movmean(x,Nmvm);
subplot(311)
plot(DATA.Ts3D,squeeze(mean(plt(f2plot),1)));
end
%% Zonal vs NZonal energies
subplot(312)
it0 = 1; itend = Ns3D;
trange = it0:itend;
plt1 = @(x) x;%-x(1);
plt2 = @(x) x./max((x(its3D:ite3D)));
toplot = sum(squeeze(plt(f2plot)).^2,1); %ST from before
% plty = @(x) x(500:end)./max(squeeze(x(500:end)));
yyaxis left
% plot(plt1(DATA.Ts3D(trange)),plt2(E_Zmode_SK(trange)),'DisplayName','$k_{zf}^2|\phi_{kzf}|^2$');
plot(plt1(DATA.Ts3D(trange)),plt2(toplot(trange)),'DisplayName','Sum $A^2$');
ylim([-0.1, 1.5]); ylabel('$E_{Z}$')
yyaxis right
plot(plt1(DATA.Ts3D(trange)),plt2(E_NZmode_SK(trange)),'DisplayName','$(1+k^2)|\phi^k|^2$');
xlim([DATA.Ts3D(it0), DATA.Ts3D(itend)]);
ylim([-0.1, 1.5]); ylabel('$E_{NZ}$')
xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]);
end
\ No newline at end of file
matlab/plot/statistical_transport_averaging.m
+ 4
− 4
View file @ f225dc1f
......@@ -5,9 +5,9 @@ Trange = options.T;
[~,it1] = min(abs(Trange(end)-data.Ts0D));
gamma = data.PGAMMA_RI(it0:it1)*scale;
dt_const = numel(unique(round(diff(data.Ts0D(it0:it1))*100)))==1;
if ~dt_const
disp('DT not const on given interval');
else
% if ~dt_const
% disp('DT not const on given interval');
% else
Ntot = (it1-it0)+1;
......@@ -29,6 +29,6 @@ else
% subplot(212)
% errorbar(N_seg,transp_seg_avg,transp_seg_std);
% xlabel('Averaging #points'); ylabel('transport value');
end
% end
end
matlab/setup.m
+ 4
− 1
View file @ f225dc1f
......@@ -51,10 +51,13 @@ COLL.mat_file = '''null''';
switch CO
case 'SG'
COLL.mat_file = '''../../../iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5''';
% COLL.mat_file = '''../../../iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5''';
% COLL.mat_file = '''../../../iCa/gk.hacked_sugama_P_4_J_2_N_75_kpm_5.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''';
% COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_75_kpm_6.0.h5''';
COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5''';
end
% Time integration and intialization parameters
TIME_INTEGRATION.numerical_scheme = '''RK4''';
......
wk/analysis_3D.m
+ 17
− 17
View file @ f225dc1f
......@@ -10,7 +10,7 @@ system(['mkdir -p ',MISCDIR]);
CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
system(CMD);
% Load outputs from jobnummin up to jobnummax
JOBNUMMIN = 00; JOBNUMMAX = 10;
JOBNUMMIN = 00; JOBNUMMAX = 20;
data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
......@@ -21,7 +21,7 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
TAVG_0 = 500; TAVG_1 = 10000; % Averaging times duration
TAVG_0 = 1000; TAVG_1 = 11000; % Averaging times duration
% chose your field to plot in spacetime diag (uzf,szf,Gx)
fig = plot_radial_transport_and_spacetime(data,TAVG_0,TAVG_1,'v_y',1);
save_figure(data,fig)
......@@ -29,7 +29,7 @@ end
if 0
%% statistical transport averaging
options.T = [5000 6000];
options.T = [1500 2500];
fig = statistical_transport_averaging(data,options);
end
......@@ -39,16 +39,16 @@ if 0
options.INTERP = 1;
options.POLARPLOT = 0;
options.NAME = '\phi';
% options.NAME = 'v_y';
% options.NAME = 'v_x';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
options.PLAN = 'xy';
% options.PLAN = 'xy';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 1;
% options.TIME = data.Ts5D;
options.TIME = 0:3:1000;
options.TIME = 900:10:2000;
data.a = data.EPS * 2000;
create_film(data,options,'.gif')
end
......@@ -56,19 +56,19 @@ end
if 0
%% 2D snapshots
% Options
options.INTERP = 0;
options.INTERP = 1;
options.POLARPLOT = 0;
options.AXISEQUAL = 1;
% options.NAME = '\phi';
% options.NAME = 'v_y';
options.NAME = 'n_e^{NZ}';
options.NAME = 'v_y';
% options.NAME = 'n_e^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
options.PLAN = 'xy';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
options.TIME = [550:20:800];
options.TIME = [000:200:1000];
data.a = data.EPS * 2000;
fig = photomaton(data,options);
save_figure(data,fig)
......@@ -92,7 +92,7 @@ options.XPERP = linspace( 0,6,64);
% options.SPAR = vp';
% options.XPERP = mu';
options.Z = 1;
options.T = 4000;
options.T = 5500;
options.CTR = 1;
options.ONED = 0;
fig = plot_fa(data,options);
......@@ -104,7 +104,7 @@ if 0
% options.TIME = 'avg';
options.TIME = 1000:4000;
options.P2J = 1;
options.ST = 1;
options.ST = 0;
options.NORMALIZED = 0;
fig = show_moments_spectrum(data,options);
save_figure(data,fig)
......@@ -112,15 +112,15 @@ end
if 0
%% 1D spectrum
options.TIME = 5000:10:5200;
options.TIME = 5000:10:5050;
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 = 1;
options.OK = 0;
options.COMPXY = 'max';
options.COMPXY = 'sum';
options.COMPT = 'avg';
options.PLOT = 'semilogy';
fig = spectrum_1D(data,options);
......@@ -147,7 +147,7 @@ if 0
%% Mode evolution
options.NORMALIZED = 0;
options.K2PLOT = 0.01:0.01:1.0;
options.TIME = 5000:1:5050;
options.TIME = 4900:1:5100;
options.NMA = 1;
options.NMODES = 20;
fig = mode_growth_meter(data,options);
......
wk/analysis_header.m
+ 11
− 2
View file @ f225dc1f
......@@ -120,15 +120,18 @@ outfile ='';
% outfile = 'nu_0.1_transport_scan/colless_kn_1.7_to_2.0';
% outfile = 'nu_0.1_transport_scan/colless_kn_2.1_to_2.5';
% outfile = 'nu_0.1_transport_scan/LB_kn_2.0';
% outfile = 'nu_0.1_transport_scan/DG_kn_1.8_to_2.1';
% outfile = 'nu_0.1_transport_scan/DG_kn_2.2_to_2.5';
% outfile = 'nu_0.1_transport_scan/DG_conv_kN_1.9';
% outfile = 'nu_0.1_transport_scan/SG_kn_1.7_to_2.0';
% outfile = 'nu_0.1_transport_scan/SG_10x5_conv_test';
% outfile = 'nu_0.1_transport_scan/SG_kn_2.2_to_2.5';
% outfile = 'nu_0.1_transport_scan/LD_kn_2.0_to_2.5';
% outfile = 'nu_0.1_transport_scan/LD_kn_1.7_to_2.5';
% outfile = 'nu_0.1_transport_scan/LR_kn_1.7_to_2.0';
% outfile = 'nu_0.1_transport_scan/LR_kn_2.1_to_2.5';
......@@ -139,10 +142,16 @@ outfile ='';
% outfile = 'nu_0.1_transport_scan/colless_kn_1.6_HD';
% outfile = 'nu_0.1_transport_scan/large_box_kN_2.1_nu_0.1';
% outfile = 'nu_0.1_transport_scan/large_box_kN_2.0_nu_0.1';
% outfile = 'predator_prey/SG_Kn_1.7_nu_0.01';
outfile = 'predator_prey_nu_scan/DG_Kn_1.7_nu_0.01';
outfile = 'ZF_damping_DK/SG_4x2_nu_0.1';
% outfile = 'ZF_damping_linear_nu_0_20x10_kn_1.6_GK/LR_4x2_nu_0.1';
% outfile = 'ZF_damping_nu_0_20x10_kn_1.6_GK/HSG_4x2_nu_0.1';
% outfile = 'ZF_damping_nu_0_5x3_kn_2.5_GK/LR_4x2_nu_0.1';
% outfile = 'hacked_sugama/hacked_B_kn_1.6_200x32_L_120x60_nu_0.1';
% outfile = 'shearless_cyclone/200x32x24_5x4_Lx_120_Ly_60_q0_1.4_e_0.18_kN_2.22_kT_6.9_nu_0';
% else% Marconi results
% outfile ='';
% outfile ='';
......
wk/benchmark_HeLaZ_gene_transport_zpinch.m
+ 16
− 18
View file @ f225dc1f
......@@ -83,20 +83,14 @@ fig = figure; set(gcf,'Position',[250 250 600 300]);
% SGGK_200x32x11x06 = [1.5e-2 1.3e+0 0.0e+0];
% DGGK_200x32x11x06 = [2.4e-3 3.8e-1 0.0e+0];
% upward scan
KN = [ 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5];
LDGK_200x32x05x03 = [2.6e-2 2.6e-1 6.3e-1 1.3e+0 2.3e+0 3.6e+0 6.3e+0 9.4e+0 1.5e+1 2.3e+1 3.8e+1];
LRGK_200x32x05x03 = [0.0e+0 4.0e-1 7.2e-1 1.2e+0 1.7e+0 2.3e+0 3.5e+0 6.0e+0 4.8e+0 6.0e+0 2.9e+1];
SGGK_200x32x05x03 = [0.0e-9 1.5e-2 1.0e-1 3.8e-1 8.4e-1 1.4e+0 2.5e+0 4.2e+0 8.0e+0 1.6e+1 3.5e+1];
DGGK_200x32x05x03 = [2.0e-4 3.6e-3 1.3e-2 4.0e-2 1.2e-1 1.8e-1 2.0e-1 2.2e+0 6.6e-1 1.3e+1 3.2e+1];
NOCO_200x32x05x03 = [0.0e+0 1.2e-2 2.6e-2 9.0e-2 2.3e-1 3.1e-1 7.0e-1 3.9e+0 6.9e+0 1.2e+1 2.4e+1];
% downward scan
% KN = [KN 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5];
% LDGK_200x32x05x03 = [LDGK_200x32x05x03 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
% LRGK_200x32x05x03 = [LRGK_200x32x05x03 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
% SGGK_200x32x05x03 = [SGGK_200x32x05x03 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
% DGGK_200x32x05x03 = [DGGK_200x32x05x03 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
% NOCO_200x32x05x03 = [NOCO_200x32x05x03 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 3.9e+0 6.9e+0 1.2e+1 0.0e+0];
KN = [ 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5];
LDGK_200x32x05x03 = [2.6e-2 2.6e-1 6.2e-1 1.2e+0 2.3e+0 3.6e+0 7.4e+0 1.1e+1 1.5e+1 2.3e+1 3.8e+1]; % validate with two box sizes (120 and 240)
LRGK_200x32x05x03 = [0.0e+0 4.0e-1 7.2e-1 1.2e+0 1.7e+0 2.3e+0 4.0e+0 6.0e+0 1.1e+1 1.9e+1 2.9e+1];
SGGK_200x32x05x03 = [0.0e-9 1.5e-2 9.5e-2 3.5e-1 8.0e-1 1.4e+0 2.5e+0 4.2e+0 8.0e+0 1.6e+1 3.5e+1];
hacked_sugama = [0.0e+0 7.0e-1 0.0e+0 2.5e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
DGGK_200x32x05x03 = [2.0e-4 5.0e-3 9.2e-3 3.5e-2 1.3e-1 3.5e-1 6.0e-1 6.0e+0 1.1e+1 1.3e+1 3.2e+1]; % validate with two box sizes (120 and 240)
LBGK_200x32x05x03 = [0.0e+0 0.0e-1 0.0e+0 0.0e+0 0.0e+0 2.8e-1 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0]; % validate with two box sizes (120 and 240)
NOCO_200x32x05x03 = [0.0e+0 1.2e-2 4.0e-2 9.0e-2 2.3e-1 3.1e-1 7.0e-1 3.9e+0 6.9e+0 1.2e+1 2.4e+1];
% X = Ginf_GENE;
......@@ -109,6 +103,7 @@ semilogy(X,plt(DGGK_200x32x05x03),'v','DisplayName','Dougherty','MarkerSize',msi
semilogy(X,plt(SGGK_200x32x05x03),'s','DisplayName', 'Sugama','MarkerSize',msize,'Color',clr_(1,:)); hold on;
semilogy(X,plt(LDGK_200x32x05x03),'d','DisplayName', 'Landau','MarkerSize',msize,'Color',clr_(5,:)); hold on;
semilogy(X,plt(LRGK_200x32x05x03),'^','DisplayName', 'Lorentz','MarkerSize',msize,'Color',clr_(3,:)); hold on;
semilogy(X,plt(hacked_sugama),'s','DisplayName', 'hacked Sugama','MarkerSize',msize,'Color',clr_(6,:)); hold on;
semilogy(X,plt(NOCO_200x32x05x03),'*k','DisplayName', '$\nu = 0$','MarkerSize',msize); hold on;
X = [ 1.6 2.0 2.5];
......@@ -122,11 +117,14 @@ xlim([ 1.55 2.55]);
% saveas(fig,'/home/ahoffman/Dropbox/Applications/Overleaf/Paper 1/figures/coll_transport_benchmark.eps')
%% Burst study Kn = 1.7
Kn = 1.7;
NU = [0.0e+0 1.0e-2 2.0e-2 3.0e-2 4.0e-2 5.0e-2 6.0e-2 7.0e-2 8.0e-2 9.0e-2 1.0e-1];
SGGK_200x32x05x03 = [1.0e-2 2.4e-2 2.8e-2 3.5e-2 4.5e-2 5.5e-2 6.5e-2 7.5e-2 8.3e-2 8.8e-2 1.0e-1];
SGGK_200x32x10x05 = [1.0e-2 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
LRGK_200x32x05x03 = [0.0e+0 2.9e-2 0.0e+0 0.0e-2 0.0e-2 0.0e+0 0.0e+0];
DGGK_200x32x05x03 = [1.0e-2 7.0e-3 7.3e-3 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
LDGK_200x32x05x03 = [1.0e-2 8.5e-2 1.4e-1 2.0e-1 0.0e-2 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
% NOCO_200x32x05x03 = [0.0e+0 0.0e+0 0.0e+0 0.0e+0];
figure
plot(NU,SGGK_200x32x05x03,'x'); hold on
grid on; xlabel('$\nu$'); ylabel('$\Gamma^\infty_x$');
plot(NU,SGGK_200x32x05x03/Kn,'s','Color',clr_(1,:)); hold on
plot(NU,DGGK_200x32x05x03/Kn,'v','Color',clr_(2,:)); hold on
plot(NU,LDGK_200x32x05x03/Kn,'d','Color',clr_(5,:)); hold on
grid on; xlabel('$\nu$'); ylabel('$\Gamma_x^\infty/\kappa_N$');
wk/linear_1D_entropy_mode.m
+ 9
− 9
View file @ f225dc1f
......@@ -8,12 +8,12 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%% PHYSICAL PARAMETERS
NU = 0.1; % Collision frequency
TAU = 1.0; % e/i temperature ratio
K_N = 1.7; % Density gradient drive
K_N = 1.6; % Density gradient drive
K_T = 0.25*K_N; % Temperature '''
K_E = 0.0; % Electrostat '''
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
%% GRID PARAMETERS
NX = 64; % real space x-gridpoints
NX = 40; % real space x-gridpoints
NY = 1; % '' y-gridpoints
LX = 120; % Size of the squared frequency domain
LY = 1; % Size of the squared frequency domain
......@@ -23,8 +23,8 @@ SHEAR = 0.0; % magnetic shear
EPS = 0.0; % inverse aspect ratio
SG = 1; % Staggered z grids option
%% TIME PARMETERS
TMAX = 100; % Maximal time unit
DT = 2e-2; % Time step
TMAX = 200; % Maximal time unit
DT = 1e-2; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
SPS3D = 1; % Sampling per time unit for 2D arrays
......@@ -37,11 +37,11 @@ LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
KIN_E = 1;
% Collision operator
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
CO = 'DG';
CO = 'LD';
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))
CLOS = 0; % Closure model (0: =0 truncation, 1: gyrofluid 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= 'mom00'; % Start simulation with a noisy mom00/phi/allmom
......@@ -72,11 +72,11 @@ CURVB = 1.0;
if 1
% Parameter scan over PJ
% PA = [4 6 8 10];
% JA = [2 3 4 5];
PA = [4];
JA = [2];
Nparam = numel(PA);
% Parameter scan over KN
PA = [4]; JA = [2];
% PA = [4]; JA = [2];
% PMAXE = PA(1); PMAXI = PA(1);
% JMAXE = JA(1); JMAXI = JA(1);
% KNA = 1.5:0.05:2.5;
......
wk/mode_growth_meter.m
+ 1
− 1
View file @ f225dc1f
......@@ -69,7 +69,7 @@ end
subplot(2,3,1+3*(i-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;
% pclr = imagesc(t,fftshift(k,1),abs(plt(fftshift(DATA.PHI,MODES_SELECTOR),1:numel(k))));
set(gca,'YDir','normal')
% xlim([t(1) t(end)]); %ylim([1e-5 1])
xlabel(['$',kstr,'\rho_s$']); ylabel('$t c_s /\rho_s$');
title(MODESTR)
......
wk/quick_gene_load.m
+ 4
− 1
View file @ f225dc1f
......@@ -24,6 +24,9 @@ path = '/home/ahoffman/gene/linear_zpinch_results/';
% fname ='GENE_LIN_Kn_1.5_KT_0.325_nuSG_0.0047_64x32.txt';
% fname ='GENE_LIN_Kn_1.4_KT_0.35_nuSG_0.0047_64x32.txt';
% fname ='GENE_LIN_Kn_1.6_KT_0.4_nuSG_0.0047_64x32.txt';
% fname ='GENE_LIN_Kn_1.6_KT_0.4_nuSGDK_0.0047_64x32.txt';
fname ='GENE_LIN_Kn_1.6_KT_0.4_nuLDDK_0.0047_64x32.txt';
% fname ='GENE_LIN_Kn_1.8_KT_0.4_nuLDDK_0.0047_64x32.txt';
% fname ='GENE_LIN_Kn_1.8_KT_0.45_nu_0_32x16.txt';
% fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSG_0.0235_64x32.txt';
% fname ='GENE_LIN_Kn_1.9_KT_0.475_nuSG_0.047_64x32.txt';
......@@ -36,7 +39,7 @@ path = '/home/ahoffman/gene/linear_zpinch_results/';
% fname ='GENE_LIN_Kn_2.0_KT_0.5_nu_0_32x16.txt';
% fname ='GENE_LIN_Kn_2.0_KT_0.5_nuSGDK_0.0235_32x16.txt';
% fname ='GENE_LIN_Kn_1.6_KT_0.4_nu_0_32x16.txt';
fname ='GENE_LIN_Kn_2.5_KT_0.625_nu_0_32x16.txt';
% fname ='GENE_LIN_Kn_2.5_KT_0.625_nu_0_32x16.txt';
data_ = load([path,fname]);
figure
......
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